Assessing co-expression module level signatures of adaptive evolution
Local adaptation is often influenced by trait covariation (O’Brienet al ., 2019), thus modules formed by sets of strongly correlated transcripts effectively capture the architecture of local adaptation. We used the estimated genetic trait values for maternal trees from eqn. 1 to construct garden-specific co-expression networks using WGCNA v1.70 (Langdelfer & Horvath, 2008) in R. This is different from the standard practice of constructing co-expression networks using raw expression values across tissues, individuals or environments. Since we are using genetic values and not raw expression levels the co-expression network can be treated as a proxy for variance-covariance matrix (G -matrix), which has a long history of use in understanding the genetic architecture of local adaptation (summarized by MacPhersonet al ., 2015). Each co-expression network consists of modules formed by sets of strongly correlated transcripts. This modularity permits selection to act on only certain functional categories. Modules formed within each network were visualised using the prefuse force-directed edge-weighted layout in Cytoscape v.3.8.2 (Shannonet al ., 2003). We filtered out weakly connected transcripts and edges with weights below the 75th quantile for the module for clarity (Fig. S3).
Transcripts classified into each module were used to perform GO enrichment analyses using the same approach as described above for the three Q ST categories of interest. To address our hypotheses concerning adaptive evolution (H1 & H2) at the module level, we evaluated among-population differences in the co-expression module using a linear model with the first eigengene of each module as the response variable and population of origin as the predictor variable. We also assessed the relative representation (RR ) of individual adaptive transcripts (identified using the Q ST-F ST approach above, Fig. 2) in each module. For modules detected at the cold garden, we used the Cold-condA as theQ ST category, while for modules detected at the warm garden we used the Warm-condA category. In all cases, RR was estimated as:
\(RR\ \ =\ \frac{\frac{N_{Qst-module}}{N_{\text{Qst}}}}{\frac{N_{\text{module}}}{N_{\text{network}}}}\ \), (3)
where NQ st -module is the number of transcripts in a module that were also classified under the respective Q ST category,NQ st is the number of transcripts classified under the Q ST category (Fig. 2),N module is the total number of transcripts in a module and N network is the total number of transcripts used in WGCNA. Statistically significant over- or under-represented modules (two tailed test, p < 0.05) were identified by permuting the transcript labels and calculating 10,000 null RR estimates for each module. For significant modules, we evaluated whether RR was associated with module size using a Kruskal-Wallis test as implemented with the stats v. 4.0.1 package in R.