Evolution and Plasticity of the Transcriptome Under Temperature Fluctuations in the Fungal Plant Pathogen Zymoseptoria tritici

Abstract

Most species live in a variable environment in nature. Yet understanding the evolutionary processes underlying molecular adaptation to fluctuations remains a challenge. In this study we investigate the transcriptome of the fungal wheat pathogen Zymoseptoria tritici after experimental evolution under stable or fluctuating temperature, by comparing ancestral and evolved populations simultaneously. We found that temperature regimes could have a large and pervasive effect on the transcriptome evolution, with as much as 38% of the genes being differentially expressed between selection regimes. Although evolved lineages showed different changes of gene expression based on ancestral genotypes, we identified a set of genes responding specifically to fluctuation. We found that transcriptome evolution in fluctuating conditions was repeatable between parallel lineages initiated from the same genotype for about 60% of the differentially expressed genes. Further, we detected several hotspots of significantly differentially expressed genes in the genome, in regions known to be enriched in repetitive elements, including accessory chromosomes. Our findings also evidenced gene expression evolution toward a gain of robustness (loss of phenotypic plasticity) associated with the fluctuating regime, suggesting robustness is adaptive in changing environment. This work provides valuable insight into the role of transcriptional rewiring for rapid adaptation to abiotic changes in filamentous plant pathogens.

Introduction

One remaining challenge in evolutionary genetics is to understand the role of regulatory variation in adaptation. In particular, in the context of rapid climate change with growing evidence that species experience more abiotic fluctuations (Fischer and Knutti, 2015), understanding the role of regulatory variants in adaptation to changing environments is increasingly attracting interest. Since early arguments in favor of regulatory variation explaining phenotypic differences between organisms (Britten and Davidson, 1971; King and Wilson, 1975), over the past two decades or so empirical studies have found abundant variation of gene expression within and between species, e.g., in Primates (Enard et al., 2002) including Humans (Rockman and Wray, 2002), in teleost fish (Oleksiak et al., 2002; Whitehead and Crawford, 2006), Drosophila (Wittkopp et al., 2004), and Saccharomyces (Fay et al., 2004). Heritable variation of gene expression within and between populations is widespread and much of which is thought to be adaptive (Brem and Kruglyak, 2002; Leder et al., 2015; Noormohammad et al., 2017; Glaser-Schmitt and Parsch, 2018).

Extensive body of work has shown that environmental changes cause rapid gene expression response [e.g., Caenorhabditis elegans (Li et al., 2006); Escherichia coli (Feugeas et al., 2016); Drosophila melanogaster (Zhou et al., 2012)]. Among all of the significant environmental factors, temperature is an important one. Many studies have reported large and pervasive effect of temperature on the transcriptome [e.g., using constant temperature stress (Podrabsky, 2004; Zhou et al., 2012; Chen et al., 2015; Jovic et al., 2017) or fluctuating temperature stress (Sørensen et al., 2016)].

Phenotypic plasticity is the ability for a genotype to produce different phenotypes in different environments (Pigliucci et al., 2006). Likewise a change in gene expression in response to temperature variation is considered as plasticity. Plasticity may play an important role for species adaptation to new environments by offering phenotypic variation prior to mutation accumulation. However our current understanding of the role of phenotypic plasticity in adaptation remains limited and it is still a matter of debate whether plasticity speeds up or impedes adaptation (Ghalambor et al., 2007; Levis and Pfennig, 2016; Fox et al., 2019). This question remains open as illustrated by recent work published in landmark papers either demonstrating adaptation to new environment by mean of plasticity – adaptive plasticity – (Ghalambor et al., 2015; Kenkel and Matz, 2016; Corl et al., 2018), or on the contrary supporting that plasticity does not facilitate adaptation (Ho and Zhang, 2018). In contrast to empirical studies insights from theoretical studies favor adaptive evolution of plasticity in fluctuating environment (Lande, 2009).

One powerful approach to study how gene expression and gene expression plasticity evolve in response to environmental variation is experimental evolution, which consists of laboratory-controlled evolution that lasts over several generations (Garland and Rose, 2009). The potential of experimental evolution in understanding adaptation to environmental fluctuations has been proven since the long term experimental evolutions in Escherichia coli under changing temperatures (Bennett et al., 1992; Leroi et al., 1994) or variable pH (Hughes et al., 2007). More recently, experimental evolution was also used to measure the effect of environmental variation on gene expression evolution [e.g., in D. melanogaster exposed to diet fluctuations (Huang and Agrawal, 2016; Zandveld et al., 2017); in Saccharomyces cerevisiae alternatively exposed to salt and oxidative stresses (Dhar et al., 2013); in Caenorhabditis remanei exposed to heat and oxidative stress (Sikkink et al., 2019)].

In contrast, there are fewer studies in genera deprived of historical model species, such as for filamentous fungi (Fisher and Lang, 2016). To fill this gap, we investigated the transcriptome evolution of the wheat fungal pathogen Zymoseptoria tritici. Toward this goal we analyzed the effect of temperature fluctuations on gene expression evolution and gene expression plasticity. Like other microorganisms, Z. tritici is an interesting model for experimental evolution due to its short generation time, small genome and the ease to maintain it in the laboratory. Z. tritici is a filamentous fungus, an ascomycete of the Mycosphaerellaceae family, and is the main causal agent of the Septoria Tritici Blotch disease of wheat (O’Driscoll et al., 2014; Fones and Gurr, 2015). Z. tritici is a haploid species that multiply asexually in vitro by budding (Steinberg, 2015). Nowadays this species becomes a good fungal model with growing interest to study its genome evolution since the publication of a complete reference genome (Goodwin et al., 2011). The genome of the reference strain has 21 chromosomes, including 13 gene-rich core chromosomes (CCs) and 8 accessory chromosomes (ACs) carrying less genes and more repetitive elements. While much attention has been paid toward genes underlying the ability of the pathogen to overcome the immune system of its host (through QTL linkage mapping (Meile et al., 2018; Stewart et al., 2018), or genome wide association studies (Hartmann et al., 2017; Zhong et al., 2017), we still have a poor understanding of the potential ability of Z. tritici to adapt to abiotic changes. Very few studies have found contrasted temperature sensitivity among natural population samples (Zhan and McDonald, 2011) or QTLs for thermal adaptation (Lendenmann et al., 2016). Transcriptome studies in the fungal pathogen Z. tritici are fairly recent and aimed for the most part at characterizing the waves of up- and down-regulated genes in association with symptom development inside the host plant (Kellner et al., 2014; Rudd et al., 2015; Palma-Guerrero et al., 2017).

Beside some evidence that there is a substantial amount of genetic variation among natural populations of Z. tritici, the contribution of regulatory variation to adaptation for this fungal species is virtually unknown. Here, we present the results of an experimental evolution to address several fundamental issues about fungal adaptation to thermal fluctuations: how does the transcriptome evolve in fluctuating temperature conditions, and to what extent is evolution repeatable between independent replicates and between genotypes? In particular, we test the assumption that phenotypic plasticity should evolve as an adaptation to a fluctuating environment.

We compared the transcriptome of evolved lineages maintained under stable or fluctuating temperature regimes, making use of 40 RNA samples obtained at two temperature treatments. We identified a pervasive effect of the selection regime on the gene expression level with a strong influence of the ancestral genotype. Results showed a few hotspots of significantly differentially expressed genes in the genome, in regions known to be enriched in repetitive elements including accessory chromosomes. Our results also showed a gain of gene expression robustness associated with the fluctuating regime as opposed to stable regime, suggesting robustness is adaptive in changing environment.

Materials and Methods

Biological Material

The two laboratory clones of Z. tritici MGGP01 and MGGP44 were collected on wheat in 2010 from the same population sample in the south of France at the location of Auzeville-Tolosane (43°53′ N – 1°48′W). For each clone, single colony was isolated and grown in the laboratory using Potato Dextrose Broth (PDB) medium at 17°C, 140 rpm and 70% humidity, for a few weeks. Strains were thus acclimated to our laboratory conditions prior to the experimental evolution. Strains were then multiplied for storage in PDB medium containing 30% glycerol at −80°C prior to the experimental evolution.

Experimental Evolution Design

Experimental evolution was conducted for two clonal strains of Z. tritici (MGGP01 and MGGP44), by multiplying large asexual populations in the lab in vitro, using 500 mL of PDB at 140 rpm in a shaking incubator (Infors HT, Inc.). Every week 10⁷–10⁹ cells were transferred into fresh medium by pipetting 20 mL of cell suspension (Figure 1A). Three different selection regimes were used: two stable temperatures (17°C and 23°C) and a third condition with fluctuating temperature between 17 and 23°C rapidly every 52–64 h (Figure 1B). These two temperatures were non-stressful and allow sufficient cell growth by budding, as opposed to hyphal growth and mycelium formation induced by high temperature, nutrient-poor medium or oxidative stress (Francisco et al., 2019). In these laboratory conditions, we estimated that 1.5 mitotic division per day occurs for both genetic backgrounds, at 17 and 23°C. We thus estimate that 48 weeks of evolution correspond approximately to 500 generations. For each clone and abiotic condition three lineages were evolved. At the end of the experimental evolution after 48 weeks, large amounts of cell suspension were collected to provide liquid stocks in 30% glycerol at −80°C.

FIGURE 1

RNA Isolation and Sequencing

RNA samples were isolated from the two ancestor clones and 8 evolved lineages (Figure 1B). As many samples as possible were sequenced under our budget constraints: 2 fluctuating, 1 stable at 17°C and 1 stable at 23°C for each ancestral genotype (MGGP01 and MGGP44). The selection of the strains within each regime was random. This design allowed to analyze (1) the transcriptome evolution under fluctuation using stable lines as controls for high and low temperatures (2) gene expression plasticity evolution of each line when compared to the plasticity of the ancestor (see section “Evolution and Plasticity of Gene Expression”). The 10 samples (8 evolved lineages plus two ancestors) were grown in the same incubator for 1 week either at 17 or 23°C prior to collect the cell suspensions. Experiments were performed in duplicate at these two temperature conditions, in order to create two independent biological replicates. The full list of 40 RNA-seq samples is described in Additional File 1: Supplementary Table S1. Prior to RNA extraction of the 40 samples, 50 mL of snap frozen cell suspension were lyophilized for 3 days (LYOVAC^TM GT 2-E freeze dryer, Steris) and ground in liquid nitrogen to a fine powder with a mortar and pestle. Lyophilized tissue was homogenized in 1 mL of Trizol reagent (Invitrogen Inc.), and 200 μL of chloroform, vortexed, incubated at room temperature for 3 min and centrifuged at 12,000 g for 15 min at 4°C. Precipitation of RNA was done by transferring the aqueous phase to a new tube and adding 100% isopropanol. Samples were vigorously shaken and incubated at room temperature for 10 min. RNA pellets were obtained after centrifugation at 12,000 g for 10 min and washed with 1 mL of 70% ethanol. After centrifugation at 7,500 g for 5min, RNA was washed in 70% ethanol, air dried and re-suspended in RNase-free water and stored at −80°C. Sampling quality was checked using agarose gels, Qbit quantification (RNA HS assay kit, Molecular Probes Inc., United Kingdom) and Fragment Analyzer (Agilent, CA, United States). Libraries were created using TruSeq Stranded mRNA Sample Prep kit from Illumina by selecting fragments between 250 and 400 bp captured with poly-dT oligonucleotides (Integragen Inc., Evry, France). Paired-end (2 × 75 bp) sequencing was done with HiSeq4000 Illumina (Integragen Inc., Evry, France).

Transcriptome Mapping and Quantification

We examined raw reads using FastQC (version 0.11.5)¹, removed adaptors and performed quality-based trimming using Trimmomatic (version 0.32; Bolger et al., 2014) with the following options: PE -phred33 SLIDINGWINDOW:4:20 MINLEN:30 TOPHRED33. For each sample, mates of trimmed reads were mapped to the reference genome IPO-323² using HISAT2 (version 2.0.4; Pertea et al., 2016). The default parameters were used except for the intron length, the number of multi mapping sites allowed per read and the stringency of mapping score, as follow: -min-intronlen 20 -max-intronlen 15,000 -k 3 –score-min L, 0, -0.4. We quantified the number of mapped reads using the annotation published by Grandaubert et al. (2015) without annotated sequence from transposable elements, using Stringtie with the following options: -p 16 -c 7 -m 200 -e –B –G (version 2.3; Pertea et al., 2016). Reads that mapped to the predicted 18S and 28S ribosomal units were all removed from the analysis. The final number of genes utilized for all analyses was 10950.

Data Normalization

We performed read count normalization using a size factor calculated for each sample according to the RLE (Relative Log Expression) method implemented in the R package DESeq2 (Love et al., 2014). To identify expressed genes among different RNAseq samples coming from different batches of sequencing (strain IPO-323), we used Fragments Per Kilobase of transcript per Million mapped reads (FPKM) normalization.

Coding Variant Calling From RNA-Sequencing

SNP calling for each RNA-sequencing sample was performed using GTAK GenomeAnalysisTK.jar with the following parameters: –T SplitNCigarReads –rf ReassignOneMappingQuality –RMQT 60 –U ALLOW_N_CIGAR_READS and GenomeAnalysisTK.jar –T HaplotypeCaller –dontUseSoftClippedBases –stand_call_conf 20 (DePristo et al., 2011). Custom R scripts were used to compare coding SNPs among samples.

Natural Variation of Gene Expression Level Among Fungal Isolates

Gene expression levels were compared between our two ancestor clones and the reference strain IPO-323 (Dutch field strain isolated in 1984), using in vitro cultures RNA-seq data from the literature (using Yeast Malt Sucrose medium at 18°C (Kellner et al., 2014); using PDB medium at 18°C (Rudd et al., 2015); (fastq files SRR1167717, SRR1167718³; and fastq files ERR789217, ERR789218 kindly provided by J. Rudd). Raw fastq files from both studies were analyzed the same way as for our samples.

The following model was fit for each gene to compare the genotypes:

Y_ij is the expression level of the ith genotype (i = IPO-323, MGGP01, MGGP44) and the jth replicate (j = 1, 2, 3, 4), μ is the intercept, and ε_ij is the random error. Significant effects were detected by comparing the full model to its null model without the genetic effect, using LRT and contrasts tests. Correction for multiple testing was performed using a FDR at 5% (Benjamini and Hochberg, 2005).

Differential Gene Expression Analysis Among Evolved Lineages

To compare the gene expression levels among all evolved lineages the following DESeq2 model was fit for each gene:

Y_lskj is the observed expression level of the lth genetic background (l = MGGP44, MGGP01), the sth selection regime (s = fluctuating, stable 17, stable 23), the kth assay temperature (k = 17°C, 23°C) and the jth replicate (j = 1,2), μ is the intercept and ε is the random error. The significance of each term was tested using a LRT by comparing the full model with a reduced model without the considered term. Correction for multiple testing was done using FDR at 5%.

To observe the relative contribution of genetic background, temperature and selection regimes effects, we also transformed the read counts with the DESeq2 rlog function to perform a Principal Component Analysis (PCA) with the DESeq2 plotPCA function. We computed the PC scores of each sample along the first four principal components. For each axis PC scores were then analyzed with the following ANOVA model using lm function in R:

PC_lskj is the PC score of the lth genetic background (l = MGGP44, MGGP01), for the sth selection regime (s = fluctuating, stable 17, stable 23), at the kth assay temperature (k = 17°C, 23°C) and the jth replicate (j = 1, 2), μ is the intercept and ε is the random error.

Evolution and Plasticity of Gene Expression

We analyzed the transcriptome evolution and the level of plasticity of gene expression for each ancestral genotype separately. Significant differences in the interaction term between lineages and temperature could not be tested with the DESeq2 Model (2), due to model overfitting. Here, the following model which includes the gene expression data from the ancestors was fit for each gene separately with the R package DESeq2:

Y_ijk is the expression level of the ith lineage (i = ancestor MGGP01, 12F, 13F, 1217, 1323, or ancestor MGGP44, 441F, 443F, 44117, 44323), at the kth temperature (k = 17°C, 23°C), and for the jth replicate (j = 1, 2), μ is the intercept, and ε_ikj is the random error. LRT tests were used to compare the full model to reduced models, and FDR threshold at 5% was used for multiple testing correction.

To identify gene expression evolution due to fluctuation, significant results were classified based on the following rationale. Among genes that showed a significant lineage effect, we considered those with a significant contrast between evolved lineage and their ancestor. Genes exclusively evolving under fluctuation were pulled out when satisfying these two requirements: (1) significant contrasts between fluctuating and stable regimes (2) not significant contrasts between stable regime and ancestor.

To identify plastic genes, we considered genes with significant temperature effect with a fold change of read count greater than 2. To analyze the evolution of plasticity and compare plasticity between the two ancestors, fluctuating and stables lineages, we pulled genes using contrasts.

Gene Ontology Enrichment Tests

To determine whether particular classes of gene function were overrepresented within the different sets of genes showing significant effects, we performed a Gene Ontology enrichment analysis for each set separately, using the R package GoSeq (Young et al., 2010). Note that 5690 genes are associated with unknown function in the current annotation of Z. tritici.

Gene Co-expression Network Analysis

We searched for correlated patterns of expression among genes across the selection regimes. Only genes that displayed at least one significant contrast for the selection regimes from our DESeq2s Model (3) were included in this analysis (3174 and 3440 genes, for the genetic background MGGP01 and MGGP44, respectively). Scale-free co-expression networks were built using a Weighted Gene Correlation Network Analysis with the R package WGCNA (weighted correlation network analysis) (Langfelder and Horvath, 2008). Using the log2 transformed read counts, the genes were clustered into modules according to their expression profile in response to the selection regimes and the assay temperatures. Modules were summarized by the first principal component of the gene expression data of the module (the eigengene), and the hub (most highly connected node of the module). All genes within modules were ranked using their correlation coefficient to the eigengene of the module.

Statistical Analyses

All statistical tests were performed with R (version 3.3.1)⁴. In addition to the statistical models used for differential expression analyses, chi-square tests were used to compare the number of genes showing significant effects among chromosomes. Agreement testing between replicates was done using Cohen’s weighted kappa from “psych” R package. To do so we grouped raw read counts into 11 levels used as ordinal categories as follows: < 10 reads, between 10 and 25 reads, between 25 and 50 reads, between 50 and 100, between 100 and 250, between 250 and 500, between 500 and 1000, between 1000 and 2500, between 2500 and 10,000 and greater than 10,000.

Welch’s two samples T-tests were used to compare the length between up and down-regulated genes which expression level evolved after fluctuation. For the analysis of gene expression plasticity, chi-square tests with Holm’s correction for multiple testing were applied to compare changes of gene number among the different selection regimes. Following the WGCNA permutation tests (10,000 bootstraps) were done to test for the significance of the ranking position of candidate genes with functional annotation related to transcription regulation or signal transduction. Two bootstraps were done, using a sampling of 24 (genes involved in transcriptional regulation) or 39 (genes involved in signal transduction) randomly chosen genes among 4786 (total number of genes in the 8 selected modules revealed by WGCNA, Additional File 11: Supplementary Table S7).

Results

Using two founder clones (hereafter labeled as MGGP01 and MGGP44), we performed a serial transfer experiment for 48 weeks, testing three abiotic conditions: two stable regimes (temperature of 17 or 23°C), and one fluctuating regime by alternating between these two temperatures every 2.5 days (Figure 1). At the end of the experiment, the transcriptome of 8 evolved lineages and ancestral clones was sequenced after 1 week multiplication at both low (17°C) and high (23°C) temperatures (see list of samples in Additional File 1: Supplementary Table S1).

Mapping Quality and Normalization

The mapping rate on the genome from the reference strain IPO-323 (Goodwin et al., 2011), was greater than 83.5% for each sample. Average count of mapped reads was 35 million after relative log expression normalization, and similar distributions of normalized read count was observed (Additional File 2: Supplementary Figures S1, S2). Correlation coefficients between our biological replicates (Kendall τ) were on average 0.89, and agreement between those replicates was good, with Cohen’s weighted κ ranging from 0.86 to 0.97 (McIntyre, 2011), making the repeatability suitable for further analysis (Additional File 2: Supplementary Figure S3). Only 10% of the genes changed in expression by more than 2-fold between biological replicates. We tested our statistical models with and without this set of genes and we found very strong correlations of P-values. We thus included the full set of genes in our analyses presented here. Notably, in comparison to growth at 17°C, growth at 23°C resulted in a greater variation of gene expression across replicates (CV = 0.021 versus CV = 0.049, at 17 and 23°C, respectively).

Genome-Wide Analysis of the Relative Abundance of Transcripts Between Acclimated Laboratory Strains

We checked the quality of our experiment by comparing our RNA-seq data with two previously published studies on the reference strain IPO-323 (Kellner et al., 2014; Rudd et al., 2015). For this purpose we used the transcriptomes from the two ancestor clones MGGP01 and MGGP44. In total, 61% of the genes were transcribed in all three genotypes (with 6663 common gene transcripts on the core chromosomes and 39 common gene transcripts on the accessory chromosomes among 10,950 genes (Table 1). Correlations between genotypes were examined using FPKM data (Additional File 3: Supplementary Figure S4). As expected, the highest correlation was found between the two ancestor clones coming from the same population (Pearson’s coefficient of 0.76). Comparisons of gene expression on CCs show congruent patterns between pairs of genotypes (Additional File 3: Supplementary Figure S5). Nevertheless, 57% (3768 out of 6663) of common transcripts were differentially expressed between genotypes (DESeq2 model 1) (Table 1). Some accessory chromosomes showed no transcriptional activity (ACs 16 and 21 for the genetic background MGGP01 and ACs 18 and 21 for MGGP44). Overall, ACs possessed a significantly lower proportion of expressed genes than CCs (Chi-2 test using a correction for the number of genes per chromosomes, P < 2.2 × 10^–16). Among CCs chromosome 7 was singular. This chromosome contained significantly less gene transcripts than the other CCs, except for the smallest core chromosome 13 (pairwise Chi-2 tests with a correction for the gene number, P < 0.025, except for CC7 vs. CC13, P = 0.057). A large fragment of 800 kb on the distal region of the left arm of chromosome 7 was without any transcriptional activity in all three genotypes (Additional File 3: Supplementary Figure S6). In conclusion, in vitro transcriptome between strains is similar for the identity of gene transcripts, but highly variable for the level of gene expression.

Contrasted Transcriptome Evolution Between Genetic Backgrounds

To evaluate whether there was any differences of gene expression among the three selection regimes, we first compared the transcriptomes without the ancestor data. The genetic background and the assay temperature were the two main factors contributing to the total variance among evolved lineages (Table 2) (Principal Component Analysis, with 51% (98% of PC1) and 11.3% (66% of PC2 and 20% of PC3), respectively). The contribution of the selection regime (stable 17°C, stable 23°C, or fluctuating) was smaller with 6% of total variance (30% of PC3 and 66% of PC4) but significant. We then used our DESeq2 model (2) to identify the set of genes differentially expressed between both genotypes, between selection regimes, and between temperature treatments. We found 7022 genes that showed significant differences of expression between the two genetic backgrounds, 4146 genes that differed across selection regimes and 5637 genes that differed due to the assay temperature (Figure 2). In addition, we detected 1728 differently expressed genes for the interaction term between genetic backgrounds and selection regimes. More specifically, many of those genes showed opposite direction of differential expression between the two genetic backgrounds (Additional File 4: Supplementary Figure S7). In sum, the expression level of a large number of genes evolved differently in both genetic backgrounds. As a consequence, to gain full insight in gene expression changes in response to the selection regime, we considered the two ancestral genetic backgrounds separately in further analyses.

FIGURE 2

Evolution of Gene Expression During Adaptation to Fluctuation

We investigated the evolution of gene expression by including the data from the ancestor (DESeq2 model 3). Evolution of gene expression in response to selection was inferred from significant differences between ancestral and evolved transcriptomes.

Thirty two percent of the genes (3420 and 3659 genes for MGGP01 and MGGP44, respectively) were differentially expressed across selection regimes (Table 3). The assay temperature also affected the level of gene expression, for up to 39% of the transcriptome. The distribution of the fold-change of gene expression for all genes showing a significant effect is given in the Additional File 5: Supplementary Figure S8. The number of genes affected by fluctuation was approximately the same in both genetic backgrounds. Nonetheless, the number of up- and down-regulated genes was asymmetric between the two genetic backgrounds, with a majority of evolution toward down-regulation for MGGP01, and an opposite pattern for MGGP44 (Table 3, Additional File 5: Supplementary Figures S8, S9). As the transcriptome was sequenced in two independent fluctuating lines for each background, we were able to measure the amount of parallel evolution in our experimental setting. Significant change in expression compared to the ancestor was highly correlated in parallel lineages for each genetic background with a large proportion (nearly 60%) of common genes between the two independently evolved lines (568 out of 907 Differentially Expressed Genes (DEGs) for MGGP01, and 366 out of 638 DEGs for MGGP44 when considering genes with a fold change greater than two) (Figure 3). We focused our attention on this set of genes, assuming that parallel evolutionary change could be the consequence of adaptive evolution in fluctuating environments. When comparing the two genetic backgrounds we detected an overlap of 70 DEGs due to fluctuation, with 45 genes evolving in the same direction. Notably, the gene length between up- and down-regulated genes in fluctuating conditions was biased toward a smaller size for down-regulated genes (Welch t-tests, P = 3.8 × 10^–9 and P = 0.003, for MGGP01 and MGGP44, respectively).

FIGURE 3

Among genes significant for the fluctuating regime, a large number were also significant for the stable regimes, to a greater extent for the background MGGP44 (with 90 and 61% of genes in common between fluctuating and stable regimes at 17°C and 23°C, respectively, versus 21 and 41%, for the background MGGP01) (Figure 4). While 20 transcripts were specific to the fluctuation regime for the genetic background MGGP44, 291 genes specifically evolved under fluctuations for the background MGGP01 (gene lists are given in Additional File 6: Supplementary Table S2). This strong contrast between the two genetic backgrounds illustrates the genetic diversity of this fungal species, contributing to the complex genetic basis of adaptation. In this study there is evidence for high level of genetic variation between two ancestors (Figure 5). First, we noticed a presence-absence polymorphism of the mini chromosomes (Figure 5A). While AC 21 seemed to be missing for both strains, AC 16 had no transcripts for MGGP01, and AC 18 had no transcripts for MGGP44. In addition, we detected many SNPs in common transcripts between the two ancestors, with 62112 SNPs between MGGP01 and MGGP44, and 22384 SNPs between IPO-323 and the two ancestors (Figure 5B).

FIGURE 4

FIGURE 5

Gene Ontology Tests

Gene Ontology (GO) term enrichment tests identified different functional terms between the two genetic backgrounds. First, when considering the 568 genes having evolved in the fluctuation regime for the wild-derived background MGGP01, we did not find any evidence for GO term enrichment. However, when considering only the 284 genes that evolved specifically under thermal fluctuations, we detected functional categories related to Cell cycle control, cell division and chromosome partitioning (P = 0.0019) and Cytoskeleton (P = 0.027). These functional categories include proteins known to be essential components of the cytoskeletal mitotic apparatus and essential regulators of cell cycle progression, such as the gene ID 59920, ortholog of the conserved regulator of mitosis Wee1, which is known to control the timing of entry into mitosis in yeast (Kellogg, 2003; Additional File 6: Supplementary Table S2). For the lineages from the background MGGP44, when considering all significant genes for the fluctuating regime (which also appeared significant for the stable regimes), the GO term defense mechanisms was found (9 genes, P = 0.013, with 6 genes up-regulated and 3 genes down-regulated, see Additional File 6: Supplementary Table S2).

Genome Wide Distribution of DEGs

We further examined the genomic distribution of differentially expressed genes due to the selection regimes. First, we contrasted the number of differentially expressed genes between CCs and ACs. For the lineages derived from the background MGGP01 only, the proportion of DEGs normalized by the total number of genes per chromosome was higher on ACs than on CCs (Additional File 7: Supplementary Table S3). Next, we surveyed the density of DEGs along each chromosome using a non-overlapping window approach (Figure 6). The distribution of the genes with significant change of expression during adaptation to stable and fluctuating regimes was not uniform, as we found several regions where the gene density was high. Unexpectedly, these hotspots with a higher density of DEGs –normalized by the number of known annotated genes– associated with thermal adaptation were sometimes located near the tips of chromosomes: 2, 6, 7, 8, and 13 for the CCs or across the ACs with mainly chromosomes 15, 17, 18, and 19. In contrast, plastic genes that showed a significant temperature effect using the same model (DESeq model 3) were heterogeneously distributed along the genome with no clear evidence of hotspots (Figure 6). The enrichment of accessory chromosomes in genes which expression changed during experimental evolution suggests that these chromosomes could be involved in adaptation to abiotic environment. Overall, about 40 genes had their expression turned off after experimental evolution (only 2 genes on ACs). None of these genes were located in hotspot regions, thus excluding the hypothesis that the hotspot regions could correspond to deletion events. As a matter of fact, 80% of the genes located in hotspots were up-regulated compared to the ancestral gene expression level. These hotspots often collocate with regions known to be enriched in transposable elements and are often near subtelomeric regions of core chromosomes (Schotanus et al., 2015). This result suggests that various evolutionary scenarios, including a rapid turn-over of gene regulation in these regions where transposition during mitosis, could impact gene expression. Such hypotheses remain difficult to test without further experimental evidence.

FIGURE 6

Gene Network Rewiring

We performed a weighted correlation network analysis to search for clusters of co-regulated genes (details in “Materials and Methods” section). A total of 16 and 15 co-expression modules were detected, respectively, for MGGP01 and MGGP44 lineages (Additional File 8: Supplementary Tables S4, S5). We retained 8 of the largest co-expression modules that include most of the genes which were significant for the fluctuating regime from our model (3) (63 and 92%, for MGGP01 and MGGP44, respectively) (Additional File 9: Supplementary Figure S10). The goal of this approach was to identify hub genes within co-expression modules, defined as the closest gene from the module average (Additional File 10: Supplementary Table S6). Several molecular functions were involved, including histone modification and phospholipase activity. Interestingly, among the 8 modules, we found 23 genes involved in transcriptional regulation and chromatin remodeling (mainly transcription factors) and 39 genes involved in signal transduction, including 11 genes with a kinase activity (Additional File 11: Supplementary Table S7). However, the ranking position of these genes within modules was not significantly biased toward the most connected genes (Wilcoxon tests with empirical permutation testing, P = 0.98 and P = 0.99, for transcription factors and signal transduction genes, respectively). When considering the relationship between the physical (genomic) distance between genes within modules with coexpression, no large cis- or trans-regulated clusters of genes were detected. Co-expressed transcripts within modules were rather scattered across the whole genome (Additional File 9: Supplementary Figure S11).

Evolution of Gene Expression Plasticity

The effect of the assay temperature in our laboratory growth conditions – the gene expression difference between 17°C and 23°C – was large and pervasive across the genome (DEseq Model 3, see Table 3). In order to make sure that differences in gene expression between temperatures could be attributed to plasticity, we verified that the genetic composition of the populations was identical at both temperatures for the genes influenced by temperature. Indeed, in a polymorphic population, where several clones could be competing at the end of the experimental evolution, we can assume that there is a confounded effect between clone competition at the two temperatures and actual gene expression plasticity. In other words, if the clone fitness varies between the two temperatures, different gene expression level could be due to changes in genotype frequency rather than single genotype differential expression. We observed new mutations in transcribed regions from the RNA sequences and measured mutation frequencies at both temperatures. We found no evidence for a change in haplotype frequency between samples (Additional File 12: Supplementary Table S8). Based on these results, we made the assumption that changes of gene expression were associated with phenotypic plasticity.

The proportion of plastic genes between the two ancestors was significantly different. For a total of 7804 plastic genes, we report 653 versus 1348 plastic genes for MGGP01 and MGGP44, respectively (P = 2.2 × 10^–16) (Figure 7). When considering genes with a raw fold change greater than two, we detected 501 versus 699 plastic genes for MGGP01 and MGGP44, respectively (P = 4.4 × 10^–8). Among those plastic genes, 92 were common between the two ancestors. Interestingly, the distribution of the expression fold-change between the two temperatures was biased toward up-regulated genes at 17°C for MGGP01. Opposite result was found for MGGP44, with more up-regulated genes at 23°C (Additional File 13: Supplementary Table S9).

FIGURE 7

We identified 3.8 and 7.8% of the transcriptome displaying a significant Temperature-by-Regime interaction in at least one of the three selection regime, for wild-derived background MGGP01 and MGGP44, respectively. These genes were thus associated with an evolutionary change in phenotypic plasticity, which was further classified as a gain or loss of plasticity in comparison with the ancestral plasticity. In both genetic backgrounds, the fluctuating lines acquired less plastic genes than stable lines (Chi-square tests with Holm’s correction for multiple comparisons, all P < 2.10^–7), with a significant loss of plastic genes when compare to stable lines in the background MGGP01 (Chi-square tests with Holm’s correction for multiple comparisons, all P < 2.10^–3) (Figure 7). These findings suggest that in our laboratory conditions, adaptation to thermal fluctuations is not associated with the evolution of enhanced plasticity of gene expression. More surprisingly perhaps, adaptation to fluctuations tends to be associated with the loss of plasticity for a substantial number of genes that were plastic in the ancestor, suggesting that plasticity could even be maladaptive in a fluctuating environment.

Discussion

Intraspecific Variation of Gene Expression

By comparing three genotypes of Z. tritici we found pervasive gene expression variation most likely due to the contribution of numerous variants (small or large scale mutations) including variation of histone modification. Substantial amount of standing genetic variation is present in populations of this species (nucleotide diversity pi = 0.01) (Hartmann et al., 2018). Although transcriptomic comparisons are often made using a limited number of strains, there is accumulating evidence for a considerable variation of gene expression within and between populations in Z. tritici (Palma-Guerrero et al., 2017; Haueisen et al., 2019). In this study we evidenced that 57% of the transcriptome varies among genotypes. This level of variation leaves room for adaptation and is comparable to the amount of natural variation already reported within eukaryotic species, mostly using eQTL mapping (e.g., Schadt et al., 2003; Huang et al., 2015; Leder et al., 2015; Zan et al., 2016; Osada et al., 2017). Interestingly, we also found a large silenced region on chromosome 7 in all three wild genotypes. This region was previously identified for the strain IPO-323 by Kellner et al. (2014); Rudd et al. (2015), and Schotanus et al. (2015). This region is known to carry histone H3K27me3 modifications that mediate transcriptional silencing, also present on the ACs. This chromosomal fragment most likely originated from a fusion of a whole accessory chromosome onto the original chromosome 7 (Schotanus et al., 2015), and is in fact conserved among the three wild genotypes compared in this study. Likewise, Haueisen et al., 2019 recently confirmed the conservation of this silenced region among several other wild isolates.

The Effect of Fluctuating Selection on the Transcriptome

Evolutionary change in gene expression was pervasive, with as much as 38% of Z. tritici transcripts displaying different expression levels among evolved lineages. In our experimental conditions evolved lineages had on average several hundreds of genes which expression level changed by at least a factor 2. Transcriptome evolution in fluctuating conditions was shown to be repeatable between parallel lineages initiated from the same clone. However, our results clearly show that evolution was not parallel when using different genetic backgrounds, with little overlap between differentially expressed genes. Interestingly, different functional categories were over-represented between the two genetic backgrounds, suggesting different evolutionary routes can be explored to respond abiotic changes. In contrast to MGGP01, the background MGGP44 showed a very few genes responding solely to the fluctuating regime. Beside the fact our strains were acclimated for several weeks prior to the experiment, we cannot exclude that many of those genes indeed evolved in response to the laboratory selection. Whether the ancestor MGGP44 was already adapted to fluctuating environment (e.g., through its initial large number of plastic genes), or that adaptation to fluctuations may not be based on a different set of genes than adaptation to stable conditions, remains to be determined. To sum up these results illustrate that this fungal species can evolve fast by means of gene expression variation operating on many functional categories, but the evolutionary outcome is different between genetic backgrounds.

Genome-Wide Distribution of Regulatory Changes

The genome-wide distribution of differentially expressed genes due to selection suggests that a few genomic regions are more prone to rapid evolution. Compartmentalization of filamentous fungal genomes is a well-known feature, characterized by the presence of regions of distinct evolutionary rates (Dong et al., 2015; Frantzeskakis et al., 2019). In our experiment, regions near subtelomeres on a few core chromosomes and some accessory chromosomes were more prone to rapid evolution of gene expression. This result has important implications for understanding the evolutionary potential of filamentous fungi. Accessory chromosomes are often qualified as a cradle for adaptive evolution, due to their high mutation rates and structural rearrangements (Croll and McDonald, 2012). We provide here the first experimental evidence of rapid evolution of gene expression in this genome “accessory” compartment. However the molecular causes of gene expression changes remained to be understood. The regions enriched in genes showing significant effect of the selection regimes (tips of CCs and large portion of ACs), tend to contain more transposable elements than the genome average (Schotanus et al., 2015). The role of transposable elements in regulating gene expression evolution in Eukaryotes has been well documented (see Schrader and Schmitz, 2019 for recent review). These results suggest that further investigations are needed to identify the causal factors –TE themselves, epigenetic regulation, or SNPs/indels– that are associated with the evolution of gene expression level.

Genome Stability During Asexual Multiplication

Loss of accessory chromosome has been recently reported for 4-week in vitro grown cells in yeast-malt-sucrose medium in Z. tritici (Möller et al., 2018). In our study we found no evidence for a loss of accessory chromosomes (among genes significant for selection the regimes, two were turned off on accessory chromosomes). However, we cannot exclude the hypothesis that some of the cells in evolved lineages lost accessory chromosomes, thus falsely leading to a down-regulated effect on gene expression. Likewise, up-regulation of gene expression could obscure a loss of ACs, and TE mobility could influence the genome structure and falsely lead to an increase of gene expression that is in fact associated with copy number variation.

Fighting Changes by Staying Stable

We observed significant differences of plasticity between both ancestral genotypes, highlighting genetic variation for transcription plasticity. Genotype-by-environment interaction for gene expression is common in many organisms (e.g., yeast (Landry et al., 2006; Smith and Kruglyak, 2008; Li et al., 2010), Arabidopsis genus (He et al., 2016). Plasticity is often thought to be an adaptive response to environmental heterogeneity, that is likely to be lost in stable environments where it is likely to be unnecessary (Lande, 2009). In this context, our observation that gene expression plasticity tends to be lost in all lineages exposed to fluctuating selection regimes appears at odds with theory. This result indicates that in our experimental conditions robustness of gene expression is favored under environmental fluctuations. The evolutionary consequences of robustness are important. Robustness can be adaptive and lead to the accumulation of cryptic genetic mutations that can increase the evolvability of the species facing new environments (Le Rouzic and Carlborg, 2008; Cuypers et al., 2017; Payne and Wagner, 2019). Several theoretical approaches have explained a gain of robustness, which is dependent on the rate of fluctuation and the strength of selection (Siegal and Bergman, 2002; Le Rouzic et al., 2013). Full understanding of the molecular basis of robustness remains largely unexplored, with a clear lack of empirical support. Robustness could rely on functional redundancy arising from gene duplications, or alternatively on specific network topological features. Unfortunately, functional annotation of the genome of Z. tritici is not advanced enough yet to gain full insights into these questions. A recent empirical study in Arabidopsis thaliana did not support the above hypothesis and demonstrated that a single gene, rather than gene network connectivity or functional redundancy, was associated to trait robustness (Lachowiec et al., 2018). In contrast, a former study evidenced robustness by means of functional redundancy in yeast (Li et al., 2010).

Conclusion

Our findings demonstrate substantial evolution of gene expression under fluctuating temperature with striking differences between genetic backgrounds. Notably, we identified a set of genes involved in gene expression, chromatin regulation and signal transduction genes. The ecological relevance of those genes in the adaptation to temperature variation remains to be elucidated. In addition our results suggest an important role of gene expression evolution on the accessory chromosomes, which has not been reported before for this species. Last, we found that under fluctuating temperature the transcriptome evolved toward a loss of plasticity. These findings suggest more work is needed to further characterize the role of gene expression plasticity in adaptation.

References

• 1

  BenjaminiY.HochbergY. (2005). Controlling the false discovery rate: a practical and powerful approach to multiple testing.J. R. Statist. Soc. B.57289–300. 10.1111/j.2517-6161.1995.tb02031.x

  • CrossRef
  • Google Scholar

• 2

  BennettA. F.LenskiR. E.MittlerJ. E. (1992). Evolutionary adaptation to temperature. I. Fitness responses of Escherichia coli to changes in its thermal environment.Evolution4616–30. 10.1111/j.1558-5646

  • CrossRef
  • Google Scholar

• 3

  BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: a flexible trimmer for illumina sequence data.Bioinformatics302114–2120. 10.1093/bioinformatics/btu170

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  BremR. B.KruglyakL. (2002). Genetic dissection of transcriptional regulation in budding yeast.Science296752–755. 10.1126/science.1069516

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  BrittenR. J.DavidsonE. H. (1971). Repetitive and non-repetitive DNA sequences and a speculation on the origins of evolutionary novelty.Q. Rev. Biol.46111–138. 10.1086/406830

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  ChenJ.NolteV.SchlöttererC. (2015). Temperature-related reaction norms of gene expression: regulatory architecture and functional implications.Mol. Biol. Evol.322393–2402. 10.1093/molbev/msv120

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  CorlA.BiK.LukeC.ChallaA. S.SternA. J.SinervoB.et al (2018). The genetic basis of adaptation following plastic changes in coloration in a novel environment.Curr. Biol.282970–2977. 10.1016/j.cub.2018.06.075

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  CrollD.McDonaldB. A. (2012). The accessory genome as a cradle for adaptive evolution in pathogens.PLoS Patho.8:e1002608. 10.1371/journal.ppat.1002608

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  CuypersT. D.RuttenJ. P.HogewegP. (2017). Evolution of evolvability and phenotypic plasticity in virtual cells.BMC Evol. Biol.17:60. 10.1186/s12862-017-0918-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  DePristoM.BanksE.PoplinR.GarimellaK.MaguireJ.HartlC.et al (2011). A framework for variation discovery and genotyping using next-generation DNA sequencing data.Nat. Genet.43491–498. 10.1038/ng.806

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  DharR.SägesserR.WeikertC.WagnerA. (2013). Yeast adapts to a changing stressful environment by evolving cross-protection and anticipatory gene regulation.Mol. Biol. Evol.30573–588. 10.1093/molbev/mss253

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  DongS. M.RaffaeleS.KamounS. (2015). The two-speed genomes of filamentous pathogens: waltz with plants.Curr. Opin. Genet. Dev.3557–65. 10.1016/j.gde.2015.09.001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  EnardW.KhaitovichP.KloseJ.ZollnerS.HeissigF. (2002). Intra- and interspecific variation in primate gene expression patterns.Science296340–343. 10.1126/science.1068996

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  FayJ. C.McCulloughH. L.SniegowskiP. D.EisenM. B. (2004). Population genetic variation in gene expression is associated with phenotypic variation in Saccharomyces cerevisiae.Genome Biol.5:R26. 10.1186/gb-2004-5-4-r26

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  FeugeasJ. P.TourretJ.LaunayA.BouvetO.HoedeC.DenamurE.et al (2016). Links between transcription, environmental adaptation and gene variability in Escherichia coli: correlations between gene expression and gene variability reflect growth efficiencies.Mol. Biol. Evol.332515–2529. 10.1093/molbev/msw105

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  FischerE. M.KnuttiR. (2015). Anthropogenic contribution to global occurrence of heavy-precipitation and high-temperature extremes.Nat. Clim. Chang.46560–564. 10.1038/NCLIMATE2617

  • CrossRef
  • Google Scholar

• 17

  FisherK. J.LangG. I. (2016). Experimental evolution in fungi: an untapped resource.Fungal Genet. Biol.9488–94. 10.1016/j.fgb.2016.06.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  FonesH.GurrS. (2015). The impact of Septoria tritici Blotch disease on wheat: an EU perspective.Fungal Genet Biol.793–7. 10.1016/j.fgb.2015.04.004

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  FoxR. J.DonelsonJ. M.SchunterC.RavasiT.Gaita ìn-EspitiaJ. D. (2019). Beyond buying time: the role of plasticity in phenotypic adaptation to rapid environmental change.Phil. Trans. R. Soc. B37420180174. 10.1098/rstb.2018.0174

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  FranciscoC. S.MaX.ZwyssigM. M.McDonaldB. A.Palma-GuerreroJ. (2019). Morphological changes in response to environmental stresses in the fungal plant pathogen Zymoseptoria tritici.Sci. Rep.9:9642.

  • Google Scholar

• 21

  FrantzeskakisL.KuschS.PanstrugaR. (2019). The need for speed: compartmentalized genome evolution in filamentous phytopathogens.Mol. Plant Pathol.203–7. 10.1111/mpp.12738

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  GarlandT.RoseM. R. (2009). Experimental Evolution: Concepts, Methods, and Applications of Selection Experiments. Berkeley, CL: University of California Press.

  • Google Scholar

• 23

  GhalamborC. K.HokeK. L.RuellE. W.FischerE. K.ReznickD. N.HughesK. A. (2015). Non-adaptive plasticity potentiates rapid adaptive evolution of gene expression in nature.Nature525372–375. 10.1038/nature15256

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  GhalamborC. K.McKayJ. K.CarrollS. P.ReznickD. N. (2007). Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments.Funct. Ecol.21394–407. 10.1111/j.1365-2435.2007.01283.x

  • CrossRef
  • Google Scholar

• 25

  Glaser-SchmittA.ParschJ. (2018). Functional characterization of adaptive variation within a cis-regulatory element influencing Drosophila melanogaster growth.PLoS Biol.16:1–28. 10.1371/journal.pbio.2004538

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  GoodwinS. B.M’BarekS.DhillonB.WittenbergA. H. J.CraneC. F.HaneJ. K.et al (2011). Finished genome of the fungal wheat pathogen Mycosphaerella graminicola reveals dispensome structure, chromosome plasticity, and stealth pathogenesis.PLoS Genet.7:e1002070. 10.1371/journal.pgen.1002070

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  GrandaubertJ.BhattacharyyaA.StukenbrockE. H. (2015). RNA-seq-based gene annotation and comparative genomics of four fungal grass pathogens in the genus Zymoseptoria identify novel orphan genes and species-specific invasions of transposable elements.G351323–1333. 10.1534/g3.115.017731

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  HartmannF. E.McDonaldB. A.CrollD. (2018). Genome-wide evidence for divergent selection between populations of a major agricultural pathogen.Mol. Ecol.272725–2741. 10.1111/mec.14711

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  HartmannF. E.Sánchez-ValletA.McDonaldB. A.CrollD. (2017). A fungal wheat pathogen evolved host specialization by extensive chromosomal rearrangements.ISME J.111189–1204. 10.1038/ismej.2016.196

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  HaueisenJ.MöllerM.EschenbrennerC. J.GrandaubertJ.SeyboldH.AdamiakH.et al (2019). Highly flexible infection programs in a specialized wheat pathogen.Ecol. Evol.9275–294. 10.1002/ece3.4724

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  HeF.ArceA. L.SchmitzG.KoornneefM.NovikovaP.BeyerA.et al (2016). The footprint of polygenic adaptation on stress-responsive cis-regulatory divergence in the Arabidopsis Genus.Mol. Biol. Evol.332088–2101. 10.1093/molbev/msw096

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 32

  HoW.-C.ZhangJ. (2018). Evolutionary adaptations to new environments generally reverse plastic phenotypic changes.Nat. Comm.9:350. 10.1038/s41467-017-02724-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  HuangW.CarboneM. A.MagwireM. M.PeifferJ. A.LymanR. F.StoneE. A.et al (2015). Genetic basis of transcriptome diversity in Drosophila melanogaster.PNAS112E6010–E6019. 10.1073/pnas.1519159112

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  HuangY.AgrawalA. F. (2016). Experimental evolution of gene expression and plasticity in alternative selective regimes.PLoS Genet.12:e1006336. 10.1371/journal.pgen.1006336

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  HughesB. S.CullumA. J.BennettA. F. (2007). An experimental evolutionary study on adaptation to temporally fluctuating pH in Escherichia coli.Physiol. Biochem. Zool.80406–421. 10.1086/518353

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  JalletA. J.Le RouzicA.And GenisselA. (2019). Evolution and plasticity of the transcriptome under temperature fluctuations in the fungal plant pathogen Zymoseptoria tritici.bioRxiv725010.

  • Google Scholar

• 37

  JovicK.SterkenM. G.GrilliJ.BeversR. P. J.RodriguezM.RiksenJ. A. G.et al (2017). Temporal dynamics of gene expression in heat-stressed Caenorhabditis elegans.PLoS ONE12:e0189445. 10.1371/journal.pone.0189445

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  KellnerR.BhattacharyyaA.PoppeS.HsuT. Y.BremR. B.StukenbrockE. H. (2014). Expression profiling of the wheat pathogen Zymoseptoria tritici reveals genomic patterns of transcription and host-specific regulatory programs.Genome Biol. Evol.61353–1365. 10.1093/gbe/evu101

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  KelloggD. R. (2003). Wee1-dependent mechanisms required for coordination of cell growth and cell division.J. Cell Sci.1164883–4890. 10.1242/jcs.00908

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  KenkelC. D.MatzM. V. (2016). Gene expression plasticity as a mechanism of coral adaptation to a variable environment.Nat. Ecol. Evol.7:14. 10.1038/s41559-016-0014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  KingM.-C.WilsonA. C. (1975). Evolution at two levels in humans and chimpanzees.Science188107–116. 10.1126/science.1090005

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  LachowiecJ.MasonG. A.SchultzK.QueitschC. (2018). Redundancy, feedback, and robustness in the Arabidopsis thaliana BZR/BEH gene family.Front. Genet.9:3389. 10.3389/fgene.2018.00523

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  LandeR. (2009). Adaptation to an extraordinary environment by evolution of phenotypic plasticity and assimilation.J. Evol. Biol.221435–1446. 10.1111/j.1420-9101.2009.01754.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  LandryC. R.OhJ.HartlD. L.CavalieriD. (2006). Genome-wide scan reveals that genetic variation for transcriptional plasticity in yeast is biased towards multi-copy and dispensable genes.Gene366343–351. 10.1016/j.gene.2005.10.042

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  LangfelderP.HorvathS. (2008). WGCNA: An R package for weighted correlation network analysis.BMC Bioinform.9:559. 10.1186/1471-2105-9-559

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  Le RouzicA.Álvarez-CastroJ. M.HansenT. F. (2013). The evolution of canalization and evolvability in stable and fluctuating environments.Evol. Biol.40317–340. 10.1007/s11692-012-9218-z

  • CrossRef
  • Google Scholar

• 47

  Le RouzicA.CarlborgÖ (2008). Evolutionary potential of hidden genetic variation.Trends Ecol. Evol.2333-37. 10.1016/j.tree.2007.09.014

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  LederE. H.McCairnsR. J. S.LeinonenT.CanoJ. M.ViitaniemiH. M. (2015). The evolution and adaptive potential of transcriptional variation in sticklebacks–signatures of selection and widespread heritability.Mol. Biol. Evol.32674–689. 10.1093/molbev/msu328

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  LendenmannM. H.CrollD.Palma-GuerreroJ.StewartE. L.McdonaldB. A. (2016). QTL mapping of temperature sensitivity reveals candidate genes for thermal adaptation and growth morphology in the plant pathogenic fungus Zymoseptoria tritici.Heredity116384–394. 10.1038/hdy.2015.111

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  LeroiA. M.LenskyR. E.BennettA. (1994). Evolutionary adaptation to temperature. III. Adaptation of Escherichia coli to a temporally varying environment.Evolution.481222–1229.

  • Google Scholar

• 51

  LevisN. A.PfennigD. W. (2016). Evaluating “plasticity-first” evolution in nature:key criteria and empirical approaches.Trends Ecol. Evol.31563–574. 10.1016/j.tree.2016.03.012

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  LiJ.YuanZ.ZhangZ. (2010). The cellular robustness by genetic redundancy in budding yeast.PLoS Genet.6:e1001187. 10.1371/journal.pgen.1001187

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  LiY.ÁlvarezO. A.GuttelingE. W.TijstermanM.FuJ.RiksenJ. A. G.et al (2006). Mapping determinants of gene expression plasticity by genetical genomics in C. elegans.PLoS Genet.2:2155–2161. 10.1371/journal.pgen.0020222

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.Genome Biol.15:550. 10.1186/s13059-014-0550-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  McIntyreL. M. (2011). RNA-seq:technical variability and sampling.BMC Genomics12:293. 10.1186/1471-2164-12-293

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  MeileL.CrollD.BrunnerP. C.PlissoneauC.HartmannF. E.McDonaldB. A.et al (2018). A fungal avirulence factor encoded in a highly plastic genomic region triggers partial resistance to septoria tritici blotch.New Phytol.2191048–1061. 10.1111/nph.15180

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  MöllerM.HabigM.FreitagM.StukenbrockE. H. (2018). Extraordinary genome instability and widespread chromosome rearrangements during vegetative growth.Genetics210517–529. 10.1534/genetics.118.301050

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  NoormohammadA.RambeaumJ.HeldT.KovakovaV.BergJ.LassigM. (2017). Adaptive evolution of gene expression in Drosophila.Cell Rep.201385–1395. 10.1016/j.celrep.2017.07.033

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  O’DriscollA. O.KildeaS.DoohanF.SpinkJ.MullinsE. (2014). The wheat – Septoria conflict: a new front opening up?Trends Plant Sci.19602–610. 10.1016/j.tplants.2014.04.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 60

  OleksiakM. F.ChurchillG. A.CrawfordD. L. (2002). Variation in gene expression within and among natural populations.Nat. Genet.32261–266. 10.1038/ng983

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  OsadaN.MiyagiR.TakahashiA. (2017). Cis- and trans-regulatory effects on gene expression.Genetics2062139–2148. 10.1534/genetics.117.201459

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 62

  Palma-GuerreroJ.MaX.TorrianiS. F. F.ZalaM.FranciscoC. S.HartmannF. E.et al (2017). Comparative transcriptome analyses in Zymoseptoria tritici reveal significant differences in gene expression among strains during plant infection.MPMI.30231–244. 10.1094/MPMI-07-16-0146-R

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  PayneJ. L.WagnerA. (2019). The causes of evolvability and their evolution.Nat. Rev. Genet.2024–38. 10.1038/s41576-018-0069-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 64

  PerteaM.KimD.PerteaG. M.LeekJ. T.SalzbergS. L. (2016). Transcript-level expression analysis of RNA-seq experiments with HISAT. StringTie and Ballgown.Nat Protoc.111650–1667. 10.1038/nprot.2016.095

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 65

  PigliucciM.MurreC. J.SchlichtingC. D. (2006). Phenotypic plasticity and evolution by genetic assimilation.J. Exp. Biol.2092362–2367. 10.1242/jeb.02070

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 66

  PodrabskyJ. E. (2004). Changes in gene expression associated with acclimation to constant temperatures and fluctuating daily temperatures in an annual killifish Austrofundulus limnaeus.J. Exp. Biol.2072237–2254. 10.1242/jeb.01016

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  RockmanM. V.WrayG. A. (2002). Abundant raw material for cis-regulatory evolution in humans.Mol. Biol. Evol.191991–1994. 10.1093/oxfordjournals.molbev.a004023

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 68

  RuddJ. J.KanyukaK.Hassani-PakK.DerbyshireM.AndongaboA.DevonshireJ.et al (2015). Transcriptome and metabolite profiling of the infection cycle of Zymoseptoria tritici on wheat reveals a biphasic interaction with plant immunity involving differential pathogen chromosomal contributions and a variation on the hemibiotrophic lifestyle definition.Plant Physiol.1671158–1185. 10.1104/pp.114.255927

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 69

  SchadtE. E.MonksS. A.DrakeT. A.LusisA. J.CheN.ColinayoV.et al (2003). Genetics of gene expression surveyed in maize, mouse and man.Nature422292–297. 10.1038/nature01434

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 70

  SchotanusK.SoyerJ. L.ConnollyL. R.GrandaubertJ.HappelP.SmithK. M.et al (2015). Histone modifications rather than the novel regional centromeres of Zymoseptoria tritici distinguish core and accessory chromosomes.Epigenetics Chromatin.8:41. 10.1186/s13072-015-0033-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  SchraderL.SchmitzJ. (2019). The impact of transposable elements in adaptive evolution.Mol. Ecol.281537–1549. 10.1111/mec.14794

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 72

  SiegalM. L.BergmanA. (2002). Waddington’s canalization revisited: Developmental stability and evolution.PNAS9910528–10532. 10.1073/pnas.102303999

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 73

  SikkinkK. L.ReynoldsR. M.ItuarteC. M.CreskoW. A.PhillipsP. C. (2019). Environmental and evolutionary drivers of the modular gene regulatory network underlying ophenotypic plasticity for stress resistance in the nematode Caenorhabditis remanei.G3.9969–982. 10.1534/g3.118.200017

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 74

  SmithE. N.KruglyakL. (2008). Gene-environment interaction in yeast gene expression.PLoS Biol.6:e83. 10.1371/journal.pbio.0060083

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 75

  SørensenJ. G.SchouM. F.KristensenT. N.LoeschckeV. (2016). Thermal fluctuations affect the transcriptome through mechanisms independent of average temperature.Sci Rep.630975. 10.1038/srep30975

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 76

  SteinbergG. (2015). Cell biology of Zymoseptoria tritici: Pathogen cell organization and wheat infection.Fungal Genet. Biol.7917–23. 10.1016/j.fgb.2015.04.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 77

  StewartE. L.CrollD.LendenmannM. H.Sanchez-ValletA.HartmannF. E.Palma-GuerreroJ.et al (2018). Quantitative trait locus mapping reveals complex genetic architecture of quantitative virulence in the wheat pathogen Zymoseptoria tritici.Mol Plant Pathol.9201–216. 10.1111/mpp.12515

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 78

  WhiteheadA.CrawfordD. L. (2006). Neutral and adaptive variation in gene expression.PNAS1035425–5430. 10.1073/pnas.0507648103

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 79

  WittkoppP. J.HaerumB. K.ClarkA. G. (2004). Evolutionary changes in cis and trans gene regulation.Nature43085–88. 10.1038/nature02698

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 80

  YoungM. D.WakefieldM. J.SmythG. K.OshlackA. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias.Genome Biol.11R14. 10.1186/gb-2010-11-2-r14

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 81

  ZanY.ShenX.ForsbergS. K. G.CarlborgÖ (2016). Genetic regulation of transcriptional variation in natural Arabidopsis thaliana accessions.G362319–2328. 10.1534/g3.116.030874

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 82

  ZandveldJ.van den HeuvelJ.MulderM.BrakefieldP. M.KirkwoodT. B. L.ShanleyD. P.et al (2017). Pervasive gene expression responses to a fluctuating diet in Drosophila melanogaster: the importance of measuring multiple traits to decouple potential mediators of life span and reproduction.Evolution712572–2583. 10.1111/evo.13327

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 83

  ZhanJ.McDonaldB. A. (2011). Thermal adaptation in the fungal pathogen Mycosphaerella graminicola.Mol. Ecol.201689–1701. 10.1111/j.1365-294X.2011.05023.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 84

  ZhongZ.MarcelT. C.HartmannF. E.MaX.PlissonneauC.ZalaM.et al (2017). A small secreted protein in Zymoseptoria tritici is responsible for avirulence on wheat cultivars carrying the Stb6 resistance gene.New Phytol.214619–631. 10.1111/nph.14434

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 85

  ZhouS.CampbellT. G.StoneE. A.MackayT. F. C.AnholtR. R. H. (2012). Phenotypic plasticity of the Drosophila transcriptome.PLoS Genet.8:e1002593. 10.1371/journal.pgen.1002593

  • Pubmed Abstract
  • CrossRef
  • Google Scholar