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 et al. 2007;Takeda and Matsuoka 2008), predicting response to climate change (Aitken et al. 2008;Franks and Hoffmann 2012;Bay et al. 2017), and conservation genetics (Funk et al. 2012). One commonly used approach to identify local adaptation is Q ST À F ST , which tests for trait divergence (Q ST ) that exceeds neutral expectations based on sequence divergence (F ST ) (Prout and Barker 1993;Spitze 1993;Whitlock 2008). However, while previous work has used Q ST À F ST 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 et al. 2002;Gibson and Weir 2005;Gilad et al. 2006;Rockman and Kruglyak 2006;Roelofs et al. 2006;Whitehead and Crawford 2006;Groen et al. 2020). Q ST À F ST has previously identified local adaptation for gene expression in Drosophila melanogaster and salmon (Roberge et al. 2007;Kohn et al. 2008) and a study has identified genes that showed relatively high or low Q ST in Populus tremula (Ma ¨hler et al. 2017). Other studies have used an extension of Q ST À F ST developed by Ovaskainen et al. (2011) to identify genes showing evidence of local adaptation in expression (Leder et al. 2015;Ravindran et al. 2019). 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 Q ST À F ST that detects adaptation of continuous traits in large diversity panels that do not have clear subpopulations (Josephs et al. 2019).
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 et al. 1997;Wang et al. 1999) 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 et al. 2014). In addition, expression variation in domesticated maize is often associated with phenotype (Kremling et al. 2019) . However, deleterious mutations are important contributors to expression variation in maize (Kremling et al. 2018), 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.
Divergence between populations for a quantitative trait can be predicted by divergence at neutral genetic markers and additive genetic variation (V A ), assuming the trait evolves neutrally and the trait value is made up of an additive combination of allelic effects (Henderson 1950(Henderson , 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 A . We briefly explain a test for selection using gene expression divergence measured across genetic PCs. More details on the test (Q PC Þ are available in Josephs et al. (2019).
Gene expression for a specific gene in M individuals is described by Z ¼ ½Z 1 ; Z 2 ; . . . Z m¼M . If the gene expression levels described by Z evolve neutrally, we can describe the distribution of Z as follows: Z $ MVNðl; V A KÞ;
(1) where l is the mean expression value across individuals, V A 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, K ¼ UKU T , where U is an n x n matrix where the columns are eigenvectors of K and K is a diagonal matrix of corresponding eigenvalues. The eigenvectors of K are the genetic PCs of the population. We define Ũm as the m th eigenvector and k m as the m th eigenvalue. The amount of trait variation explained by the m th PC, standardized by how much neutral genetic variation is explained by that PC, is
Under neutrality, C m $ Nð0; V A Þ. If selection contributes to trait divergence along the m th PC, C m 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 A . To test for selection, we use a test statistic (Q PC ).
For a focal PC i,
Intuitively, these ratios of variances are similar to a standard measure of Q ST in that the numerator describes betweenpopulation expression-level variance and the denominator describes within-population expression-level variance. Genes with a high value of Q PC will have expression levels that are the most divergent at the between-population level compared to the neutral expectation. An important feature of Q PC is that environmental variation in phenotype will increase variation at the lower PCs used to estimate V A 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.
Expression and genotype data came from a subset of a maize diversity panel generated by Flint-Garcia et al. (2005). 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 et al. 2003). Temperate lines have less genetic variation than tropical lines, likely due to bottlenecks that occurred during breeding (Kremling et al. 2018). 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 et al. 2003;Flint-Garcia et al. 2005). 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 et al. 2013).
Whole-genome sequence (Bukowski et al. 2018) and RNAseq data for seven tissues (Kremling et al. 2018) 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 ij of the n  n matrix is the genotypic covariance between the i th and j ith lines following the procedure described in Josephs et al. (2019). 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).
We tested for local adaptation in the expression of gene coexpression modules. Walley et al. (2016) 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.
We tested for enrichment of signals of selection in genes that show expression changes in response to cold and drought. Coldresponse genes were identified by Avila et al. (2018), 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 et al. (2020), who measured transcript abundance in young leaves of the inbred line B73 and calculated differential expression between well-watered and droughtstressed (10 days) treatments. Forestan et al. (2020) 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).
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 et al. 2000;Gene Ontology Consortium 2019;Mi et al. 2019). 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. Forestan et al. (2020)
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, andS4). 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 andC). 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).
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 et al. 2016) 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 et al. (2016), 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.
The spread of maize into North America required adaptation to different climatic factors (Swarts et al. 2017), so we investigated selection specifically on genes that were differentially expressed in response to cold (Avila et al. 2018) and in response to drought (Forestan et al. 2020).
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 et al. (2018). 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 andS4). 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).
Systematically identifying genes important for local adaptation is crucial for understanding how local adaptation shapes trait variation. Here, we used an extension of Q ST À F ST 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 et al. 2015;Kremling et al. 2018;Glassberg et al. 2019) and that there is weak stabilizing selection on gene expression levels in the field (Groen et al. 2020). Alongside evidence for negative selection, Q ST À F ST and related analyses have demonstrated that local adaptation shapes between-population divergence in expression for some genes (Whitehead and Crawford 2006;Roberge et al. 2007;Kohn et al. 2008;Jueterbock et al. 2016;Ravindran et al. 2019). This is the first study to use Q PC , a Q ST À F ST -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 et al. 2010) 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. Stressinduced 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 et al. 2007). 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 et al. 2005;Hufford et al. 2012;Swarts et al. 2017;Wang et al. 2020), 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 PC loses power when there is high environmental variation (V E ) for a trait. V E 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 E will increase the amount of expression variation expected under neutrality (Supplementary Figure S6). This pattern means that high V E will reduce power to detect selection (Josephs et al. 2019) and the reduction in power due to V E 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 PC approach is that we were only able to investigate genes that were expressed in all individuals for a given tissue type. Q PC models phenotypes as additive combinations of allelic effects (Josephs et al. 2019), 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 et al. 2020). 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 et al. 2014;Zhou et al. 2019). 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 PC 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.
J. Blanc et al. | 3 Downloaded from https://academic.oup.com/g3journal/article/11/2/jkab004/6114460 by guest on 04 August 2021