Rapid evolutionary changes in gene expression in response to climate fluctuations.

There is now abundant evidence of rapid evolution in natural populations, but the genetic mechanisms of these changes remain unclear. One possible route to rapid evolution is through changes in the expression of genes that influence traits under selection. We examined contemporary evolutionary gene expression changes in plant populations responding to environmental fluctuations. We compared genome-wide gene expression, using RNA-seq, in two populations of Brassica rapa collected over four time points between 1997 and 2014, during which precipitation in southern California fluctuated dramatically and phenotypic and genotypic changes occurred. By combining transcriptome profiling with the resurrection approach, we directly examined evolutionary changes in gene expression over time. For both populations, we found a substantial number of differentially expressed genes between generations, indicating rapid evolution in the expression of many genes. Using existing gene annotations, we found that many changes occurred in genes involved in regulating stress responses and flowering time. These appeared related to the fluctuations in precipitation and were potentially adaptive. However, the evolutionary changes in gene expression differed across generations within and between populations, indicating largely independent evolutionary trajectories across populations and over time. Our study provides strong evidence for rapid evolution in gene expression, and indicates that changes in gene expression can be one mechanism of rapid evolutionary responses to selection episodes. This study also illustrates that combining resurrection studies with transcriptomics is a powerful approach for investigating evolutionary changes at the gene regulatory level, and will provide new insights into the genetic basis of contemporary evolution.

INTRODUCTION

We now have abundant evidence for rapid phenotypic evolution in contemporary natural populations (Franks, Weber, & Aitken, 2014; Parmesan, 2006; Thompson, 2013). Thanks to recent work in genomics, we also know the genetic basis of many ecologically important adaptive traits (Anderson, Willis, & Mitchell‐Olds, 2011; Ehrenreich & Purugganan, 2006; Hughes, 1999; Stinchcombe & Hoekstra, 2007). However, our understanding of the genetic basis of contemporary evolution remains relatively modest (Ehrenreich & Purugganan, 2006; Franks & Hoffmann, 2012). Phenotypic traits can rapidly evolve through shifts in the frequency of alleles that affect coding sequences and/or patterns of gene expression (Hancock & Rienzo, 2008; Romero, Ruvinsky, & Gilad, 2012). Prior research has examined changes in allele frequencies by using genomic information to look for signatures of selection (Bakker, Toomajian, Kreitman, & Bergelson, 2006; Burke et al., 2010; Guo et al., 2018; Turner, Bourne, Von Wettberg, Hu, & Nuzhdin, 2010), or by directly comparing allele frequencies in ancestral and descendant populations (Burke et al., 2010; Franks, Kane, O'Hara, Tittes, & Rest, 2016). However, information on the role of gene expression changes in rapid contemporary evolution remains scarce (McCairns & Bernatchez, 2010).

We know that changes in gene expression can alter phenotypes and play a key role in adaptive divergence among species (Fay & Wittkopp, 2008; King & Wilson, 1975; Wray, 2003). The same principles may be applied to intraspecific population adaptation (López‐Maury, Marguerat, & Bähler, 2008). For example, phenotypic variation in a trait under selection could be due to variation in gene expression, so that selection on the trait results in rapid evolutionary changes in gene expression that underlie the adaptive phenotypic change. However, it remains particularly challenging to link changes in gene expression with adaptive evolution, in part because gene expression levels differ among cell types, across developmental stages, in response to numerous environmental factors, and because populations often present pronounced levels of genetic variation that influences gene expression (Babbitt, Haygood, Nielsen, & Wray, 2017). Yet, a growing number of studies have established a direct causal link between phenotypic evolution and changes in gene expression. Some of these studies have linked contemporary evolution to changes in expression of a single gene (Abzhanov, 2004; Doebley, 2004; Shapiro et al., 2004). For example, studies have shown that the famous evolutionary changes in beak morphology in Darwin's finches are due to evolutionary changes in expression during development of the gene Bmp4 (Abzhanov, 2004). Additionally, a modest number of studies have used transcriptomics to observe genome‐wide variation in gene expression and to assess evolutionary changes in expression patterns. For example, the transcriptome of guppies whose ancestors were transplanted from streams with cichlid predators to streams without predators underwent regulatory changes, and after selection, reflected that of the native predator‐free population, suggesting rapid adaptive evolution of gene expression (Ghalambor et al., 2015). In a study in which expression changes were examined over time, Campbell‐Staton et al. (2017) compared phenotypes and transcriptomes in populations of green anole lizards sampled before and after an extreme winter storm and found rapid adaptive changes in gene expression linked to increased cold tolerance. One experimental evolution study examined regulatory changes between Escherichia coli lineages evolved in a glucose‐limited medium for 20,000 generations. This study identified parallel regulatory changes in 59 genes (Cooper, Rozen, & Lenski, 2003), demonstrating not only that evolution in gene expression can occur, but that in this case, the changes were parallel among replicated populations. However, the degree of parallelism in evolution in general remains a central question in evolutionary biology (Blount, Lenski, & Losos, 2018; Grant & Grant, 2002; Stern & Orgogozo, 2009). Such studies provide unique insights into the genetic basis of contemporary evolution and suggest that rapid evolution can potentially occur through changes in gene expression. However, we generally still lack direct evidence that natural selection can lead to adaptation through contemporary evolutionary changes in gene expression.

A particularly powerful approach to studying contemporary evolution of traits, genes or gene expression is the “resurrection approach”, where ancestral and descendant generations from natural populations are grown together under common conditions (Franks, Sim, & Weis, 2007). For organisms with seeds or other propagules that can be stored, this approach allows for a direct assessment of genetically‐based evolutionary change over time (Franks, Hamann, & Weis, 2018). Such resurrection studies can be combined with increasingly accessible genomics tools to compare DNA and/or RNA sequences between ancestors and descendants to uncover the genetic basis and mechanisms contributing to adaptive responses to climate change (Franks et al., 2018; Franks & Hoffmann, 2012). In this approach, differences in allele frequencies, levels of gene expression, or phenotypic traits between ancestors and descendants provide direct evidence of evolutionary changes in these characteristics (Franks et al., 2018).

One particularly well‐studied case of evolutionary responses to fluctuations in climatic conditions in natural populations is the evolution of Californian Brassica rapa populations in response to drought (Franks, 2011; Franks et al., 2007, 2015, 2016; Franks & Weis, 2008; Hamann, Weis, & Franks, 2018). Previous studies in this system demonstrated that variation in flowering time has a genetic basis and is heritable (Franks et al., 2007), is under directional selection by drought (Weis, Wadgymar, Sekor, & Franks, 2014), and evolves rapidly following dry (Franks, 2011; Franks et al., 2007) and wet periods (Hamann et al., 2018). Individual flowering time candidate genes differed in expression levels between early and late flowering individuals (Franks et al., 2015). At the genome scale, many genes showed evolutionary shifts in frequencies between ancestral and descendant lines. Yet, few of these genes showed similarities between two geographically distinct populations, suggesting independent evolutionary trajectories rather than parallel evolutionary change in these populations (Franks et al., 2016). This work has started to shed light on the complexity of the genetic basis underlying phenotypic changes in B. rapa following climatic fluctuations. However, we do not yet know to what extent evolutionary changes in gene expression might underlie adaptive phenotypic evolution, either in this system or more broadly (Franks & Hoffmann, 2012; Richards et al., 2009; Todd, Black, & Gemmell, 2016).

In this study, we used seeds collected from two Californian B. rapa populations over 18 years of fluctuating precipitation patterns (Hamann et al., 2018). We compared these populations, resurrected under common conditions, in genome‐wide analyses of gene expression using RNA sequencing (RNA‐seq) to investigate whether evolutionary changes in gene expression contribute to rapid adaptive evolution. Specifically, the study aimed to determine: (a) evolutionary changes in gene expression; (b) whether those changes were adaptive; and (c) if evolutionary changes in expression were similar between populations and over different drought episodes. We assessed evolutionary changes in expression by comparing expression profiles of ancestors and descendants grown under common conditions and determining which genes were significantly differentially expressed between ancestors and descendants. To determine if the changes in expression were adaptive, we examined the function of the differentially expressed genes using existing annotations, with the hypothesis that differentially expressed genes would be linked to functions such as stress response and flowering, which are likely to be adaptive in the context of drought. Finally, we examined changes in gene expression over the two drought episodes to determine if changes in expression were repeated over time, and we compared evolutionary changes between the populations to determine if these changes were parallel at the two sites.

MATERIALS AND METHODS

Study species

Brassica rapa (L.) Brassicaceae, commonly known as field mustard, is an annual, herbaceous plant. It comprises numerous economically important varieties of vegetable and oil crops cultivated around the world. Its genome has been sequenced and annotated (Wang et al., 2011; Zhang et al., 2018).

In coastal California, B. rapa is a winter annual, where germination is triggered by the winter rains and the growing season is terminated by the onset of annual summer droughts (Franke et al., 2006; Franks et al., 2007). Two populations, Arboretum (ARB) and Back Bay (BB), located in Orange County, California, have been extensively studied for their adaptive phenotypic changes in response to drought episodes (Franks, 2011; Franks et al., 2007, 2008; Hamann et al., 2018), and were used for this current study. These populations are located within a 3 km radius around the UC Irvine campus and experience similar climatic conditions. However, the BB population grows on a more coastal site, with sandier and consistently drier soil relative to the more variable ARB site (Franke et al., 2006; Franks, 2011; Hamann et al., 2018).

Plant material and experimental setup

The detailed seed sampling and full experimental design has been described previously (Hamann et al., 2018). Briefly, we made use of a unique set of resources consisting of seeds from ancestral and descendant lines from two B. rapa populations (ARB and BB) collected before and after two consecutive drought episodes (drought 1:1999–2004, drought 2:2012–2014, with below‐average precipitation relative to patterns over a 100‐year period, Table S1). The first baseline seed sampling was done in 1997 after a series of years with normal precipitation patterns. Seeds were then collected in 2004, after five years with below‐average precipitation during the germination‐ and growing‐season (Franke et al., 2006; Franks et al., 2007). A series of years with fluctuating precipitation followed, and seeds were again collected in 2011 after two intermediate wet seasons. Finally, seeds were sampled in 2014 after three record‐breaking drought seasons. In total, seeds were collected from eight accessions (four time points for two populations) over an 18‐year period: ARB'97, BB'97, ARB'04, BB'04, ARB'11, BB'11, ARB'14, and BB'14.

In the greenhouse, a refresher generation was grown from 100 randomly selected seeds from each accession following best practices for resurrection studies (Franks et al., 2018). Sixty maternal lines were then randomly selected from each accession, from which two seeds were grown under well‐watered conditions (8 accessions × 60 maternal lines × 2 seeds = 960 plants). At 30 days after germination, leaf tissue was sampled from 6 randomly selected biological samples among a subset of maternal lines per accession for whole transcriptome analysis by RNA‐seq (6 biological replicates × 4 generations × 2 populations = 48 samples). At this time point, all of the plants had started flowering, so they were at the same developmental stage. The tissue type and collection time was previously identified as optimal for comparative analysis of gene expression in B. rapa (Franks et al., 2015). Fresh leaf material was sampled from fully developed basal leaves, and the leaves collected were similar among all plants. The leaf collection took place during the course of 1 hr (10.00–11.00 hr, 4 hr after dawn), and each sample was collected in chilled 2 ml safe‐lock Eppendorf tubes, immediately flash‐frozen in liquid nitrogen, and stored at –80°C.

RNA extraction and sequencing

Leaf material was finely ground in liquid nitrogen and total RNA was extracted using the RNeasy Plant Mini kit (Qiagen). The integrity and quality of the RNA was confirmed with a NanoDrop 1000 spectrophotometer. Only one sample was excluded as extraction was insufficient, thus we retained only five biological samples (instead of six) for the BB’04 accession. For library preparation, we used 0.4 μg of total RNA per sample. Library preparation was conducted at the Center for Genomics and Systems Biology, New York University, using the Agilent Bravo Automated Liquid‐Handling Platform and following their NEBNext Ultra RNA‐adapted protocol. The NEBNext Ultra RNA Library Kit for Illumina (Cat. No. E7530L) was used to perform cDNA barcoded library preparation with NEBNext Multiplex Oligos for Illumina Index Primers (Cat. No. E7335L). Fragment size estimates were confirmed on an Agilent TapeStation using High Sensitivity D1000 ScreenTape, and final concentrations of sequencing libraries were quantified by qPCR using the KAPA Library Quantification Kit (Kapa Biosystems). The final cDNA libraries were sequenced using the Illumina HiSeq 2500 v.3 high‐output system, with a 1 × 100 bp configuration at the Center for Genomics and Systems Biology, New York University. RNA‐seq libraries were initially checked using MultiQC (Ewels, Magnusson, Lundin, & Käller, 2016). On average, 12 million reads were sequenced per sample, and all bases had a quality score above 30 (Phred 33), indicating that the sequencing was reliable.

De novo assembly, annotation and alignment

Reads were processed with trimmomatic v.0.36 to remove the adapter sequences, and low‐quality bases (Bolger, Lohse, & Usadel, 2014). Using trinity v.2.5.1, a reference transcriptome was assembled de novo (Haas et al., 2013) from the combined data of the sample with the largest library size from each of the accessions (see details in S2). Although a reference genome is available for a domesticated accession of B. rapa (Cheng et al., 2011; Wang et al., 2011), we used a de novo transcriptome assembly to preserve the diversity of the wild Californian populations.

To annotate this newly assembled reference, we used two Universal Protein Resource (UniProt) reference protein databases from Brassica and the closely related Arabidopsis (Brassicaceae), and identified homologous proteins encoded by the transcripts. Using blast v.2.5.1 (Camacho et al., 2009), and the “add_blastx_hit_to_trinity_id.pl” script, we sequentially added annotations, first from Brassica and then from Arabidopsis, to find the most similar proteins. Next, the BLAST output was filtered by amino acid length (>50 aa), percent identity (>70%), and e‐value (cutoff < 1e‐20) to annotate the transcripts with the best‐matching UniprotIDs.

To estimate gene expression levels for each individual per condition, we ran Trinity's “align_and_estimate_abundance.pl” script (see details in S2) with RSEM estimation (Li & Dewey, 2011; Li et al., 2009) using bowtie 2 v.2.3.4.1 (Langmead & Salzberg, 2012). To generate un‐normalized, as well as TPM‐ and TMM‐normalized count tables for each de novo‐assembled transcript, we used the Trinity “abundance_estimates_to_matrix.pl” script, specifying the gene‐to‐transcript map previously generated.

Evolutionary changes in gene expression

In our resurrection experimental design, evolutionary changes in gene expression are indicated by statistically significant differences in expression between ancestors and descendants. In other words, for any genes that show significant differences in expression between ancestors and descendants, there is evidence for evolutionary changes in expression of those genes. These significant differences in gene expression between ancestors and descendants were identified as differentially expressed genes (DEGs) by comparing generations using edgeR implemented in Trinity (Love, Anders, & Huber, 2014; Robinson, McCarthy, & Smyth, 2010; Zhou, Xia, & Wright, 2011). We used Trinity's default filters to remove lowly expressed genes before analysing both populations (ARB and BB) separately. Exact‐Tests were performed to compare gene expression between generations within populations using the “run_DE_analysis.pl” script in Trinity. A cutoff of 1% FDR was used to identify DEGs. In addition, we extracted the strongest DEGs for each comparison (0.1% FDR cutoff and 4‐fold expression change) by running the “analyse_diff_expr.pl”. Finally, we partitioned these genes into clusters according to their patterns of differential expression across samples using hierarchically clustered gene trees cut at 60% as recommended in the “define_clusters_by_cutting_tree.pl” Trinity script (Haas et al., 2013).

After filtering out lowly expressed genes, 29,647 genes were retained for analysis, out of which 17,419 (58.7%) were annotated with UniProt. Populations were analysed separately, since previous phenotypic and genotypic studies revealed pronounced divergence between ARB and BB populations (Franks et al., 2016; Hamann et al., 2018). Specifically, we analysed changes in gene expression in four pairwise generation comparisons. First, we compared the 2004 post‐drought generation to the 1997 predrought generation (’04 vs. ’97), representing the initial transition from normal wet years to abnormally dry seasons. We then compared the 2011 generation, collected after two intermediate wet years, to the 2004 post‐drought generation (’11 vs. ’04). We further examined the 2014 post‐drought generation, collected during a record‐breaking drought year, relative to the intermediate 2011 generation (’14 vs. ’11). These chronological comparisons of generations stemming from before and after successive drought episodes allowed us to evaluate whether both drought episodes induced similar changes in gene expression that may be reversed during the intermediate wet seasons. Finally, in addition to assessing successive changes in gene expression in response to fluctuating precipitation patterns, we considered long‐term changes between the most recent post‐drought generation from 2014 and the earliest predrought generation from 1997 (’14 vs. ’97).

For each comparison, we recorded the number of significant DEGs: the number of up‐ and down‐regulated genes, and the proportion of annotated genes. Venn diagrams were drawn using venny 2.1 (Oliveros, 2015) to visualize genes unique to or shared evolutionary changes in gene expression between generation comparisons within and across populations.

To determine if the evolutionary changes in gene expression were adaptive, we explored the functional annotations of all significant DEGs that were annotated. While many genes lacked annotations, there was still a substantial sample of DEGs with functional information. To assess gene function, we examined gene ontology (GO) terms for biological processes assigned from the UniProt database, and GO terms were grouped under functional categories using david 6.8 (Huang, Sherman, & Lempicki, 2009). We hypothesized that if the evolutionary changes in gene expression were adaptive, then the functional categories would include biological processes related to stress response and flowering, which are linked to drought stress.

RESULTS

A large number of genes showed significant differential expression across generations of plants grown under common conditions, indicating rapid evolutionary changes in gene expression accompanying fluctuations in precipitation levels (Figure 1). The number of DEGs (genes showing evolutionary changes in expression) differed between populations and pairwise generation comparisons. For the ARB population, the number of genes showing evolutionary changes in expression was close to a thousand when the earliest (1997) and most recent (2014) generations were compared (Figure 1a). When comparing ARB generations chronologically, the largest proportion of evolutionary changes in expression occurred over the first five‐year drought episode between 1997 and 2004 (Figure 1a). For the BB population, there were fewer than 100 genes that showed significant evolutionary changes in expression between 1997 and 2014, and 139 genes evolved in expression between 2011 and 2014, the second drought episode, while fewer genes evolved in response to the first drought episode (Figure 1b).

FIGURE 1

Number of DEGs identified between generations (at 1% FDR) in (a) ARB and (b) BB for 2004 vs. 1997 ('04 vs. '97), 2011 vs. 2004 ('11 vs. '04), 2014 vs. 2011 ('14 vs. '11), and the long‐term comparison 2014 vs. 1997 ('14 vs. '97). Upregulated genes are shown in orange (above the zero line) and downregulated genes in blue (below the zero‐line), with the hatched portion of the bars representing the percentage of annotated genes. The total number of up‐ and downregulated genes are detailed at the top of the bars and the percentage of annotated genes is indicated in parenthesis. On the x‐axis, the years in red (2004, 2014) represent drought periods, while the years in blue (1997, 2011) represent wet periods. The total number of genes analysed for each comparison differed depending on filtering for lowly expressed genes. Specifically, our analyses included 41,577 genes for ARB ’04 vs. ’97, 43,562 genes for ARB ’11 vs. ’04, 43,384 for ARB ’14 vs. ’11, 41,674 for ARB ’14 vs. ’97, 42,578 for BB ’04 vs. ’97, 43,388 for BB ’11 vs. ’04, 43,290 for BB ’14 vs. ’11, and 42,744 for BB ’14 vs. ’97

Regulatory evolution of genes involved in drought stress responses and flowering

The proportion of annotated genes showing evolutionary changes in expression (DEGs) ranged from 7% to 48% depending on the pairwise generation comparison (Figure 1). While less than half of DEGs were confidently annotated, we had annotations for more than 600 DEGs, representing a substantial sample of the total number of DEGs. Of those DEGs that had annotations, many had functions that appeared to be related to the fluctuations in precipitation patterns in this system. Many of the genes showing evolutionary changes in expression had annotations involving GO terms related to stress responses, including drought stress, as well the regulation of flowering time, which is a key trait involved in drought escape (Figure 2; Tables S3–S8).

FIGURE 2

Functional classification of annotated DEGs identified for the long‐term generation comparison (2014 vs. 1997) in (a) ARB and (b) BB into enriched GO terms for biological processes associated with stress responses (in dark red), regulation of flowering time (in yellow), and other functions (in blue). Asterisks indicate the five GO categories found in both populations (all others are unique to one population). For ARB, out of the 865 DEGs that were identified between 2014 and 1997, 514 were annotated, and 305 were grouped under GO terms. For BB, out of the 84 DEGs that were identified between 2014 and 1997, 29 were annotated, and 22 were grouped under GO terms

For ARB (2014 vs. 1997), we found evolutionary regulatory changes in genes involved in the response to water deprivation (Figure 2a, Table S5), including genes encoding AQUAPORINS, which facilitate the transport of water (Table S11). Numerous genes were also related to phytohormonal signalling associated with stress responses such as signalling through abscisic acid, salicylic acid, jasmonic acid, and gibberellin (Figure 2a, Table S5). Several MYB‐transcription factor‐encoding genes were represented in these functional categories (Table S11), which were consistently downregulated in the 2014 post‐drought generation. Multiple genes were related to circadian rhythm, photoperiodism and flowering (Figure 2a; Tables S5 and S11), including genes coding for the transcription factors Circadian clock associated 1 (CCA1) and GATA8, as well as the zinc finger proteins CONSTANS‐LIKE 5 and 9 (COL5 and COL9).

For BB (2014 vs. 1997), we found genes related to cell wall differentiation and responses to the hormones ethylene, auxin and brassinosteroid (Figure 2b; Table S8), including ones encoding MYB‐related transcription factors and Heat shock proteins (Table S10). Two genes associated with the term circadian rhythm, including one encoding for GATA8, were upregulated in the 2014 post‐drought generation (Table S12).

Common evolutionary changes in gene expression across populations and drought events

To identify common expression patterns across time points within and between populations, we used two additional comparisons. Within each population, we examined whether systematic changes in gene expression could be observed across generations in response to fluctuating precipitation patterns and consecutive drought events. Additionally, we examined whether common evolutionary changes in gene expression occurred across generations in both populations, which could point towards parallel evolution between populations.

Systematic changes in gene expression in response to consecutive drought events

We aimed to determine if consecutive drought events caused repeated evolutionary changes in gene expression. To do so, we assessed the proportion of DEGs that were in common in different generation comparisons within populations (Figure 3), and examined whether the direction of gene regulation was consistent across drought events. Across these comparisons, the percentage of genes showing common changes in expression across different generations and drought episodes was 1%–8% in ARB (Figure 3a), and 0.5%–12% in BB (Figure 3b). The highest number of systematic changes was found between ARB'04 versus ’97 (196 DEGs) and ARB'14 versus ’97 (865 DEGs). Here, we identified 86 common genes (out of a total of 1,086 DEGs) that showed similar expression patterns in the two post‐drought generations relative to the initial predrought generation (Figure 3a). Interestingly, 12 of these common DEGs belonged to one of ten identified gene clusters, where genes were downregulated after the first drought episode and stayed downregulated thereafter, even during the intermediate wet season (Figure 4a). Among this cluster, we found three genes involved in the regulation of flowering time that were consistently downregulated relative to the initial predrought generation. Eleven other genes were grouped into another cluster, where gene expression was upregulated after the first drought episode and stayed upregulated thereafter. The remaining gene clusters had less regular expression patterns across generations (but see complete heatmap in S9). Across BB generation comparisons, we found the highest number of common genes between BB'11 versus ’04 (116 DEGs) and BB'14 versus ’11 (139 DEGs), with 30 common genes out of a total of 255 DEGs (Figure 3b). Here, most genes were unannotated. Six gene clusters, with an average of 18 genes, were identified as having similar expression profiles across generations. In contrast to ARB, BB expression profiles indicate a pronounced change in gene expression during the intermediate wet season followed by a reversal to the initial state during the second drought episode (see complete heatmap in S10). One gene cluster in particular reflects this regulation pattern, where eight genes shared similar expression profiles in both post‐drought generations (Figure 4b).

FIGURE 3

Comparison of common and unique evolutionary changes across time within each population. Number of DEGs at 1% FDR, and percentages in parentheses, between pairwise generation comparisons (i.e., ’04 vs. ’97, ’11 vs. ’04, ’14 vs. ’11, and ’14 vs. ’97) within populations for (a) ARB and (b) BB

FIGURE 4

Expression patterns of one representative gene cluster across generations in (a) ARB and (b) BB. The gene cluster represented for ARB includes 12 genes, and 8 for BB. The expression of single genes within each cluster are plotted in grey, in addition to their mean expression profile in blue. Gene expression is shown as log2‐transformed, median‐centred reads per kilobase of target transcript length per million reads (RPKM). Biological samples on the x‐axis are clustered based on expression similarities of all highly significant DEGs (0.1% FDR cutoff and 4‐fold expression change) identified across generations within population (see full heatmaps for both generation in S9–10). Predrought generations are shown in blue (1997 and 2011), and post‐drought generations are shown in red (2004 and 2014)

Parallel evolutionary changes across populations

We examined whether common evolutionary changes in gene expression occurred across generations in both populations to examine the degree of parallel evolutionary change between populations. We found that 1%–3% of genes showing evolutionary changes in expression evolved in parallel in the two populations (Figure 5a–d). The highest number of parallel evolutionary changes between populations was found for the long‐term comparison between ARB'14 versus ’97 (865 DEGs) and BB'14 versus ’97 (84 DEGs; Figure 5d). Out of this total of 949 DEGs, 24 genes evolved in parallel across populations, and 15 of those were annotated. The annotated genes that showed common evolutionary changes in both populations were associated with major plant hormonal pathways, defence response and cell differentiation (Figure 5e). Particularly, one GATA and two MYB transcription factor‐encoding genes were commonly downregulated in the 2014 post‐drought generations of both populations (Figure 5f, Table S13).

FIGURE 5

Common and unique evolutionary changes between the ARB and BB populations across generations (DEGs at 1% FDR) for (a–d) all generation comparisons; (e) functional classification of the 15 annotated (out of the 24 total) common DEGs between ARB'14 vs. ’97 and BB'14 vs. ’97 into GO terms for biological processes associated with stress response (in red), regulation growth and flowering time (in yellow) and other processes (in blue); and (f) common evolutionary changes in expression profiles of transcription factors found in the ARB and BB populations across predrought (1997) and post‐drought generations (2014)

DISCUSSION

This study directly documented rapid evolutionary changes in gene expression in response to contemporary shifts in climatic conditions. By combining a resurrection approach using four generations of Brassica rapa collected before and after two major drought episodes in California with transcriptomic comparisons, we identified a large number of genes that showed differences in expression patterns between ancestors and descendants grown under common conditions, indicating evolutionary changes in gene expression. This finding suggests that gene expression can rapidly evolve at a contemporary time scale, and that changes in expression may underlie some cases of phenotypic evolution. Many of the genes were not annotated, so it is not possible to determine the functional characteristics of every gene that evolved differential expression. However, a substantial fraction of the genes was annotated, giving us a sizable sample of the transcriptomic changes that evolved to examine. Many of the annotated genes that showed evolutionary changes in expression were associated with functions related to drought stress responses and the regulation of flowering time, a trait key to drought escape (Franks, 2011; Franks et al., 2007), suggesting that evolutionary changes in gene expression probably contributed to adaptation. Moreover, few common changes in gene expression were found across generations and populations, indicating that evolutionary changes were largely independent. Here, we discuss our results in light of the potential importance of regulatory changes in contemporary evolution.

Regulation of gene expression contributes to rapid contemporary evolution

Based on the reasoning that evolutionary modifications of gene expression can be the basis for morphological divergence between species (King & Wilson, 1975), it follows that regulatory evolution can also be a mechanism of rapid adaptation within species in response to changing environments. Gene expression regulation is a flexible mechanism for adjusting rapidly to local environmental conditions, but it can also lead to long‐term evolutionary responses (López‐Maury et al., 2008; Romero et al., 2012). In fact, changes in gene expression may be one of the fastest routes to rapid evolution without the need for novel mutations in coding regions (Mäkinen, Papakostas, Vøllestad, Leder, & Primmer, 2016), or for changes in allele frequencies at important coding regions (Hancock & Rienzo, 2008). Instead, existing genes within a population may undergo changes in expression levels in response to changes in environmental conditions both in the short‐term and for long‐term evolutionary adaptation (López‐Maury et al., 2008), yet, this has rarely been explicitly studied. Our unique combination of a resurrection approach with comparative transcriptomics allowed for the direct detection of evolutionary changes in gene expression in response to changes in environmental conditions.

Although the resurrection approach we used provides strong evidence that differences in gene expression of the plants from different generations are based on evolutionary changes, it is also possible that these differences are due to variation in developmental stage or tissue composition at the time of collection (Hodgins‐Davis & Townsend, 2009). However, the plants in this study were at the same stage (right after flowering initiation), and we collected only one tissue type. Even for the same tissue type and developmental stage, genes could differ in expression levels because of differences in timing. For example, on the day of collection, some genes might show greater levels of expression in plants that had flowered more recently compared to plants that had flowered several days earlier. However, looking for such differences was exactly the point of our study. We interpret such findings not as due to error of collecting plants that differ in developmental timing, but rather as uncovering expression level differences linked with phenotypes that may be under selection. Furthermore, by raising four generations collected before and after drought episodes under common conditions, we controlled for environmental effects and minimized manipulation and perturbation before and during tissue sampling to extract the clearest possible signal for gene expression differences across generations (Hodgins‐Davis & Townsend, 2009). Because we are comparing expression levels in ancestors and descendants raised under common conditions, we can confidently state that significant changes in gene expression we observed between generations are probably evidence of evolutionary changes in expression. We suggest that the uncovered changes in gene expression are the result of phenotypic selection and represent the genetic basis of rapid adaptive evolution previously documented in phenotypic traits. Ultimately, by documenting rapid changes in gene expression across generations, our study revealed that regulatory changes are a contributing mechanism of rapid, contemporary phenotypic evolution. To our knowledge, only one other study has shown rapid, intraspecific evolutionary changes in genome‐wide gene expression following climatic shifts (Campbell‐Staton et al., 2017). Our results, along with other recent work (Campbell‐Staton et al., 2017; Cooper et al., 2003; Ghalambor et al., 2015), provide evidence that regulatory changes can rapidly occur in response to strong selective events, and that changes in gene expression can drive rapid phenotypic evolution at contemporary time scales.

Adaptive potential of regulatory changes in B. rapa's drought stress response

Because the resurrection approach that we used does not in and of itself determine the cause of evolution (Franks et al., 2018), the evolutionary shifts in gene expression documented here could potentially be due to either drift or selection. We suggest that the evolutionary changes in gene expression we observed were unlikely to be driven by genetic drift alone, and that instead many of these changes may have been caused by natural selection on phenotypic traits and be adaptive, for several reasons. Given the relatively short evolutionary timespan (17 generations or fewer), the outcrossing mode of reproduction of this species, and the relatively large population sizes of this weedy plant (>1,000 individuals per population based on field observations), we expect that selection would generally be strong relative to drift. Simulation modeling we conducted based on this data indicated that in the most cases, the evolutionary changes in gene expression we observed could not be explained by drift and sampling alone (unpublished data). We analysed a large number of genes, and it is certainly possible that genetic drift played a role in evolutionary changes at some loci. However, we also expect that the evolutionary changes in gene expression we observed were not due to drift alone and may have been driven by selection in many cases because many of the genes observed to evolve in expression had annotations related to drought stress responses and flowering time regulation (as described below), providing further evidence that many regulatory changes were probably adaptive rather than caused by random drift alone.

We also found evidence that the observed evolutionary changes in gene expression may have been adaptive based on the annotation of the genes that were differentially expressed. Because of the fact that the populations experienced repeated drought episodes, and that prior studies showed that traits such as flowering time were important in drought responses (Franks, 2011; Franks et al., 2007), we hypothesized that the DEGs may have had such biological functions as stress response and flowering phenology. Many of the genes showing evolutionary changes in expression did not have annotations, so we do not have further information on these genes, though they might be interesting for future study. The remaining genes represent a sample of the full transcriptome that we can use to functionally characterize the evolutionary changes in expression. It is important to point out that this assessment was based on existing annotations, and the DEGs were not functionally validated, which was beyond the scope of this study. But this assessment revealed that many of the genes that showed evolutionary changes in expression had functions that appeared related to drought response, including the regulation of stress response and flowering time, making them good candidates for adaptive changes affecting plant fitness, and supporting the hypothesis that some of the evolutionary changes in expression were driven by natural selection.

Under water‐deficit, plants can respond via physiological and biochemical responses, including stomatal closure, repression of cell growth, photosynthesis and activation of respiration, as well as via cellular and molecular processes, such as osmolyte accumulation (Shinozaki & Yamaguchi‐Shinozaki, 2007). Additionally, the phytohormone abscisic acid plays a critical role for plant adaptation to adverse abiotic conditions, including drought stress, by integrating various stress signals and controlling downstream responses (Shinozaki & Yamaguchi‐Shinozaki, 2007). In our study, we identified evolutionary changes in functional and regulatory genes related to responses to water deprivation and abscisic acid, the oxidation‐reduction process, cellular water homeostasis, and stomatal control, suggesting functions in the initial stress response, as well as in establishing drought resistance in B. rapa. Notably, in the ARB population, regulatory changes occurred in several genes coding for Aquaporins, which control plant‐water relations (Afzal, Howton, Sun, & Mukhtar, 2016), suggesting their important role in this population's drought stress response. In both populations, we also found regulatory changes of several MYB transcription factors associated with the regulation of stomatal movement, the production of cuticular waxes, and with the regulation of flower development (Baldoni, Genga, & Cominelli, 2015). Their systematic downregulation in post‐drought generations of both populations suggests their important role in B. rapa's evolutionary response to drought.

Our study also aimed at detecting regulatory changes in genes linked to flowering time, to explain the rapid evolutionary advances documented in B. rapa post‐drought generations (Franks et al., 2007; Hamann et al., 2018). Flowering is triggered by internal and external signals and controlled by four genetically defined pathways, including over 100 genes that code for transcription factors, receptors and other proteins (Bäurle & Dean, 2006; Engelmann & Purugganan, 2006; Flowers, Hanzawa, Hall, Moore, & Purugganan, 2009; Srikanth & Schmid, 2011). Certain proteins play key roles in the flowering time network; for example, the MADS‐box transcription factor FLOWERING LOCUS C (FLC) acts to suppress flowering (Yan, Liang, Liu, & Zheng, 2010), while suppressor of overexpression of CONSTANS 1 (SOC1) and CONSTANS (CO) act as flowering promoters (Srikanth & Schmid, 2011). In our transcriptomic comparisons of pre‐ and post‐drought generations, we identified several DEGs with functions related to gibberellin signalling, circadian rhythm, photomorphogenesis, flower development and flowering time, confirming that evolutionary changes occurred along the flowering time gene network. Most noteworthy, two genes encoding CONSTANS‐like protein flowering promoters were upregulated in the post‐drought descendant generation of the ARB population, which could be responsible for their advanced flowering time relative to predrought ancestors (Hamann et al., 2018). In both populations, regulatory changes occurred in the transcription factor GATA8, which can promote the expression of SOC1 (Richter, Bastakis, & Schwechheimer, 2013). The close association to SOC1, and the common regulation in both populations, suggests the adaptive role of GATA8 in the regulation of flowering time. Finally, in BB populations we identified downregulation of AGL70, coding for a flowering repressing FLC‐like protein, which may be responsible for the consistently earlier flowering time in BB relative to ARB populations (Hamann et al., 2018).

While our study identified regulatory changes in genes with important functions under drought stress or for the regulation of flowering time, it is important to note that many DEGs could not be confidently associated with biological functions and much of the regulatory changes need further investigation. Additionally, it remains challenging to assertively link gene expression to fitness to confirm the adaptive nature of variation in gene expression, and to examine the type of selection that occurred at the regulatory level (Signor & Nuzhdin, 2018; Whitehead & Crawford, 2006). Recent studies are beginning to explore the relationships among genome‐wide gene expression, phenotypic traits and fitness to better understand the genetic and phenotypic basis of adaptive responses to selection (Ayroles et al., 2009; Groen et al., 2020).

Although many of the evolutionary changes in gene expression observed may have been adaptive, the mechanisms behind the evolution of gene expression were not determined in this study. Gene expression is regulated both by cis‐elements near the expressed genes, including promotors and enhancers, and by trans‐factors such as transcription factors that are encoded elsewhere in the genome (Emerson et al., 2010; Wittkopp, Haerum, & Clark, 2004). Our study focused on gene expression patterns, so we could not directly determine if the observed changes in expression were due to changes in cis‐ or trans‐factors. However, we would expect that if the changes are due to cis‐factors, we might find a correspondence between the genes that changed in expression in this study and the genes that were previously identified to show evolutionary changes in allele frequency (Franks et al., 2016). In contrast to this expectation, less than 1% of the genes that showed changes in gene expression had also previously been identified as F ST outliers corresponding to loci with allele frequency changes in response to the first drought episode (1997 vs. 2004; Franks et al., 2016). This result suggests that cis‐regulatory loci may not have been as important as trans‐factors in causing the expression changes that resulted in phenotypic evolution. However, this tentative conclusion should be regarded with substantial caution because cis‐ and trans‐factors were not compared directly. But if this is the case, this finding would be in line with a growing body of studies showing that many changes in gene expression are the result of selection on trans‐acting loci, mostly with small, but cumulatively large, effects (Groen et al., 2020), that changes in trans‐acting factors often make a larger contribution to intraspecific changes (Romero et al., 2012; Signor & Nuzhdin, 2018; Wittkopp et al., 2004; Wittkopp, Haerum, & Clark, 2008), and that trans‐acting regulation may be more important for environment‐dependent differences in gene expression (Signor & Nuzhdin, 2018).

Independence of evolutionary changes in gene expression over time and across populations

By comparing evolutionary changes in the regulatory profiles of two populations across consecutive drought events, our experimental design allowed investigating whether evolutionary changes in gene expression occurred recurrently over drought episodes and in parallel across populations, or whether regulatory changes were largely independent.

While both populations showed parallel phenotypic evolutionary changes (i.e., post‐drought advanced flowering time; Franks et al., 2007; Hamann et al., 2018), our transcriptomic results suggest that the genetic basis to these changes differed between the ARB and BB populations, and that population divergence was maintained even after the second drought episode. The natural populations in this study experienced repeated episodes of drought, and responses of phenotypic traits such as flowering time showed consistent patterns, including the evolution of earlier flowering after each drought and of later flowering following wet periods (Hamann et al., 2018). We therefore predicted that changes in gene expression might show similarly consistent patterns in response to climatic fluctuations. However, very few regulatory changes occurred repeatedly across fluctuating precipitation (Figures 3 and 4). In the ARB population, more regulatory changes evolved in response to the first drought episode than in response to the second. In contrast, fewer regulatory changes were selected for over the first drought episode in the BB population. These different responses to selective events may be related to the fact that the two populations were probably at different distances from the phenotypic optimum because of past adaptation to local soil conditions. As the ARB population was initially adapted to generally wetter and more variable soil moisture conditions (Franks, 2011), it is likely that the first drought episode imposed stronger selection on the ARB population, which we detected at the regulatory level and previously at the phenotypic level (Hamann et al., 2018). Conversely, the BB population grows on sandier and more drained soil, and thus under consistently drier conditions (Franks & Weis, 2008). Previous studies found that BB plants generally flower earlier than ARB plants, even in predrought generations, suggesting that the BB population is adapted to drier conditions (Franks, 2011; Hamann et al., 2018). Our new transcriptomic data support this hypothesis, as we found comparatively fewer regulatory changes across BB generations in response to drought episodes.

Finally, we examined whether certain changes in gene expression evolved in both populations in response to drought episodes to assess the degree of parallel evolution between populations. A prior study found that changes in allele frequencies after the first 5‐year drought event were largely independent across these two populations (Franks et al., 2016). At the regulatory level, our results corroborate this prior genomic assessment, as we found very few common evolutionary changes between generations across populations. Only a small number of genes evolved in parallel in both populations between 1997 and 2014. These genes mainly coded for transcription factors such as MYBs, and the parallel evolutionary changes observed across populations suggest their important role in B. rapa's drought stress response. However, these parallel evolutionary changes represent less than 3% of the otherwise predominantly independent regulatory changes.

Taken together, our comparisons of expression profiles between generations among and between populations yielded very little commonality overall, indicating that evolutionary changes in gene expression are largely independent over time and across populations. In contrast with other studies that have found parallel genetic responses in populations responding to similar selection pressures (Jones et al., 2012; Wood, Burke, & Rieseberg, 2005), our results are in keeping with studies showing a role for contingency and a lack of parallel evolutionary changes at the genetic level (Grant & Grant, 2002; Nachman, Hoekstra, & D'Agostino, 2003; Roda et al., 2013).

To conclude, by combining transcriptomic profiling with the resurrection approach, our study revealed rapid evolutionary changes in gene expression between ancestors and descendants of B. rapa populations that experienced consecutive drought episodes in California. Many of the genes showing regulatory changes were grouped in functional categories related to stress responses and flowering time regulation, consistent with an adaptive response to drought. However, the evolved regulatory changes were not systematic across consecutive drought episodes within populations, and the genetic basis of B. rapa's rapid adaptation to drought differed across populations, indicating rather independent evolutionary trajectories across time and between populations. Future studies should explore the relationships among genome‐wide gene expression, phenotypic traits and fitness to better understand the genetic and phenotypic basis of adaptive responses to selection.

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

AUTHOR CONTRIBUTIONS

E.H., and S.J.F. conceived and designed the study. E.H. performed the greenhouse experiment; E.H., Z.J.‐L., and S.C.G. performed laboratory work; E.H., and C.P. analysed the data with guidance from N.C.K., and J.S.R; E.H., and S.J.F wrote the manuscript and S.C.G., Z.J.L., and C.P. edited the manuscript.

Supplementary Material

Click here for additional data file.

Table S11‐13

Click here for additional data file.