Analyzing Oak Transcriptomic Responses to Elevation

Cleaning and mapping were performed as described by [115 (link)]. Briefly, we first removed low-quality reads and retained sequences with a mean quality value above 20. We mapped reads onto the 25,808 oak gene models published with the reference oak genome [116 (link)], using BWA (V.0.6.1) with the default parameters. We then identified differentially expressed genes (DEG) with the DESeq2 package, using a P-value < 0.01 after adjustment for multiple testing with a false discovery rate (FDR) of 5%. We also considered two-fold changes in expression ratio as a threshold for identifying the genes with the highest degree of differential expression. We investigated the effects of dormancy stage, elevation and their interaction in likelihood ratio tests. The dormancy and elevation effects were assessed by comparing a statistical model without interaction (M1) with two reduced models for the dormancy (M2) and elevation (M3) effects, respectively. $$\mathrm{M}1:{\mathrm{Y}}_{\mathrm{ijk}}=\mathrm{\mu }+{\mathrm{D}}_{\mathrm{i}}+{\mathrm{E}}_{\mathrm{j}}+{\mathrm{\epsilon }}_{\mathrm{ijk}}$$
where D_i is the dormancy stage (i = ”Endodormancy” or “Ecodormancy”), and E_j the elevation effect (j = “low-100 m”, “medium-800 m” or “high-1,600 m” elevation). $$\mathrm{M}2:{\mathrm{Y}}_{\mathrm{jk}}=\mathrm{\mu }+{\mathrm{E}}_{\mathrm{j}}+{\mathrm{\epsilon }}_{\mathrm{jk}}$$ $$\mathrm{M}3:{\mathrm{Y}}_{\mathrm{jk}}=\mathrm{\mu }+{\mathrm{D}}_{\mathrm{i}}+{\mathrm{\epsilon }}_{\mathrm{jk}}$$
For estimation of the interaction effect, we compared the following a complete model: $${\mathrm{Y}}_{\mathrm{ijk}}=\mathrm{\mu }+{\mathrm{D}}_{\mathrm{i}}+{\mathrm{E}}_{\mathrm{j}}+{\left(\mathrm{D},,\ast ,\mathrm{E}\right)}_{\mathrm{ij}}+{\mathrm{\epsilon }}_{\mathrm{ijk}}\left(\mathrm{M},4\right)\mathrm{to}\mathrm{M}1$$
Three gene sets were thus generated: (i) DEGs (geneset#1) corresponding to dormancy regulation (regardless of elevation), (ii) DEGs (geneset#2) corresponding to differences in elevation (regardless of dormancy stage), and (ii) DEG (geneset#3) displaying a significant dormancy-by-elevation interaction. Annotations for each DEG were recovered from the published pedunculate oak genome sequence [116 (link)]. Below, we focus particularly on geneset#2 and #3, which should include the key molecular players involved in response to temperature (and potentially to local adaptation) to temperature for geneset#2, and should reveal differences in the strategies of oak stands analysed across bud phenological stages to cope with temperature variation for geneset#3. The elevation term used here does not allow to disentangle the effect of phenotypic plasticity from those of genetic differentiation. Indeed, [117 ] reported that the adaptive response of populations to divergent selection pressures along the elevation systematically results from a combination of non-optimal phenotypic plasticity and genetic differentiation. Thus, other experimental design such as reciprocal transplantations are needed to separate these two effects in the stands analyzed.
DEGs from geneset#2 and #3 were analyzed with EXPANDER software [118 (link)], which clusters genes according to their expression profile, using a Kmeans algorithm [118 (link)]. For both gene sets, we set the number of clusters to 5 (k = 5) to maximize the homogeneity of each cluster. The genes from each cluster were then used for independent gene set and subnetwork enrichment analysis (see below).