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.