1. Introduction

Forest ecosystems around the world can be affected by extreme climatic events, such as heatwaves and severe droughts (Lindner et al., 2010), especially in boreal (Gauthier et al., 2015) and temperate forests (Ning et al., 2022). Trees respond to such occurrences in a multidimensional way, including various mechanisms, to overcome the negative effects and restore the forest to their optimal state. These abilities are referred to collectively as resilience (Fuller and Quine, 2016). The processes involved in resilience include optimization of water transport efficiency in the xylem (Jupa and Pokorná, 2024), implementation of osmotic adjustment strategies (Shaffique et al., 2024), initiation of stomatal regulation processes (Pernicová et al., 2024) and stimulation of antioxidant metabolism (Cui et al., 2024). In situations of water scarcity, trees may reduce or even temporarily halt their growth to conserve resources. This is reflected in the ring width size, in which reduced growth or even no growth can be observed in extremely dry years (Fritts, 1976).

Examination of tree rings, especially in conifers, is essential for analysis of drought resilience (Depardieu et al., 2024). During interannual growth episodes, trees undergo cell differentiation in earlywood and latewood (Babushkina et al., 2024), with unique characteristics in cell size and wall thickness driving the formation of ring width features. These can be measured to determine the dendrophenotypes (Heer et al., 2018).

Dendrogenomics is a new interdisciplinary field of research that integrates dendrochronology, dendroecology, dendroclimatology, genetics and genomics (Krutovsky, 2022). Analysis of dendrophenotypes, combined with genetic studies, has become an essential tool for identifying genotypes resilient to dry conditions (Depardieu et al., 2020). By examining the marks left by past droughts in the annual rings, it is possible to determine the trees' responses to extreme climatic events and also to associate these responses with the expression of genes related to drought resilience (Krutovsky, 2022). Such associations may reveal specific genetic variants that influence the ability of trees to survive and adapt to water-scarce environments (Candido-Ribeiro and Aitken, 2024). Combining dendrochronological and genomic data will not only yield a better understanding of how trees adapt to changing environments, but will also open up new avenues for the selection of more suitable genotypes to face challenges related to climate change, such as increased severity of droughts (Krutovsky, 2022).

The two fragmented, isolated, rare, endemic and clearly genetically distant (Lockwood et al., 2013) Picea martinezii Patterson and P. mexicana Martínez are recognized as relict species from the last glacial period. Both species are also strongly susceptible to the impacts of climate change and classified as endangered according to the IUCN Red List (2001, 2004) and NOM-059-SEMARNAT-2010 (SEMARNAT, 2019). In total, about 89, 266 individuals, including seedlings (of height 30 cm), saplings and mature trees of P. martinezii have been distributed over a total area of only 157 ha (ha) in four isolated populations at medium–high elevations of between 1820 and 2515 m in the Sierra Madre Oriental. By contrast, the total number of P. mexicana trees was estimated to be about 39, 059, covering 173 ha in three isolated populations at the highest elevations of 3096 and 3528 m in the Sierra Madre Occidental and Sierra Madre Oriental, respectively (Mendoza-Maya et al., 2022; Wehenkel et al., 2022, 2023). Like many other rare and endemic plant species, these conifer species are subject to conservation measures. They may therefore be useful as model species to examine the resilience of rare and endemic plants to severe ecosystem disturbances including climate change (Mendoza-Maya et al., 2022).

The study aimed to determine the dendrogenomic resilience mechanisms in Picea martinezii and P. mexicana for the first time through by analysis of annual rings and their association with single nucleotide polymorphism (SNP) markers derived from genotyping by sequencing (GBS). We hypothesized that combined analysis of dendrochronological and genomic data will enable identification of the most important genes involved in resilience processes (Krutovsky, 2022).

2. Materials and methods 2.1. Tree sampling

For annual-growth ring analyses, 72 increment cores were collected at a height of 130 cm from the base of the stem. These samples were obtained from 10 Picea martinezii trees, aged between 109 and 201 years, in the populations of La Encantada (four samples), El Butano (three), and Agua Fría (three), as well as from 14 P. mexicana trees, aged 74–232 years, in the populations of La Marta (10 samples) and El Coahuilón (four) (Table 1). The number of increment core samples was restricted due to the protected status and the limited number of mature trees of both species. For GBS analysis, needle samples were collected from 114 trees located in the four recorded P. martinezii populations: La Encantada (28 samples), El Butano (30), Agua Fría (27), and Agua de Alardín (29). Additionally, samples were collected from 87 individuals in three documented P. mexicana populations: La Marta (28 samples), El Coahuilón (29), and El Mohinora (30) (see details in Mendoza-Maya et al., 2024). This needle collection also included the trees from which increment core data were obtained.

Once mounted, the year of formation of each ring was assigned using the cross-dating procedure (Stokes and Smiley, 1968). Using the mark obtained, the ring widths (earlywood, latewood and total ring width) were measured to an accuracy of 0.001 mm by using the VELMEX sliding stage measurement system (Velmex, Inc., Bloomfield, NY, USA). Measurements were made by microscopic examination of the samples (in a Fotgear X4 1600× digital microscope and a Zeiss Stemi 2000-c binocular microscope) with a 0.1 mm calibration micrometer scale (Nikinmaa et al., 2020).

The tracheids were quantified in a digital microscope (Fotgear X4 1600×) at magnifications ranging from 40 to 60×. For each annual ring, tracheids were counted manually across the ring width from earlywood to latewood. This procedure was replicated three times per annual ring within the selected core area. Mean values were calculated from the values of the three repetitions.

The resilience analysis, based on the modified formula by Lloret et al. (2011), distinguished between resistance, recovery, and resilience, which is challenging with the data resolution provided by annual tree rings. Specifically, resistance (the immediate reduction in growth during the drought) and recovery (the speed of growth increase afterward) are highly sensitive to year-to-year fluctuations and anomalies. This can lead to numerical instability and limits the interpretability, especially if the growth stops completely (e.g., in the case of missing or extremely narrow rings). By contrast, resilience offers a more robust and integrative measure. Therefore, this was one of the reasons why we only examined resilience.

Peltier et al. (2016) reported changes in drought sensitivities after droughts persisting in some tree species for five years. Instead of three years (Lloret et al., 2011), therefore, the selected core area comprised the four annual rings before each drought event and the four rings produced after the drought event (Serra-Maluquer et al., 2018; DeSoto et al., 2020). Nevertheless, one-year measurements were also calculated as they can sometimes be relevant (DeSoto et al., 2020). Four resilience indices were calculated from the ring width and the number of linear tracheids per ring, providing a comprehensive view of the impact of stress on the growth and internal structure of trees.

The following resilience indices were used in the study (Table 1):

The ratio of the sum of the four annual ring widths after (W_a) and the sum of the four annual ring widths before (W_b) the severe drought, expressed as RWM = W_a/W_b.

The ratio of the sum of the four annual tracheid widths after (WT_a) and the sum of the four annual tracheid widths before (WT_b) the severe drought, expressed as TWM = WT_a/WT_b.

The ratio of the first annual ring width after (R_a) and the last annual ring width before (R_b) the severe drought, expressed as RAB = R_a/R_b.

The ratio of the first annual tracheid width after (T_a) and the last annual tracheid width before (T_b) the severe drought, expressed as TAB = T_a/T_b.

2.2. DNA extraction, sequencing, genotyping and filtering of SNPs

The cetyltrimethylammonium bromide (CTAB) protocol described by Doyle and Doyle (1987) was employed to isolate DNA from 15.0 mg of needle tissue obtained from collected samples. The modifications introduced aimed primarily to utilize readily available in-stock reagents serving equivalent functions (e.g., polyvinylpyrrolidone instead of polyethylene glycol, or LiCl in place of NaCl), as well as to enhance DNA isolation and purity through the addition of 2-mercaptoethanol and Proteinase K. Although a direct comparison between the original method described by Doyle and Doyle (1987) and our modified protocol was not conducted, we considered these methodological adjustments in a routine protocol noteworthy, given that the DNA quality and concentrations required for sequencing were successfully achieved. These modifications allowed the determination of DNA concentrations in volumes suitable for mini-preparations (i.e., small-scale DNA extractions performed in 1.5 ml microcentrifuge tubes). Specifically, the following changes, as implemented by Mendoza-Maya et al. (2024), were applied: polyethylene glycol (PEG) was replaced with 1.0% polyvinylpyrrolidone (PVP40); LiCl was substituted for 5 M NaCl; and 5.0 μl of 2-mercaptoethanol and 50 μl of Proteinase K were added.

DNA was processed following the Genotyping by Sequencing (GBS) reduced representation method, as described by Abed et al. (2019). Briefly, 200 ng of DNA per individual was double-digested with PstI (R0106L) and MspI (R3140L; both from New England Biolabs). Barcodes were ligated with T4 DNA ligase (New England Biolabs) in the same plate. Samples were pooled, purified, and amplified at the "Plateforme d'analyses génomiques" at the "Institut de biologie intégrative et des systems-IBIS (Université Laval)". Paired-end sequencing (2 × 150 bp) was finally performed on an Illumina NovaSeq 6000 S4 (1 lane) System at the "Centre d'expertise et de services, Génome Québec" by implementing 330 cycles of sequencing including paired-end reads and index reads.

After processing raw sequencing data (initial quality check, filtering low quality sequences and trimming possible presence of adapters used for sequencing), reads were de novo assembled using Stacks and the references generated for Picea martinezii and P. mexicana were subsequently used to perform the variant calling (Catchen et al., 2013) for each species. The resulting VCF files were filtered by means of VCFTOOLS v0.1.17 (Danecek et al., 2011), using the following parameters: –mac 2, –min-meanDP 10, –max-meanDP 120. Final VCF files included 19, 057 and 31, 678 SNPs for P. martinezii and P. mexicana, respectively. For more details, see Mendoza-Maya et al. (2024).

2.3. Identification of outlier SNPs

Four approaches were used to identify outlier SNPs in Picea martinezii and P. mexicana: i) the package PCAdapt v.4.3.3 (Luu et al., 2016; Privé et al., 2020) for R (R Core Team, 2022); ii) the genotype environment association (GEA) method by using latent factor mixed models (LFMM2), included in the R package LEA (Caye et al., 2019; Gain and François, 2021); iii) the Local Outlier Factor (LOF) algorithm of the Python package scikit-allel (http://scikit-allel.readthedocs.org) (Breunig et al., 2000; Buitinck et al., 2013); and iv) the R package OutFLANK v.0.2 (Lotterhos and Whitlock, 2015).

PCAdapt provides advantages over other genomic scanning software, such as BayeScan, hapflk, OutFLANK and sNMF, due to its speed, ability to handle missing data and individualized approach that eliminates the need to group individuals into populations (Luu et al., 2016). The package includes a method for detecting outlier SNPs based on principal component analysis (PCA), accounting for population structure (Luu et al., 2016). PCAdapt is based on Mahalanobis distance (D) statistics, from which a z-score vector is derived by regressing each SNP against K (K = 2) principal components (Li et al., 2017). The number K of principal components was chosen by applying Cattell's rule (1966). P-values were obtained by transforming D based on a chi-square distribution (Cattell, 1966). Bonferroni correction was applied to adjust the obtained p-values, considering SNPs with Bonferroni-corrected p-values < 0.05 as candidates under selection.

LFMM2 was used to evaluate GEAs. Prior to analysis, non-collinear bioclimatic and edaphic variables were selected using the absolute Spearman rank correlation coefficient criterion (|r_s|) < 0.7, resulting in six bioclimatic variables and thirteen edaphic variables. The LFMM2 analysis was executed with four latent factors selected on the basis of the variance explained by each principal component for climatic and edaphic variables separately. The resulting p-values were adjusted using Bonferroni correction to minimize false positives, considering significant SNPs with an adjusted p-value less than 0.05.

The LOF algorithm included 20 nearest neighbors, with the choice of k = 20 aiming to balance the detection of significant deviations without sacrificing precision. The kNN-LOF algorithm employs weightings based on sequential distances to adjust the influence of neighbors, enabling k = 20 to perform effectively even in datasets with irregular distributions (Xu et al., 2022), a contamination rate of 1%, a Z threshold of 2, and Manhattan distance as the measurement metric (Liu et al., 2022).

The OutFLANK detects outlier SNPs using F statistics, with standard parameters: a mean F statistic value for neutral loci of 0.05, a proportion of neutral loci of 0.95, and minimum degrees of freedom (df_min) of 10. Considering that OutFLANK automatically removes loci with extreme Fst values (both lower and upper) to infer the neutral distribution, by default it trims 5% of loci at each tail of the Fst distribution, assuming that the central 90% corresponds to neutral loci. However, the mentioned 95% proportion of neutral loci arises from an implicit contamination rate of 5% (1–0.95), which represents the maximum fraction of non-neutral loci the algorithm can tolerate without biasing the estimation of the neutral distribution (Whitlock and Lotterhos, 2015).

2.4. Association of resilience indexes with SNP genotypes

The association analyses of the resilience indices with the SNP genotypes were only carried out for the trees for which the annual ring and GBS analyses were performed. The genotypic values of the outlier SNPs were numerically coded into classes "0", "1" and "2", for homozygotes for the major allele (AA), heterozygotes (Aa), and homozygotes for the minor or alternate allele (aa), respectively. The missing data were imputed with the median value of the dataset of the respective SNP. Nonparametric Spearman's correlation coefficient (r_s) and their p-values were calculated between each of these SNPs and the four resilience indices (RWM, TWM, RAB and TAB) in Picea martinezii and P. mexicana using the Hmisc package in R v.4.3.1. A positive r_s indicated a positive association between outlier SNP genotypic values and resilience indices, while a negative r_s suggested an inverse relationship (Whitacre, 2012; Smeeth et al., 2021; Flotildes et al., 2023; Thorogood et al., 2023). PCA was performed including outlier SNPs with significant r_s with resilience indices at α = 0.05 to calculate the first Principal Component (PC1) and to determine r_s and their p-value between PC1 and resilience indices.

To test the potential multivariate associations of resilience indices with outlier SNPs with significant r_s with resilience indices at α = 0.05, moreover, pseudo coefficient of determination (R_rf²) and its 99.999% confidence interval (CI_99.999%) were determined with Random Forest for regression (RF) (Breiman, 2001) by the randomForest package (Liaw and Wiener, 2002) in R.

2.5. Identification and analysis of SNP-containing sequences

The genome sequences of Picea martinezii and P. mexicana and the resilience-associated SNPs were obtained from the de novo assembled references. These sequences, ranging from 120 to 200 nucleotides, were compared with the genomic reference of Picea abies (Nystedt et al., 2013) detected with the Basic Local Alignment Search Tool (BLAST) (Altschul et al., 1990) implemented in http://PlantGenIE.org (Sundell et al., 2015). The resulting sequences were aligned using AliView (v.1.28) (Larsson, 2014) to detect potential deletions.

3. Results 3.1. Association between SNP and resilience

The mean values of the four resilience indices studied in Picea martinezii and P. mexicana, respectively, were as follows: RWM = 0.98 and 0.95, TWM = 1.31 and 1.05, RAB = 0.85 and 1.00 as well as TAB = 1.22 and 0.88. In P. martinezii, three indices were significantly positively correlated: RBA vs. TWM, RAB vs. RWM and TAB vs. RWM, but in P. mexicana, on the other hand, there was only a significant correlation between the TAB and RWM indices (p < 0.05) (Table 1) (see Fig. 1).

Using PCAdapt, 162 outlier SNPs were identified in Picea martinezii and 87 in P. mexicana. LFMM 2.0 identified 1678 outliers in P. martinezii and 8272 in P. mexicana, while Scikit-allel identified 483 outliers in P. martinezii and two in P. mexicana. In contrast, outflank did not detect any outliers in either species (Fig. 2).

In Picea martinezii, PCAdapt identified five outlier SNPs correlated with some of the resilience indices (p < 0.05). One of these outlier SNPs was also detected by LFMM 2.0. The PC1 of the principal component analysis (PCA) including the five outlier SNPs was statistically significantly associated with the TWM resilience index (p = 0.0007), after Bonferroni correction. The CI_99.999% of [0.16–0.79] of the pseudo coefficient of determination (R_rf²) of the multivariate association of TWM with these five outlier SNPs also shows that a robust relationship exists.

In Picea mexicana, PCAdapt identified three outlier SNPs correlated with some of the resilience indexes (p < 0.05). However, the relationship between PC1 of these three outlier SNPs was not significantly correlated with any of the four resilience indices after Bonferroni correction (minimum p = 0.015). The CI_99.999% of R_rf² of largest multivariate association (i.e., of TAB with these three outlier SNPs) was [0.00001, 0.76] suggesting that the association is not strong enough to be considered statistically significant. None of the individual outlier SNPs were significantly correlated with the resilience indices after Bonferroni correction (Table 2).

3.2. Identification of candidate genes for resilience in Picea martinezii and P. mexicana

In the absence of reference genomes for Picea martinezii and P. mexicana, the annotated P. abies genome was used as reference to perform BLAST search and provide a functional annotation of the sequences where the SNP outliers were located. With genetic similarities of almost 100%, these resiliencies associated SNPs (p < 0.05) for P. martinezii were located in genomic sequences of Picea abies genes that encode the following proteins: reticulon-like protein B22 (RTNLB22) (MA_10435443g0020), pollen-specific leucine-rich extensin (MA_10428704g0010), proline transporter 2 (MA_95628g0010), LisH/CRA/RING-U-box domain-containing protein (MA_10002g0020), and ornithine decarboxylase (MA_10427924g0010). In P. mexicana, genome sequences and the resilience associated SNPs (p < 0.05) were very similar to P. abies cyclin-dependent kinase B1-1 (MA_63670g0010), an uncharacterized GPI-anchored protein (MA_19048g0010), and an abscisic acid receptor PYL8-like (MA_9600951g0010) (Table 2).

4. Discussion

The study findings support the proposed hypothesis that a combined analysis of dendrochronological and genomic data enables the identification of key genes involved in resilience processes (Krutovsky, 2022), owing to observation of a significant association in Picea martinezii between the first principal component (derived from five outlier SNPs) and the TWM drought resilience index for the four drought events in 1962, 1989, 1998, and 2011 (Table 2). This suggests a genetic basis for TWM. As only one of these five SNPs was identified by PCAdapt and LFMM2, but not by the other outlier detection algorithms, the evidence for strong adaptive selection remains limited. Moreover, since GBS represents a reduced version of the entire genome, it could only explain a small number of putative adaptive outlier SNPs (Visscher et al., 2017). Nevertheless, these results are best explained by selection in heterogeneous environments, where different alleles are favored in different trees due to local adaptation to varying spatial and temporal conditions (Alberto et al., 2013; Collevatti et al., 2019; Vieira et al., 2022). This contrasts with classical balancing selection, which generally maintains genetic diversity within a single, homogeneous population (Delph and Kelly, 2014).

Recent studies have also used dendrogenomics to find relationships in other tree species, such as Fasanella et al. (2021), reporting individual-based dendrogenomic correlations of forest dieback in Nothofagus dombeyi (Mirb.) Blume caused by extreme droughts. Using the example of relict forests of Abies marocana Trab., Méndez-Cea et al. (2023) showed the potential of dendrogenomic analysis to improve our understanding of the adaptive capacity of drought-sensitive forests. Additionally, genetic adaptation mechanisms have been demonstrated in response to climatic stress in Larix sibirica Ledeb (Novikova et al., 2023), and in Pinus sibirica Du Tour in response to the climatic and pest outbreak stresses (Novikova et al., 2024).

In addition, the following possible causes of the presence of outlier SNPs associated with TWM without selection pressure cannot be ruled out: i) linkage disequilibrium, where the observed SNPs could be linked to genes that influence TWM (Wills, 2007); ii) pleiotropic effects causing some SNPs to influence multiple traits, considered a constraint in the evolution of adaptive phenotypes (Auge et al., 2019); iii) epigenetic mechanisms influencing gene expression related to TWM (Bird, 2007); iv) genetic–drought interactions affecting the expression of some genes that regulate tracheid development (Gillespie, 2004); and v) recurrent mutation, where certain loci or regions of the genome have a higher probability or tendency to undergo changes in their DNA sequence. This may be due to specific characteristics of the DNA in those regions, such as repetitive sequences or secondary structures and these mutations have occurred repeatedly over time and have been maintained (Gillespie, 2004). Epigenetic mechanisms are non-permanent modifications of the DNA sequence and may modulate gene expression related to TWM in response to environmental factors, promoting phenotypic plasticity without changing the underlying genetic variation (Bird, 2007).

The 160–200 bp DNA sequences of the five outlier SNPs in P. martinezii were homologous with P. abies genes encoding: i) reticulon-like protein B22 (RTNLB22), which is involved in endoplasmic reticulum configuration (Howell, 2013; Azmat et al., 2024); ii) pollen-specific leucine-rich extension, which is involved in cellular signalling and abiotic stress responses in plants (Hui et al., 2024); iii) the LisH/CRA/RING-U-box domain, which assists in protein ubiquitination (Hatakeyama et al., 2001); iv) isoform X1 of the proline transporter, which connects metabolic pathways through ornithine and glutamate and promotes growth (Verslues and Sharma, 2010; Chen et al., 2024); and v) ornithine decarboxylase (ODC), a key enzyme in polyamine and nicotine biosynthesis (Bajguz and Piotrowska-Niczyporuk, 2023) (Table 2). All these proteins have been associated to drought stress tolerance metabolism (Singh et al., 2022), which is a crucial part of a resilience response in plant species (Movahedi et al., 2023).

The RTNLB22 encoding gene was the only one linked to a sequence with a SNP found by two out of the five outlier detection algorithms used in the present study (Table 2). Reticulon (RTN) is a ubiquitous protein found in all eukaryotes and is involved in cellular processes such as apoptosis and intracellular trafficking (Nziengui and Schoefs, 2009). Under stress conditions, RTN responses, including the unfolded protein response (UPR) and endoplasmic reticulum-associated degradation (ERAD), are activated to help plants adapt to long-term environmental changes (Thibault and Brandizzi, 2024; Azmat et al., 2024). The reticulon-like protein B22 may play a role in homeostasis in response to climatic variations and has been observed to induce high curvature of endoplasmic reticulum tubules, thus facilitating cellular transport (Hu et al., 2008; Shao and Shun, 2024). Additionally, RTNs mediate stress responses in various species, including Hibiscus sabdariffa (Yong and Atheeqah-Hamzah, 2024) and Pseudotsuga menziesii (Compton et al., 2023), and are present in agricultural crops such as maize, in which they contribute to drought resistance (Tian and Qin, 2023).

In Picea mexicana, the association between three outlier SNPs and their PC1 and resilience indices was not significant after Bonferroni correction (Table 2). However, homology analysis indicated that these three SNPs were located in genes coding for CDKB 1-1, GPI-anchored protein and abscisic acid receptor PYL8. CDKB 1-1 may influence resilience by affecting cell cycle regulation and recovery from drought-induced growth cessation (Valles et al., 2024; Romeiro Motta et al., 2024; Farooq et al., 2024). The GPI-anchored is a class of proteins involved in stress responses (Takahashi et al., 2016; Li et al., 2017; Chang et al., 2021). The abscisic acid receptor PYL8 is crucial for drought stress signaling, and variations in PYL8 could modulate stomatal closure and other mechanisms involved in drought tolerance (Finkelstein et al., 2002; Robert-Seilaniantz et al., 2007; Ton et al., 2009; Fujii et al., 2009; Fujita et al., 2009; Nakashima et al., 2009). Consequently, although these associations were not significant according to Bonferroni's strict criteria, the location of outlier SNPs detected in genes involved in resilience metabolism suggests a possible biological relevance (Table 2).

The high similarity between the genomimc sequences containing resilience-related outlier SNPs from the two Mexican spruces studied here and Picea abies, despite the 20–25 million years diversification between these species (late Oligocene/early Miocene), can be explained by a combination of evolutionarily conserved genes, parallel evolution and the relatedness of the species. Although the species are geographically distant from each other, similar ecological selection conditions have probably shaped their resilience and stress response genes, leading to their conservation over million years (Lockwood et al., 2013).

Although our ratios of the first annual ring or tracheid width after (Ra or Ta) and the last annual ring or tracheid width before (Rb or Tb) the severe drought (RAB and TAB) should be more influenced by random factors and less reliable (Lloret et al., 2011) than ratios based on multi-year measurements, we found similarly high correlations with some outlier SNPs as with the four-year indices. These SNPs were different from the SNPs that were correlated with four-year indices and were associated with important resilience-related proteins (Table 2). Therefore, RAB and TAB should not be excluded in future resilience studies, as also stated by DeSoto et al. (2020) in a study on the growth resilience of trees under drought.

5. Conclusions

This study reveals some genetic mechanisms associated with drought resilience in the threatened species Picea martinezii and P. mexicana by integrating dendrochronology and genomics. The conservation status and limited number of old trees of both species limited the number of the increment core samples (Table 1) and, thus, the probability of detecting true SNP-resilience associations. Despite this sample size limitation, this study shows that genetic markers for resilience can be detected. The study also illustrates the potential of genome analysis to provide important insights even under difficult conditions and lays the basis for further studies with larger and more diverse samples.

The discovery of genetic markers associated with drought resilience highlights the importance of preserving genetic diversity in the isolated populations of both endangered species. Given the high risks posed by climate change and human activity, targeted conservation measures/strategies must be developed to protect genetic diversity and the species' ability to adapt to future climate conditions. Identifying resilience genes could contribute to breeding programs or reforestation strategies that help these species withstand.

Acknowledgments

This study was conducted thanks to the funding from the Mixed Fund of the National Council of Humanities, Sciences, and Technologies of Mexico and the National Forestry Commission (CONACYT-CONAFOR-2017-4-292615), awarded to Christian Wehenkel. Additionally, SECIHTI provided a graduate scholarship to Carlos Alberto Segura Sánchez (776540).

The present work is the result of the joint efforts of the Working Group on Forest Genetic Resources of the North American Forestry Commission. We express our sincere gratitude to the landowners and ejidatarios of the communities and private properties where the fir populations are located. Their invaluable hospitality, support during fieldwork, and access to the populations were fundamental for the development of this study. We are thankful to Oscar Alfredo Diaz-Carrillo and Anton Christian Wehenkel-Lara for their assistance in data collection.

CRediT authorship contribution statement

Carlos Alberto Segura-Sanchez: Writing – review & editing, Writing – original draft, Visualization, Resources, Formal analysis, Data curation, Conceptualization. Javier Hernández-Velasco: Writing – review & editing. José Villanueva-Díaz: Writing – review & editing. Víctor Chano: Writing – review & editing. José Ciro Hernández-Díaz: Writing – review & editing. Eduardo Mendoza-Maya: Writing – review & editing. Artemio Carrillo-Parra: Writing – review & editing. Christian Wehenkel: Writing – review & editing, Writing – original draft, Supervision, Resources, Project administration, Formal analysis, Data curation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.