A drought stress experiment was carried out using Atlas cedar seedlings germinated from seeds sampled in two afforested populations of southern Spain (Fiñana, F and Dornajo, D, respectively; Figures 1A, C ). The reforestation project was performed with C. atlantica individuals from natural stands of northern Morocco during the last century. The study sites show contrasting climate conditions due to a longitudinal gradient, where aridity increases eastward, and due to elevation, where aridity increases downward. F site (Almería), is characterized by higher mean temperature and lower total precipitation, which determined increased summer drought, compared to D site ( Figures 1B, D ). During the period 2011-2015, air temperature (T) was registered hourly in both locations using a Hobo H8 data logger. Data loggers were placed inside open-bottom PVC cylinders situated at 1.5m from the ground and covered by sealing foam to prevent direct heating from the sun. Further, to obtain long-term climate data, monthly total precipitation (P) and T were downloaded from the EOBS database v23.1e for the period 1971–2021 (Haylock et al., 2008) using the Climate Explorer webpage (https://climexp.knmi.nl/). Local data of monthly total precipitation were also used to account for elevation gradients.

Seeds were obtained in the field from mature cones ( Figure 1E ) and germinated ( Figure 1F ) in the common garden of the University Pablo de Olavide (Seville, Spain). Both populations showed similar germination percentage (80.1% and 82.1%, for F and D, respectively), while seed dry weight was significantly higher in D site ( Figure 1G ). Seedlings grown in 0.77-liter pots with 38.5 cm² of surface exposed to evaporation, filled with sterile sand and perlite substrate in 3:1 ratio for about one year in the common garden. Environmental conditions in the common garden (temperature, relative humidity and PAR) were monitored using Onset Computing HOBO Pro v2 data loggers (HOBO U23-series). Mean annual air temperature was 17.9°C, 27.4°C in the hottest month (August), and 10.3°C in the coldest month, (January). Mean annual relative humidity was 68.3%, 55.9% in the drier month (July), and 81.2% in the wettest month (November). Mean PAR was 445.8 lux. Irrigation was provided to saturation. The seedlings were fertilized in pots with solid NPK synthetic fertilizer, containing 12.0% of total Nitrogen (N), of which 6.5% was nitric Nitrogen (N-NO₃), and 5.5% was ammoniacal Nitrogen (N-NH₄); Phosphorus pentoxide (P₂O₅) soluble in neutral ammonium citrate and water 12.0%, of which 6.0% was Phosphorus pentoxide (P₂O₅) water soluble; Sulfur trioxide (SO₃) total 15.0%, of which 12.0% was Sulfur trioxide (SO₃) water soluble; Magnesium Oxide (MgO) Total 2.0%, of which 1.5% was Magnesium oxide (MgO) water soluble; Potassium oxide (K₂O) water soluble 17.0%, Boron (B) total 0.02%, and Zinc (Zn) total 0.01%. Then, a subset of seedlings were randomly selected and introduced in controlled growth chambers (CGC) under constant conditions of humidity (70% relative humidity RH) and temperature (27°C). CGC air vapour pressure deficit (VPD, kPa) was calculated using temperature and relative humidity in the following equations (Sanchez-Salguero et al., 2015):

Saturated vapour pressure (svp) was calculated as: s v p = 0.66 * e ( 0.06 * T )

Then, VPD was obtained as: V P D = ( ( 100 − R H ) 100 ) * s v p where RH is the CGC air relative moisture (%), and T is the CGC air temperature (°C).

Independent of the following drought stress treatment, all seedlings stayed 15 days in CGC for acclimation, until the variables used to monitor control seedlings hydric status (net photosynthesis and stomatal conductance) showed high and stable values. Seedlings (both controls and those subjected to drought treatments) were monitored weekly by means of a gas exchange system provided with an infrared gas analyzer (IRGA) (LI-COR, LI-6400XT, EEUU). During this acclimation period, watering was carried out by adding 10 mL of water per day. After the acclimation, the variables used to monitor control seedlings (A and Gs) were relatively steady and averaged (± standard error) (A = 6.5 ± 0.5 μmol CO2·m^-2·s^-1, Gs = 110 ± 13 mmol H₂O cm^-2 s^-1). Photoperiod was stablished between 7:00 AM and 19:00 PM (12 hours light/12 hours darkness, light at 600 lux), while physiological measurements and further sampling for RNA extraction were performed between 10:00 AM and 1:00 PM, again, for both controls and seedlings subjected to the drought treatments.

Drought treatment ( Figure 2 ) was applied to 38 seedlings by reducing the RH of the CGC up to 30%, which is equivalent to a VPD of 2.3 kPa, whereas 22 remained as control (RH=70%), which is equivalent to a VPD of 1.0 kPa. Both CGC maintained 27 ° C of temperature and 10 mL per day of watering. Immediate drought response, extended drought response and drought resilience were assessed using gas exchange monitoring and gene expression by RNA extraction ( Figure 2 ). This experimental design allows a quantification of drought effects at a molecular resolution for seedling subjected to experimental drought. Drought resilience was defined here as the capacity of a seedling to reach gas exchange rates (A and Gs) similar to those prior to drought treatment. Defined this way, drought resilience encompasses the capacity to cope with drought, that is, resistance, and the ability to return to prior gas exchange levels after drought, that is recovery (Lloret et al., 2011). We defined three sample types within the experiment ( Figure 2 ): (i) Controls (C, n=11 from Dornajo and n=11 from Fiñana), seedlings were raised in CGC at 70% RH; (ii) immediate drought response treatment (I, n=8 from Dornajo and n=13 from Fiñana), seedlings were raised in CGC at 30% RH and needles for RNA extraction were collected 24 hours after; and (iii) extended drought response treatment (E, n=9 from Dornajo and n=8 from Fiñana), seedlings were raised in CGC at 30% RH and needles were collected for RNA extraction 20 days after the drought treatment started ( Figure 2 ).

The recovering drought response treatment introduces the concept of drought resilience (Lloret et al., 2011). E treatment seedlings were returned to control conditions after 20 days of drought treatment ( Figure 2 ). Then, recovery (regarding A and Gs values observed after drought) was monitored to quantify drought resilience (sensu Lloret et al., 2011). According to that, we defined two groups: drought resilient seedlings (R, n=4 from Fiñana; n=5 from Dornajo), defined as those E seedlings that reached, at least, the 30% of their A and Gs values in control conditions in less than 3 days; and non-resilient seedlings (N, n=4 from Fiñana, n=5 from Dornajo), defined as those E seedlings that died or did not manage to reach, at least, the 30% of their A and Gs values in control conditions in less than 3 days ( Figure 2 ). We set this threshold (30%) based on our own experience on previous essays, where we determined that, overall, seedlings that are able to recover about 1/3 of their initial rates are thereafter able to survive after the drought conditions end. All collected needles were stored in a freezer at -80°C until RNA extraction was performed.

Four individuals belonging to each treatment/response (C, I, R, N) from each location (Fiñana and Dornajo) were selected for the RNA extraction to analyze the same quantity of biological replicates. Spectrum Plant Total RNA Kit (Sigma-Aldrich, Saint Louis, MO, USA) was used to carry out a total RNA extraction from the collected needles. The manufacturer’s protocol was followed with one modification. 2% polyvinylpyrrolidone (PVP) K30 (Sigma-Aldrich) was added to the lysis buffer to enhance the extraction efficiency. Bioanalyzer 2100 (Agilent Technologies, Inc, Santa Clara, CA, USA) was used to test the quality of the extracted RNA. A ND-1000 Spectrophotometer (NanoDrop™, Thermo Fisher Scientific, Waltham, MA, USA) was used to measure its concentration. 1.8 ≤ OD260/280 ≤ 2.2, OD260/230 ≥ 1.8 and RNA Integrity Number (RIN) ≥ 7 were applied as quality standards. Then, all samples were sent to the National Center for Genomic Regulation (CNAG-CRG, Barcelona, Spain) for sequencing. There, quality control (QC) was performed by fluorescent-based quantification (RiboGreen, Invitrogen™, Thermo Fisher Scientific) and integrity evaluation (2100 Bioanalyzer, Agilent Technologies, Santa Clara, CA, USA), followed by a cDNA library preparation (Illumina Stranded mRNA Prep) and RNA-sequencing by Illumina HiSeq2500 (Illumina, Inc, San Diego, CA, USA).

To confirm the results from the RNA-seq analysis, some DEGs were selected, and their expression levels were confirmed by qRT-PCR analysis. For the quantitative real-time PCR (qRT-PCR), a TB Green Advantage qPCR Premix (Takara, Japan) was used and run on a QuantStudio 7 Flex Real-Time PCR System (Thermofisher, USA). For each gene, the measurements were performed in four replicates, and the average cycle thresholds (Ct) were used to determine the fold-change. The expression ratio (Log2FC) change was calculated by the 2-ΔΔCT method. 18S rRNA gene was used as a reference gene to normalize the qRT-PCR data. All of the primers used for the qRT-PCR are listed in Table S1 .

Net photosynthesis (A) and stomatal conductance (Gs) showed non-significant differences between populations, while drought treatment yielded a significant gas exchange decline ( Figure 3 ; Supplementary Material, Table S2 ). During the post-drought recovery, R individuals showed significantly higher values of A and Gs, compared to N individuals, in both populations ( Figure 3 , Table S2 ; see ANOVA results at Appendix 1).

We obtained statistically significant DEGs during immediate treatment when analyzing only Fiñana, and Fiñana and Dornajo individuals together ( Table 1 ). Nevertheless, given the lack of gene modules showing an expression pattern of interest when analyzing both populations ( Figure 8B ), we will focus the discussion of the C. atlantica seedling transcriptional profile after 24 hours of drought on Fiñana population.

Despite the low number of significantly up and down regulated genes in the two immediate drought gene modules ( Figure 8A ), they showed biologically meaningful functions ( Figure 9A, B ). C. atlantica is an isohydric tree regarding its physiological response to drought, characterized by its sensitivity to the atmospheric-moisture demand, which leads to a water-saving strategy by means of stomatal closure (Spinelli et al., 2011; Sanchez-Salguero et al., 2015; Lechuga et al., 2019). This stomatal closure is accomplished by increasing abscisic acid (ABA) concentrations (Moran et al., 2017). Consequently, it leads to a pronounced reduction of stomatal conductance (Gs) and limited carbon uptake by photosynthesis (A). This isohydric pattern has been observed in this work in C. atlantica seedlings through a progressive decrease of A and Gs values in both populations ( Figure 3 ). These findings are similar to those obtained in A. pinsapo seedlings using the same study design (Cobo-Simón et al., 2022), whose physiological response to drought is also isohydric.

During the 24 hours of drought treatment (I), C. atlantica seedling transcriptional response showed consistencies with an isohydric profile. Among the hub genes, we observed several genes related to chloroplast activity and photosynthesis among the up (isoamylase 3 chloroplastic isoform X1, high chlorophyll fluorescent 107) and down-regulated genes (probable plastid-lipid-associated 8, chloroplastic, soluble inorganic pyrophosphatase 6, chloroplastic-like), including transcription factors (transcription termination factor MTERF4, chloroplastic) ( Table S4 ), pointing to a regulation of photosynthetic activity.

Hub genes related with protein transport, protein kinase activity and protein phosphorylation pathways were observed among the up- (probable inactive serine threonine- kinase scy1, receptor kinase TMK1-like, GPI ethanolamine phosphate transferase 1, coatomer subunit delta, Tetratricopeptide repeat-containing domain) ( Table S4 ) and down-regulated genes (leucine-rich repeat receptor-like serine threonine tyrosine-kinase, Table S5 ). Protein phosphorylation has been associated with drought-priming heat stress tolerance in plants (Zhang et al., 2020). Interestingly, protein kinases were also extremely abundant during both immediate and extended drought stress and both up and down regulated in A. pinsapo (Cobo-Simón et al., 2022). Protein kinases have been reported to be involved in core-stress signaling pathways related to abiotic stress, such as drought, salinity and extreme temperatures (Zhu, 2016), including ABA production and stomatal closing (Chen et al., 2021). Finally, hub genes related with transcription regulation (lncRNA, Table S4 ) and circadian clock regulation (XAP5 CIRCADIAN TIMEKEEPER, Table S4 ) were found among the up-regulated genes. Several evidence points that circadian clock contributes to plants’ ability to tolerate different types of environmental stress, and to acclimate to them (Grundy et al., 2015), including, again, ABA production and stomatal responses under water stress (Kamrani et al., 2022).

After 20 days of drought stress (E treatment), C. atlantica seedlings showed a higher number DEGs, pointing to a stronger drought-related stress at this stage, compared to immediate drought response, in both analyses (using Fiñana and Fiñana+Dornajo individuals) ( Table 1 ). Again, given the lack of gene modules showing an expression pattern of interest when analyzing both locations together, we will focus the discussion of the gene expression profile after 20 days of drought on Fiñana population. Similar to immediate response, C. atlantica transcriptional profile was consistent with an isohydric drought response ( Figure 9C, D ). Hub genes related to stomatal closure, ABA transportation (NRT1 PTR family-like proteins, Guo et al., 2003 Chiba et al., 2015; acetate butyrate ligases AAE7, Dong et al., 2018) ( Table S6 ), and ABA homeostasis (probable E3 ubiquitin-ligase XERICO, Ko et al., 2006) ( Table S6 ) were activated. In addition, hub genes related to photosynthesis (photosystem I reaction center subunit VI-1, chloroplastic; metal transporter Nramp3-like, Lanquar et al., 2010) ( Table S7 ) and metabolism (mitochondrial adenine nucleotide transporter ADNT1 and phosphoglycerate kinase, both involved in ATP production; lysosomal Pro-X carboxypeptidase, involved in proteolysis; coatomer subunit delta, involved in Golgi function) ( Table S7 ) were inhibited.

In addition, cell wall-related hub genes (probable E3 ubiquitin- ligase HIP1 isoform X1, Sun et al., 2009, exocyst complex component EXO70A1 Synek et al., 2006) ( Table S6 ), together with carbohydrate/sugar metabolism and transportation-related genes (malate synthase, glyoxysomal, acetate butyrate ligase AAE7, peroxisomal, Allen et al., 2011, sugar transporter 7-like, which might contribute to the uptake and recycling of cell wall sugars, Rottmann et al., 2018) ( Table S6 ) were activated in C. atlantica seedlings after 20 days of drought. Other plant growth-related hub genes were, conversely, inhibited (probable purine permease 11, Qi and Xiong, 2013; kinesin KIN-14I, Cooper et al., 2003; beta-tubulin, Eckardt, 2006; probable strigolactone esterase DAD2, involved in plant shoot branching, Hamiaux et al., 2012; cellulose synthase E6, involved in synthesis of cell wall cellulose) ( Table S7 ), including plant polyamine biosynthesis-related genes (ornithine decarboxylase-like, Illingworth and Michael, 2012; Chen et al., 2019) ( Table S7 ).

Sugar transport, distribution and signaling has a critical role in plants under drought stress conditions (Kaur et al., 2021). Abiotic stress and especially drought stress-mediated injury results in reprogramming of sugar distribution across the cellular and subcellular compartments (Kaur et al., 2021). Plant cell walls are composed of carbohydrate polymers, lignin and structural proteins in variable amounts (Tenhaken, 2015). Plants exposed to drought stress reduce shoot growth while maintaining root growth to improve water intake (Miao et al., 2017), a process requiring differential cell wall synthesis and remodeling (Tenhaken, 2015). A cell wall and sugar distribution remodeling in needles, where the RNA from this work was obtained, in order to allow further growth of stressed organs (roots) to improve water uptake under drought conditions, might explain these cell wall and sugar metabolism-related genes activation, and the inhibition of growth and cell division, in aerial tissues. Metabolic pathways implied in shoot:root ratio decrease were observed in A. pinsapo seedling needles, related to their drought resilience as well (Cobo-Simón et al., 2022). We do not have any measurements from this work to test this hypothetical decrease in the shoot:root ratio. Besides, the time spam of this experiment and the limited size of the plants pose limitations to reliably detect it. However, we strongly recommend to test the possible role of a shoot:root ratio decrease in the drought resilience of conifers in future studies, suggested by the results of this and previous works (Cobo-Simón et al., 2022).

Other activated hub genes were related to transcription regulation under environmental stress (regulator of nonsense transcripts 1 homolog, Shi et al., 2012, heavy metal-associated isoprenylated plant 23-like, Braga de Abreu-Neto et al., 2013), including TFs (LHY isoform X1, Moenga et al., 2020) and histone demethylases (FLOWERING LOCUS D isoform X1, related to the epigenetic regulation of flowering in A. thaliana, Luo et al., 2015). Since C. atlantica is a conifer, this histone demethylase might have another function in this species, probably related to the epigenetic regulation of drought response, which would need further investigation. Epigenetic regulation-related genes were also found in A. pinsapo seedlings, associated to their drought resilience (Cobo-Simón et al., 2022). In addition, a gene related to alternative splicing under salt stress response in Arabidopsis was inhibited (serine arginine-rich splicing factor 4-like, Table S7 ) since its overexpression produces hypersensitivity of salt stress (Li Y. et al., 2021)

Other genes related to water stress response were activated (water-stress inducible 1, partial, probable carboxylesterase 15, Cao et al., 2019, linked to drought stress response in Abies alba Behringer et al., 2015), including, again, protein phosphorylation (receptor kinase HAIKU2, associated with drought tolerance in Sequoia sempervirens, de la Torre et al., 2022), similar to the immediate drought treatment.

A gene related to GABA biosynthesis was found among the down-regulated hub genes (glutamate decarboxylase-like, Table S7 ). GABA has a protective role on plants subjected to abiotic stress. It enhances drought stress tolerance by improving photosynthesis, inhibiting reactive oxygen species (ROS) generation, activating antioxidant enzymes, and regulating stomatal opening (Li et al., 2021). However, its potential role in C. atlantica drought response requires further investigation, given its versatile role in plants (Li et al., 2021) and its down-regulation during drought stress in this species ( Table S7 ). This down-regulation might be a consequence of the inhibition of polyamine biosynthesis-related genes, previously described, since GABA can be synthesized through the polyamine metabolic pathway (Li et al., 2021).

Flavonoid metabolism-related hub genes were inhibited after 20 days of drought stress (isoflavone reductase and flavonoid 3,5 -hydroxylase, Table S7 ). It has been suggested that flavonoids protect plants against abiotic stress, and thus, their accumulation is induced rapidly upon drought (Munné-Bosch et al., 2009; Nakabayashi et al., 2014). Nevertheless, flavonoid production has been found to be inhibited in leaves during drought stress in several plants (Salekdeh et al., 2002, Pandey et al., 2010; Rani et al., 2012) including conifers (Fox et al., 2018). The pattern of down-regulation observed in our study may suggest that the plant limits the energy-intensive biosynthesis of secondary metabolites in favor of drought-related survival mechanisms.

Finally, mildew resistance locus o (Mlo) gene (MLO6), which loss-of-function confer broad-spectrum resistance to almost all known isolates of the fungal barley powdery mildew pathogen in plants (Acevedo-Garcia et al., 2014), was inhibited ( Table S7 ). However, it is known that its function is not restricted to plant-powdery mildew interactions (Acevedo-Garcia et al., 2014). Therefore, a further investigation of its possible role in C. atlantica drought response would be advisable.

We hypothesized that differences in gene expression will cause the two contrasting post-drought recovery phenotypes (R and N). The identification of uniquely DEGs in drought sensitive individuals (N) was one of the main findings of this study ( Table 1 ) and supports our initial hypothesis. In addition, we found gene modules showing an expected expression pattern based on our experiment and biologically meaningful functions in both analyses (using only Fiñana population, Figure 8A , 9E, F ; and both populations, Figure 8B , 9G, H ), up ( Tables S9, S11 ) and down-regulated ( Tables S8, S10 ). These results contrast with those obtained in A. pinsapo seedlings using a similar experimental design (Cobo-Simón et al., 2022), where uniquely DEGs were found for R individuals instead. These results suggest different molecular strategies conferring drought resilience in these drought sensitive and biogeographically related conifers (Scaltsoyiannes et al., 1999; Qiao et al., 2007; Linares, 2011; Alba-Sánchez et al., 2018), which would probably explain the wider drought stress tolerances in C. atlantica (Abel-Schaad et al., 2018).

Drought sensitive individuals from both analyses (F and F+D) showed an activation of flavonoid metabolism compared to drought resilient ones (isoflavone reductases, Table S9 ; Transparent testa 12, Debeaujon et al., 2001, Table S11 ; and a Transcription factor TT2, Table S9 , which is involved in flavonoid metabolism when interacting with MYB TFs, Zimmermann et al., 2004). MYB TFs were found among the up-regulated genes in N individuals in both analyses as well ( Tables S9 and S11 ). These results contrast with those obtained for the extended drought response, where the flavonoid metabolism was down-regulated. Hence, they support, again, that limiting the energy-intensive biosynthesis of secondary metabolites in favor of drought-related survival mechanisms might enhance drought resilience in this species. Moreover, an activation of terpene production-related hub genes was found in N individuals from both analyses (7-deoxyloganetin glucosyltransferase-like, involved in irioid biosynthesis, Asada et al., 2013), supporting this hypothesis.

Several hub genes related to ethylene production were up-regulated in N individuals from both analyses (1-aminocyclopropane-1-carboxylate oxidase 5-like, Qin et al., 2007) ( Tables S9, S11 ), including TFs (ethylene-responsive transcription factor 2-like, ethylene-responsive transcription factor 1A-like) ( Tables S9, S11 ). Ethylene is a phytohormone involved in the regulation of gene expression by stress factors and by components of stress signal transduction pathways (Fujimoto et al., 2000). In addition, several hub genes related to antifungal response were activated: chitinases (basic endochitinase, partial; chitinase 2-like, class VII chitinase), peroxidases/oxidoreductases ((R)-mandelonitrile lyase-like; blue copper-like, De Rienzo et al., 2000; cannabidiolic acid synthase-like, Taura et al., 2007), thaumatin-like proteins (Ruiz-Medrano et al., 1992; de Jesús-Pires et al., 2020), beta 1,3 glucanase, Balasubramanian et al., 2012), and others (wall-associated receptor kinase-like 20, Delteil et al., 2016, zinc finger ZAT5-like, zinc finger CCCH domain-containing 17, Gupta et al., 2012; TIR NBS LRR disease resistance, McHale et al., 2006) ( Tables S9, S11 ). These anti-fungal proteins are usually co-expressed (Balasubramanian et al., 2012) and activated in plants by ethylene and other phytohormones under several biotic and abiotic stress conditions (e.g. Kasprzewska, 2003), like wounding, osmotic pressure, cold, heavy metal stress and salt in plants (e.g. Vaghela et al., 2022). Furthermore, ethylene has been found to induce flavonoid accumulation in plants (Lewis et al., 2011). All these findings point to ethylene as an important phytohormone in the drought sensitivity of this species, probably causing a failure in limiting the energy-intensive biosynthesis of secondary metabolites.

Other TFs and transcriptional regulation-related genes were activated in both populations (lncRNA, DRE-binding 2, Table S9 ; two-component response regulator ORR21 isoform X1, probable WRKY transcription factor 4, dynein light chain 1, cytoplasmic, Stelter et al., 2007; nuclear pore complex NUP 35, Table S11 ), including one related to epigenetic regulation (methyltransferase DDB_G0268948, Table S9 ). This methyltransferase supports, again, the role of epigenetic regulation during drought stress in this species, reported in A. pinsapo as well (Cobo-Simón et al., 2022). In addition, an E3 ubiquitin-ligase MIEL1 was up-regulated ( Table S11 ), which mediates degradation of MYB-like TFs, weakening plant defense to pathogens (Marino et al., 2013). Therefore, the expression of this gene could also explain the drought sensitiveness of these individuals.

Furthermore, detoxification proteins were found among the up-regulated genes (detoxification 40-like, Li et al., 2002) ( Tables S9, S11 ) in N individuals, which might be a consequence of their higher levels of drought stress due to their sensitiveness to drought.

Finally, photosynthesis (probable 1-deoxy-D-xylulose-5-phosphate synthase 2, chloroplastic; ORF124), and metabolism-related pathways (probable fructokinase-4, Stein and Granot, 2018, cytochrome P450 82C4, mitochondrial; mitochondrial phosphate carrier 3, mitochondrial; glutamine–fructose-6-phosphate aminotransferase [isomerizing] 2 Vu et al., 2019; hydroquinone glucosyltransferase-like, glycosyl hydrolase; glycerol-3-phosphate dehydrogenase [NAD(+)] GPDHC1, cytosolic, Minic, 2008; acyl carrier 2, mitochondrial-like) were found among the up-regulated hub genes in N individuals when analyzing Fiñana and Dornajo populations ( Table S11 ). Moreover, stomatal closure pathways were activated (sphingosine kinase, Coursol et al., 2005; probable phosphatase 2C 76 isoform X2, Yang et al., 2018; Hsu et al., 2021) ( Table S11 ). It is worth noting that a phosphatase 2C was related to drought resilience in A. pinsapo and was also up-regulated under drought stress conditions in Pinus halepensis (Fox et al., 2018). These findings suggest that drought sensitivity in C. atlantica might be caused by failing in down-regulating photosynthesis and metabolism under drought stress conditions, which will negatively affect the whole plant carbon status.

One of the main objectives of this work was the identification of candidate genes mechanistically involved in the drought resilience of C. atlantica. However, our findings provide a set of drought-sensitiveness candidate genes instead, since we found uniquely up and down regulated genes in non-resilient individuals. These drought-sensitiveness candidate genes are key for the design of conservation and management strategies, as well as breeding programs. To ensure their success, it is not only necessary to know the genes conferring tolerance, but those that we aim to avoid, that is, those conferring sensitivity.

Given that both post-drought responses were subjected to the same drought treatment, we hypothesized SNPs as polymorphisms associated to their differentially expressed genes, and thus, related to their drought resilience. In addition, we hypothesized a drought-related rapid local adaptation in both C. atlantica locations given the lack of differences in their fitness despite the contrasting precipitation rates. Both hypotheses were supported by our results, showing genomic differences between responses and locations ( Figure 11 ).

The lack of clear gene expression differences between treatments and responses, with the exception of N individuals ( Figures 4 , 5 ), in contrast with A. pinsapo results which showed a clear separation of all treatment/responses (Cobo-Simón et al., 2022), led us to hypothesize a more constitutive gene expression related to drought in C. atlantica, compared to A. pinsapo. This drought-related constitutive gene expression might lead to its wider drought tolerances (Abel-Schaad et al., 2018). The lack of genomic differences between treatments and responses, with the exception, again, of the molecular singularity of N individuals, based on the SNP results ( Figure 11 ) supports this hypothesized constitutive gene expression conferring drought resilience on this species. In addition, the genomic differences between locations ( Figure 11 ), consistent with their gene expression patterns ( Figure 4 , 5 ), support our hypothesis of a rapid local adaptation to the contrasting precipitation rates. This hypothetical rapid local adaptation might have been favored by this hypothesized drought-related constitutive genetic diversity. Evidence of local adaptation to drought has been found in other relict conifers (A. pinsapo, Cobo-Simón et al., 2021). Despite the low number of samples used in the SNP calling analysis, which prevents us to draw statistically robust conclusions based on SNPs, these findings point to enough adaptive potential of this species under a climate change scenario, although further investigation is needed to fully validate it.

As expected, few SNP genes were found among the DEGs ( Tables 3 , 4 ), since most of the causative variation of this DEGs is not located in the transcriptome, but in regulatory regions of the genome, both in cis and in trans. These results are similar to those obtained for A. pinsapo (Cobo-Simón et al., 2022). The annotation and GO term enrichment analysis of the SNP-related genes between responses ( Table S12 ) and locations ( Table S13 ) showed similar functions to the DEGs ( Figure 10 , Tables S4-S11 ). Nevertheless, the abundance of transposon-related genes in both comparisons, which were absent on the DEGs, was remarkable. These results point to a role of TEs in the local adaptation and the constitutive molecular response conferring drought resilience of C. atlantica.

TEs have been found to have a key role on genome and adaptive evolution of conifers, as an important force in shaping gene regulatory networks, their downstream metabolic and physiological outputs, and thus their phenotypes (Liu and El-Kassaby, 2019). In fact, the large genome size of conifers seems to result from the slow and steady accumulation of a diverse set of long-terminal repeat TEs, possibly due to a lack of an efficient elimination mechanism (Nystedt et al., 2013). TEs are DNA fragments that can move and amplify their copy number within a host genome (Casacuberta and González, 2013). As a result, they are a major and powerful source of genomic mutations, causing the evolution of new genes and their functionalities and generating novel regulatory sequences, such as promoters and enhancers (Liu and El-Kassaby, 2019). These in cis regulatory modifications change the regulatory and/or epigenetic environment, leading to modified gene expression (Mita and Boeke, 2016). Furthermore, TEs could generate genetic variation in response to environmental and genomic perturbations, which could facilitate adaptation (Casacuberta and González, 2013). Hence, it has been suggested that TEs may lead to adaptation to novel habitats in the face of limited genetic variation, such as invasive species (the so-called genetic paradox of invasive species, Stapley et al., 2015), which can also explain the local adaptation of relicts, such as C. atlantica. Hence, we strongly recommend to further investigate the potential role of TEs in rapid local adaptation and drought resilience of C. atlantica in particular, and conifers in general.

Here we offer the first genomic resources for the climate change threatened relict C. atlantica, including its first de novo transcriptome assembly, drought-related transcriptomic profile and candidate genes of drought sensitiveness. Our results contrast with those obtained in a similar experimental design in the relict and biogeographically related fir A. pinsapo (Cobo-Simón et al., 2022), where uniquely expressed genes were found for drought resilient individuals instead. This finding is remarkable, specially from an evolutionary perspective, since it suggests two different molecular strategies to the same abiotic stress in both conifers (drought), which would probably explain their different drought tolerances. In addition, our results suggest a constitutive molecular drought response conferring drought resilience in C. atlantica and a rapid local adaptation to drought in the two studied populations. TEs might play a key role in this rapid local adaptation and its standing genetic variation related to drought resilience, hypothesized for relict species (Kawecki, 2008). Thus, we recommend further research using a statistically robust number of biological replicates per treatment/response in Dornajo population to elucidate potential drought-related molecular and functional differences between populations. Moreover, we recommend further validation of this probable local adaptation suggested by our results, and the potential role of TEs, by using a wider sample size in both locations, by integrating multi-omic data (e.g. epigenomic, metabolomic) with other traits/environmental variables of interest (e.g. GWAS, GEA) and by studying the C. atlantica natural stands from northern Morocco.

Our findings, especially when complemented with those from A. pinsapo (Cobo-Simón et al., 2022), shed light on the adaptive potential of trees in general, and relict conifers in particular, to the projected increasing drought caused by current climate change, which is particularly concerning in Mediterranean and other semi-arid climates. This information will allow us better predict their future perspectives facing current global warming.

Finally, our results offer a diversity of drought-related genomic resources, including candidate genes of drought sensitiveness, which will help guide the design of conservation strategies and breeding programs for this climate-change threatened and drought-sensitive conifer. Furthermore, they can be used to advance future drought-resilience research, such as comparative genomic studies and functional transfer across species, which will also shed light on trees’ molecular adaptive potential to current climate change.