Local adaptation contributes to gene expression divergence in maize.

Gene expression links genotypes to phenotypes, so identifying genes whose expression is shaped by selection will be important for understanding the traits and processes underlying local adaptation. However, detecting local adaptation for gene expression will require distinguishing between divergence due to selection and divergence due to genetic drift. Here, we adapt a QST-FST framework to detect local adaptation for transcriptome-wide gene expression levels in a population of diverse maize genotypes. We compare the number and types of selected genes across a wide range of maize populations and tissues, as well as selection on cold-response genes, drought-response genes, and coexpression clusters. We identify a number of genes whose expression levels are consistent with local adaptation and show that genes involved in stress response show enrichment for selection. Due to its history of intense selective breeding and domestication, maize evolution has long been of interest to researchers, and our study provides insight into the genes and processes important for in local adaptation of maize.

Introduction

Local adaptation occurs when different optimal trait values across environments lead to phenotypic differentiation among populations (Kawecki and Ebert 2004). Identifying locally adapted traits is important for animal and crop production (Howden ; Takeda and Matsuoka 2008), predicting response to climate change (Aitken ; Franks and Hoffmann 2012; Bay ), and conservation genetics (Funk ). One commonly used approach to identify local adaptation is , which tests for trait divergence (Q) that exceeds neutral expectations based on sequence divergence (F) (Prout and Barker 1993; Spitze 1993; Whitlock 2008). However, while previous work has used and related approaches to identify specific traits showing evidence of selection, we lack broad-scale systematic investigations into the number and types of traits that are locally adapted.

Gene expression is a useful model trait for systematically investigating the evolutionary forces shaping phenotypic variation: expression is quantitative, can be heritable, and variation in gene expression can contribute to phenotypic variation and adaptation (Oleksiak ; Gibson and Weir 2005; Gilad ; Rockman and Kruglyak 2006; Roelofs ; Whitehead and Crawford 2006; Groen ). has previously identified local adaptation for gene expression in Drosophila melanogaster and salmon (Roberge ; Kohn ) and a study has identified genes that showed relatively high or low Q in Populus tremula (Mähler ). Other studies have used an extension of developed by Ovaskainen to identify genes showing evidence of local adaptation in expression (Leder ; Ravindran ). In this study, we leverage next-generation sequencing data for expression and genetic variation to test for selection on expression of the entire transcriptome. In addition, we take advantage of a recent extension of that detects adaptation of continuous traits in large diversity panels that do not have clear subpopulations (Josephs ).

In this study, we investigate the role of local adaptation in shaping gene expression in the crop species Zea mays. Selection on gene expression has previously been shown to be important for maize evolution. For example, expression of the locus tb1 (Doebley ; Wang ) is responsible for the evolution of apical dominance during domestication and, transcriptome-wide expression divergence is prevalent between domesticated maize and its wild relative teosinte (Lemmon ). In addition, expression variation in domesticated maize is often associated with phenotype (Kremling ) . However, deleterious mutations are important contributors to expression variation in maize (Kremling ), implying that not all expression variation in maize is adaptive.

Here, we aim to understand the extent to which variation in gene expression in domesticated maize is driven by divergent selection caused by local adaptation and identify which genes show evidence of selection on their expression levels. We tested for selection using a published data set of 302 diverse maize lines each with RNAseq data from approximately 37,000 genes. We investigated enrichments of selective signals in genes that were differentially expressed in response to cold stress and drought, and selection on gene expression modules identified with coexpression network analyses taken from tissue-specific expression data. We detected selection on the expression of 60 unique genes across seven different tissue types and found an enrichment of drought-response genes among genes with the strongest signal of selection. Overall, these results show that local adaptation has shaped the expression of some genes and that this method has potential to identify specific genes and processes that are important for local adaptation.

Methods

Testing for selection on gene expression

Divergence between populations for a quantitative trait can be predicted by divergence at neutral genetic markers and additive genetic variation (V), assuming the trait evolves neutrally and the trait value is made up of an additive combination of allelic effects (Henderson 1950, 1953; Thompson 2008). If a sample does not have discrete populations, the genetic principal components (PCs) that explain most of the genetic variation can be used as a measure of divergence between populations and the other PCs can be used to estimate V. We briefly explain a test for selection using gene expression divergence measured across genetic PCs. More details on the test ( are available in Josephs .

Gene expression for a specific gene in M individuals is described by . If the gene expression levels described by evolve neutrally, we can describe the distribution of as follows: where μ is the mean expression value across individuals, V is the additive genetic variation for expression, and K is the kinship matrix of the individuals. The kinship matrix K can be decomposed so that, , where U is an matrix where the columns are eigenvectors of K and Λ is a diagonal matrix of corresponding eigenvalues. The eigenvectors of K are the genetic PCs of the population. We define as the mth eigenvector and λ as the mth eigenvalue. The amount of trait variation explained by the mth PC, standardized by how much neutral genetic variation is explained by that PC, is

Under neutrality, . If selection contributes to trait divergence along the mth PC, C may fall outside the neutral distribution. For this study, we tested the first five PCs for selection and the remaining PCs were used to estimate V. To test for selection, we use a test statistic (Q).

For a focal PC i,

Intuitively, these ratios of variances are similar to a standard measure of Q in that the numerator describes between-population expression-level variance and the denominator describes within-population expression-level variance. Genes with a high value of Q will have expression levels that are the most divergent at the between-population level compared to the neutral expectation. An important feature of Q is that environmental variation in phenotype will increase variation at the lower PCs used to estimate V and thus the amount of divergence expected due to drift. This property of the test means that environmental variation reduces our ability to detect selection.

Maize genomic and transcriptomic data

Expression and genotype data came from a subset of a maize diversity panel generated by Flint-Garcia . These lines represent the diversity present in public-sector maize-breeding programs worldwide. The dataset includes both temperate and tropical lines, as well as popcorn and sweet corn lines. In general, the temperate lines, along with sweet and popcorn lines, are more closely related to each other than they are to the tropical lines (Liu ). Temperate lines have less genetic variation than tropical lines, likely due to bottlenecks that occurred during breeding (Kremling ). The temperate group contains stiff stalk and nonstiff stalk lines, which represent the major heterotic groups used in breeding to create hybrids with heterosis (Liu ; Flint-Garcia ). Members of the stiff stalk group, which includes the reference line B73, are often used as the female parent for hybrids while members of the nonstiff stalk are often used as the male parents (Romay ).

Whole-genome sequence (Bukowski ) and RNAseq data for seven tissues (Kremling ) from plants grown in a common garden are available for these lines. The majority of lines had one RNAseq sample per tissue. For the subset of lines with more than one replicate, we randomly selected a single replicate to represent the line. Subsequent analysis only included genes that were expressed in all individuals for a given tissue type; which meant that we had between 8435 and 11,555 genes per tissue type (sample sizes listed in Supplementary Table S1).

We used 78,342 randomly chosen SNPs to create a kinship matrix for each tissue type, reflecting the slightly differing set of lines present for each tissue. We arranged and standardized each kinship matrix so that each cell, K of the n × n matrix is the genotypic covariance between the ith and jth lines following the procedure described in Josephs . After testing for selection as described above, FDR adjusted P-values were calculated to correct for multiple testing with the P.adjust function in R (Benjamini and Hochberg 1995; R Core Team 2020).

Cluster enrichment

We tested for local adaptation in the expression of gene coexpression modules. Walley profiled the transcriptome and proteome of 23 tissues spanning vegetative and reproductive stages of maize development using mRNA-seq and electrospray ionization tandem mass spectrometry. They then used weighted gene coexpression network analysis (WGCNA) to identify three gene expression networks consisting of 31,447 mRNA, 13,175 proteins, and 4267 phosphoproteins, respectively. They grouped genes with similar expression patterns into coexpression modules (“clusters”) using hierarchical clustering. Each cluster was assigned to the tissue(s) in which the cluster eigengene was most highly expressed. Their analysis resulted in 66 coexpression networks containing anywhere from 4 to 9574 genes. We calculated the median expression value for the genes in the 51 clusters that had more than 100 genes and used the same method outlined above on the median expression of each cluster to identify clusters that could be locally adapted.

Environmental response genes

We tested for enrichment of signals of selection in genes that show expression changes in response to cold and drought. Cold-response genes were identified by Avila , who estimated the transcript abundance in leaves of 22,000 genes in two Z. mays inbred lines (CG60 and CG102) during and after cold temperature exposure and identified 10,549 genes differentially expressed in response to cold exposure. Drought-response genes were identified by Forestan , who measured transcript abundance in young leaves of the inbred line B73 and calculated differential expression between well-watered and drought-stressed (10 days) treatments. Forestan identified 3181 differentially expressed genes (FDR < 0.01) and 28,983 nondifferentially expressed genes.

Drought-response genes had higher daytime expression level in leaves than genes that did not show drought response (Supplementary Figure S1). To ensure that overlaps between drought-response genes and selected genes were not due to both sets of genes being biased toward high expression genes, we chose a subsample of 3500 of the non-drought-response genes with high expression to use as a comparison set (Supplementary Figure S1). There was not a significant difference in daytime leaf expression level between cold-response and non-cold response genes, so we did not adjust the test for gene expression level.

With both datasets, we used a Fisher’s exact test to compare the proportion of genes that show evidence of selection (un-adjusted P-value <0.05) in environmental-response genes compared with other genes (see Supplementary Tables S5, S6, and S7 for sample sizes). We used the un-adjusted P-value so that we had enough genes in each category to use Fisher’s exact test. We only tested for enrichment in tissue-PC combinations that had evidence of at least one selected gene at FDR < 0.1. P-values were then adjusted for multiple testing using a Bonferroni correction (n = 15).

GO Enrichment Analysis

We tested subsets of genes identified as having signals of selection on gene expression for enrichment of GO biological process terms using the GO Enrichment Analysis tool on geneontology.org. (Ashburner ; Gene Ontology Consortium 2019; Mi ). We used the genes that went into our selection analysis for a given tissue as the reference list and the genes whose expression was under selection along a specific PC in that same tissue as the analyzed list. We used Fisher’s exact test and FDR as calculated by the Benjamini–Hochberg procedure for multiple testing correction as the settings for the enrichment analysis.

Data availability statement

All data used in this study were previously published. The RNA-seq data are from Kremling and available from https://datacommons.cyverse.org/browse/iplant/home/shared/commons_repo/curated/Kremling_Nature3RNASeq282_March2018. The genomic data are from Bukowski and available from https://datacommons.cyverse.org/browse/iplant/home/shared/commons_repo/curated/Qi_Sun_Zea_mays_haplotype_map_2018/282_onHmp321. Cold-response genes were identified by Avila and are available from Additional File 7 at https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-018-5134-7. Drought-response genes were identified by Forestan and are available in Supplementary Data S3 at https://doi.org/10.1111/pce.13660. Coexpression clusters were generated by Walley and are available in Supplementary Table S9 at https://science.sciencemag.org/content/353/6301/814.full. All codes used to produce this manuscript are available at https://jgblanc.github.io/Blancetal/. The kinship matrices can be found here: https://github.com/jgblanc/Blancetal/tree/master/data/Kinship_matrices.

Supplementary material available at figshare: https://doi.org/10.25387/g3.13172018.

Results

Detecting selection on expression of individual genes

We tested for selection on gene expression of 8435 to 11,555 genes in seven tissues for 109–239 genotypes (see Supplementary Table S1 for sample sizes), along the first five PCs within each tissue type. Note that because there were different genotypes sampled in each tissue type, the genetic PCs do not always correspond across tissues (Supplementary Figures S2, S3, and S4). Across all tissues, PC 1 separated out tropical from temperate genotypes and lower PCs separated stiff stalk from nonstiff stalk genotypes, popcorns from other genotypes, or separated out genotypes within the stiff stalk and/or nonstiff stalk subpopulations (Supplementary Figures S2, S3, and S4).

Sixty unique genes show evidence of expression divergence consistent with local adaptation along one of the first 5 PCs (FDR < 0.1, Figure 1A, Supplementary Table S2). We plot an example of the signal of selection on two genes to demonstrate what expression values look like when selection is inferred along a specific PC (Figure 1, B and C). There were five genes that had evidence for selection on expression in multiple tissues and/or multiple PCs. The PC-tissue combination with the most genes under selection was PC 5 in adult leaf expression measured during the day. Genes with divergence along PC 5 in adult leaf tissue are enriched for GO biological process terms cellulose catabolic process (FDR = 0.0323), plant-type cell wall biogenesis (FDR = 0.00853), and glucan biosynthetic process (FDR = 0.0287).

Figure 1

Signals of selection on gene expression in domesticated maize. (A) The number of genes where FDR < 0.1 in each of the seven tissues for the first five PCs. (B) PC 1 plotted against the mean-centered expression level of the gene GRMZM2G152686 as expressed in adult leaves during the day. Each point represents one maize genotype and is colored by subpopulation. The solid line shows the linear regression and the dashed lines show 95% confidence intervals of the neutral expectation. (C) Similar to plot (B) except PC 5 plotted against mean-centered expression of the gene GRMZM2G069762.

Selection on expression of coexpression clusters

Gene expression is often correlated across genes, so summarizing expression across coexpression clusters could improve power to detect selection (Kliebenstein 2020). With this in mind, we calculated median expression across previously identified coexpression modules (Walley ) and tested for selection on median gene expression for each module. However, none of the clusters showed evidence of selection (FDR > 0.1). The test with the strongest evidence of selection was the “Root Meristem” cluster, which showed evidence of selection along PC 5 in leaf adult tissue measured during the day (P = 2.4 × 10−4, FDR = 0.43). While the “Root Meristem” cluster had the highest expression in root meristems in Walley , many of these genes were still expressed in adult leaves in their study. Overall, these results suggest that coexpression clusters, as identified by correlations in expression within one genotype, are not broad targets of selection.

Selection on expression of environmental response genes

The spread of maize into North America required adaptation to different climatic factors (Swarts ), so we investigated selection specifically on genes that were differentially expressed in response to cold (Avila ) and in response to drought (Forestan ).

To test for evidence of selection on genes that were differentially expressed in response to cold, we compared selection signals in 12,239 genes that showed differential expression (FDR < 0.1) after either 1 or 4 days of cold treatment to 11,379 genes that did not show evidence of differential expression using data from Avila . We only investigated the 15 tissue-PC combinations where at least one gene showed significant evidence of selection at FDR < 0.01, looking at all genes under selection at an uncorrected level of P < 0.05 (Supplementary Tables S3 and S4). The strongest signal for enrichment was for daytime expression in adult leaf tissue along PC 5, where genes whose expression changed in response to cold were more likely to have evidence of local adaptation for expression (Bonferroni P = 0.06, Supplementary Table S5, Figure S5).

We found a significant enrichment of selection signals in 560 genes that showed decreased expression in response to drought in the B73 line compared to 3500 genes with similar leaf expression levels but that were not differentially expressed in drought (Supplementary Table S6). Specifically, expression in adult leaf tissue in both day and night showed evidence of enrichment for signals of selection along PC 5. Fourteen percent of genes downregulated in drought showed evidence of selection on leaf expression during day and night, while 8.1% of genes without drought response had evidence of selection for leaf expression during the day and 6.9% had evidence for selection on leaf expression at night (Bonferroni P = 0.00363 for day Bonferroni P = 1.635 × 10−5 for night) (Figure 2). The 328 genes that had increased expression in drought did not show any enrichment for selection (Figure 2, Supplementary Table S7).

Figure 2

Enrichment for signals of selection in genes downregulated in drought. The percentage of genes that show evidence of selection along PC 5 (P < 0.05) in adult leaf expression during the day and night for genes that are downregulated in drought, upregulated in drought, and show no change in response to drought. The asterisks show significance for enrichment compared to genes that were not drought-responsive: ** P < 0.01, *** P < 0.001.

Discussion

Systematically identifying genes important for local adaptation is crucial for understanding how local adaptation shapes trait variation. Here, we used an extension of to identify genes with expression divergence consistent with local adaptation in domesticated maize. Out of a dataset of expression of ∼10,000 genes measured across seven tissue types, we identified 60 genes with expression divergence consistent with local adaptation in at least one tissue type. Additionally, we found evidence that genes involved in drought response and cold response are enriched for signals of selection.

Our results contribute to a growing body of evidence that genetic variation for gene expression is shaped by selection. Previous studies in maize and other species have shown that rare variants affecting gene expression are often under negative selection (Josephs ; Kremling ; Glassberg ) and that there is weak stabilizing selection on gene expression levels in the field (Groen ). Alongside evidence for negative selection, and related analyses have demonstrated that local adaptation shapes between-population divergence in expression for some genes (Whitehead and Crawford 2006; Roberge ; Kohn ; Jueterbock ; Ravindran ). This is the first study to use Q, a -based method that detects selection on expression in the absence of clear subpopulations. With increasing availability of large transcriptomic studies conducted on diversity panels, methods for detecting selection on expression in the absence of clear subpopulations will be useful for understanding how selection shapes expression variation.

The enrichment of signals of adaptive divergence in genes involved in environmental response provides evidence for types of environmental factors that could contribute to adaptive divergence in expression. A number of pieces of evidence suggest that genes important for drought response had expression values shaped by local adaptation. There is an enrichment for signals of selection along PC 5 in genes that have decreased expression in response to experimental drought. One gene that shows adaptive expression divergence along PC 5 in leaf tissue (FDR = 0.02 for day and FDR = 0.01 for night) codes for the protein ZmRD22B, a putative maize RD22-like protein (Phillips and Ludidi 2017). RD22 proteins are thought to play a role in drought response through the ABA (abscisic acid) signaling pathway (Xu ) and ZmRD22B itself is predicted to localize to the cell wall and is upregulated in response to drought and exogenous ABA (Phillips and Ludidi 2017). Additionally, the group of genes we detected as having significant expression divergence along PC 5 in leaf tissue, including ZmRD22B, is enriched for GO biological processes cellulose catabolic process, plant-type cell wall biogenesis, and glucan biosynthetic process. In leaf tissue, PC 5 separated out individuals in the nonstiff stalk heterotic group of maize, suggesting that further investigations into gene expression and drought response in this subpopulation may be a promising future direction.

However, the link between genes important for stress response and evidence of local adaptation for gene expression in well-watered conditions is complex. The environmental response genes used in this study were identified from studies of differential expression in a few temperate maize genotypes. Stress-induced changes in gene expression could be beneficial responses that help the individual cope with stress or deleterious responses caused by the individual’s inability to maintain function in stressful conditions (Ghalambor ). If stress responses tend to be adaptive and improve function in the stressful condition, then local adaptation for expression in nonstressful conditions could reflect constitutive changes in expression in genotypes more likely to experience the stress. In contrast, if stress responses tend to be maladaptive in the stress environment, then local adaptation for expression in nonstressed environment could reflect further selection for reduced response even in nonstressful environments. For both cases, clearly understanding selection on the expression of environment-response genes will require additional experiments that measure expression changes in different environments across a diverse panel of genotypes.

While our method was successful in identifying genes whose expression is consistent with local adaptation, we only detected selection on 60 genes. Maize domestication and improvement have involved genome-wide selection (Wright ; Hufford ; Swarts ; Wang ), so we may expect to see evidence of selection on the expression on many more than 60 genes. There are a few potential explanations for why evidence of selection on gene expression may be limited. First, transcriptomes are a snapshot in a specific developmental time and environment and this study may have missed tissues, developmental time points, or environments in which expression has been under strong selection. Second, Q loses power when there is high environmental variation (V) for a trait. V increases trait variance explained by later (“within population”) PCs and, since these later PCs are used to generate a neutral expectation of divergence along focal PCs, high V will increase the amount of expression variation expected under neutrality (Supplementary Figure S6). This pattern means that high V will reduce power to detect selection (Josephs ) and the reduction in power due to V may be especially strong in expression data, which tends to be noisy and measured in few or no replicates.

An additional limitation of this study and the Q approach is that we were only able to investigate genes that were expressed in all individuals for a given tissue type. Q models phenotypes as additive combinations of allelic effects (Josephs ), and so the model is not robust to phenotypic distributions where a large number of individuals have a phenotype of 0. However, many of the expression changes that are important for phenotypic change may involve genes being turned on and off, not quantitative expression changes (Zhou ). In addition, maize has many presence–absence variants and the expression of these genes will appear to be 0 in individuals with the absent allele (Hirsch ; Zhou ). Methods to detect adaptive divergence in traits with non-normal distributions will be useful for future progress and may be able to detect more instances of adaptation.

Altogether, our work demonstrates that Q can be used to systematically detect genes whose expression is shaped by local adaptation and has shown its effectiveness in a large dataset from domesticated maize. We not only were able to detect selection on specific genes, but on combinations of genes based on environmental response patterns. Overall, our work shows that this method has potential for use in a number of large diversity panels while suggesting ways forward for better detecting selection on gene expression.

Click here for additional data file.