De novo transcriptome sequencing and gene co-expression reveal a genomic basis for drought sensitivity and evidence of a rapid local adaptation on Atlas cedar (Cedrus atlantica)

Abstract

Introduction: Understanding the adaptive capacity to current climate change of drought-sensitive tree species is mandatory, given their limited prospect of migration and adaptation as long-lived, sessile organisms. Knowledge about the molecular and eco-physiological mechanisms that control drought resilience is thus key, since water shortage appears as one of the main abiotic factors threatening forests ecosystems. However, our current background is scarce, especially in conifers, due to their huge and complex genomes.

Methods: Here we investigated the eco-physiological and transcriptomic basis of drought response of the climate change-threatened conifer Cedrus atlantica. We studied C. atlantica seedlings from two locations with contrasting drought conditions to investigate a local adaptation. Seedlings were subjected to experimental drought conditions, and were monitored at immediate (24 hours) and extended (20 days) times. In addition, post-drought recovery was investigated, depicting two contrasting responses in both locations (drought resilient and non-resilient). Single nucleotide polymorphisms (SNPs) were also studied to characterize the genomic basis of drought resilience and investigate a rapid local adaptation of C. atlantica.

Results: De novo transcriptome assembly was performed for the first time in this species, providing differences in gene expression between the immediate and extended treatments, as well as among the post-drought recovery phenotypes. Weighted gene co-expression network analysis showed a regulation of stomatal closing and photosynthetic activity during the immediate drought, consistent with an isohydric dynamic. During the extended drought, growth and flavonoid biosynthesis inhibition mechanisms prevailed, probably to increase root-to-shoot ratio and to limit the energy-intensive biosynthesis of secondary metabolites. Drought sensitive individuals failed in metabolism and photosynthesis regulation under drought stress, and in limiting secondary metabolite production. Moreover, genomic differences (SNPs) were found between drought resilient and sensitive seedlings, and between the two studied locations, which were mostly related to transposable elements.

Discussion: This work provides novel insights into the transcriptomic basis of drought response of C. atlantica, a set of candidate genes mechanistically involved in its drought sensitivity and evidence of a rapid local adaptation. Our results may help guide conservation programs for this threatened conifer, contribute to advance drought-resilience research and shed light on trees' adaptive potential to current climate change.

Keywords: Atlas cedar; RNA-Seq; adaptive capacity; climate change; conifers; drought sensitiveness; eco-physiology; phenotypic diversity.

Figure 1

Atlas cedar (Cedrus atlantica) trees at Fiñana, F (A) and Dornajo, D (C) sites (Sierra Nevada range, south Spain). Sampling location (latitude N, longitude W, and elevation m a.s.l.) and mean climate for the period 1971-2021 are indicated for F (B) and D (D) sites, respectively. The blue lines indicate the monthly total precipitation, the red lines indicate the monthly mean temperature, the shaded orange area indicates the water deficit period; total annual precipitation (P, mm) and mean annual temperature (T, °C) are also noted. Mature cones with spreading seeds (E) and germinating seedlings (F) were highly viable in both sites. Mean seed dry weight (G); the inset indicates the p value for ANOVA.

Figure 2

Drought treatment experimental design and sampling procedure. Net photosynthesis (A) and stomatal conductance (Gs) were monitored in controls (C), immediate drought (I, after 24 hours of drought treatment), and extended drought (E, after 20 days of drought treatment).The E group was divided into drought resilient seedlings (R) and non-resilient seedlings (N), based on recovery capacity after the extended drought. T, temperature; RH, relative humidity; VPD, vapor-pressure deficit; n indicates the number of individuals from Fiñana (F) and Dornajo (D), respectively (n=F+D).

Figure 3

Net photosynthesis, A (A) and stomatal conductance, Gs (B) measured in seedlings from Fiñana (F) and Dornajo (D) subjected to treatments: control (C), immediate drought (I, after 24 hours of drought treatment), and extended drought (E, after 20 days of drought treatment).The E group was divided into drought resilient seedlings (R) and non-resilient seedlings (N), based on recovery capacity after the extended drought. Groups are noted by the treatment/response (C, I, E; R, N) and population codes (F, D). Different letters indicate significant diferenteces (p<0.05) by ANOVA.

Figure 4

Principal Component Analysis (PCA) among the 17 normalized samples belonging to the four drought treatments/responses (Control C, immediate I, drought resilient R and non-resilient N) from Fiñana (F) and Dornajo (D) locations.

Figure 5

Dendrogram representing the clustering analysis of normalized samples (voom) subset for high variance with the objective of detecting outliers. Samples with a “D” after the number belong to Dornajo and with a F, to Fiñana. The last letters represent the treatment/response (C control, I immediate, R drought resilient, N non-resilient).

Figure 6

Hierarchical cluster dendrogram showing co-expression modules using weighted gene co-expression network analysis (WGCNA) of the counts to identify gene modules underlying drought stress at four different treatments and responses (C, I, R and N) in Fiñana (A) and Fiñana+Dornajo populations (B). (A) 466 modules corresponding to branches are labelled with colours indicated by the colour bands underneath the tree. With 0.25 threshold merging, 194 modules were generated. (B) 64 modules corresponding to branches are labelled with colours indicated by the colour bands underneath the tree. With 0.25 threshold merging, 15 modules were generated.

Figure 7

Heatmap showing the -log10(P-value) of the Fisher’s exact test calculated between the genes belonging to the 194 obtained modules in Fiñana (A) and to the 15 obtained modules in Fiñana + Dornajo samples (B), and the pairwise DE genes up and down-regulated (IvsC up, IvsC down, EvsI up, EvsI down, RvsN up, RvsN down). Given the low number of modules in Fiñana+Dornajo, it was possible to plot the number of overlapping genes in the heatmap (B). 20 modules were found to be statistically significant in Fiñana (A), and 8 modules in Fiñana + Dornajo (B).

Figure 8

Heatmaps showing gene expression levels of the genes within the modules across three treatments (C, I, E) and responses (R, N). Warm colours represent up-regulated genes whereas cold colours represent down-regulated genes. The six modules that showed expression patterns related to the treatment/response are represented for Fiñana (A), and the two modules related to the response for Fiñana + Dornajo (B).

Figure 9

GO term enrichment results for the hub genes of the selected modules visualized using REVIGO. (A) Aquamarine module, I vs C up comparison (Fiñana samples). (B) Lightblue3 module, I vs C down comparison (Fiñana). (C) Red module, E vs I up comparison (Fiñana). (D) Beige module, E vs I down comparison (Fiñana). (E) Paleturquoise module, R vs N up comparison (Fiñana). (F) Blue module, R vs N down comparison (Fiñana). (G) Turquoise module, R vs N up comparison (Fiñana + Dornajo). (H) Greenyellow module, R vs N down comparison (Fiñana + Dornajo).

Figure 10

GO term enrichment results for genes containing fixed SNPs between F and D locations (A) and R and N responses (B) visualized using REVIGO.

Figure 11

Principal component analysis (PCA) of the SNPs of the 17 samples from Fiñana and Dornajo populations.