Signatures of natural selection in a foundation tree along Mediterranean climatic gradients.

Temperature and precipitation regimes are rapidly changing, resulting in forest dieback and extinction events, particularly in Mediterranean-type climates (MTC). Forest management that enhance forests' resilience is urgently required, however adaptation to climates in heterogeneous landscapes with multiple selection pressures is complex. For widespread trees in MTC we hypothesized that: patterns of local adaptation are associated with climate; precipitation is a stronger factor of adaptation than temperature; functionally related genes show similar signatures of adaptation; and adaptive variants are independently sorting across the landscape. We sampled 28 populations across the geographic distribution of Eucalyptus marginata (jarrah), in South-west Western Australia, and obtained 13,534 independent single nucleotide polymorphic (SNP) markers across the genome. Three genotype-association analyses that employ different ways of correcting population structure were used to identify putatively adapted SNPs associated with independent climate variables. While overall levels of population differentiation were low (FST  = 0.04), environmental association analyses found a total of 2336 unique SNPs associated with temperature and precipitation variables, with 1440 SNPs annotated to genic regions. Considerable allelic turnover was identified for SNPs associated with temperature seasonality and mean precipitation of the warmest quarter, suggesting that both temperature and precipitation are important factors in adaptation. SNPs with similar gene functions had analogous allelic turnover along climate gradients, while SNPs among temperature and precipitation variables had uncorrelated patterns of adaptation. These contrasting patterns provide evidence that there may be standing genomic variation adapted to current climate gradients, providing the basis for adaptive management strategies to bolster forest resilience in the future.

INTRODUCTION

Climate change is a key pressure on ecosystem persistence and function (Brondizio et al., 2019; Urban, 2015). The shift in climate trends will have an impact on ecosystem structure, potentially making organisms more susceptible to the effects of extreme climate events (Harris et al., 2018; Pacifici et al., 2015). Precipitation patterns are changing in heterogenous ways, with some areas becoming wetter and others drier. While global surface temperature is predicted to rise by 1–4°C on average by the end of the current century, the level of temperature rise is also heterogeneous depending on various factors (e.g., latitude, elevation). In addition, the frequency of extreme events such as heatwaves and droughts have increased over recent decades in several regions of the world (IPCC, 2021). Because these changes are spatially assorted, understanding broad patterns of adaptation across landscapes is often challenging.

Mediterranean‐type climates (MTC) are defined by reliable precipitation and temperature regimes, with predictable summer periods of low rainfall and hot temperatures, and winter periods of high rainfall and moderate temperatures. Ecological studies in the Mediterranean basin consistently identify more frequent drought periods, together with warmer temperatures, as main drivers for declines in oaks (Quercus spp.) (Corcobado et al., 2014; Gentilesca et al., 2017) and pines (Pinus spp.) (Camarero et al., 2018). In the South‐west Western Australia biodiversity hotspot, a region with MTC, the 2010–2011 extreme drought and heatwave conditions resulted in large‐scale forest collapses in eucalypts (Matusick et al., 2013). While some variation in climatic factors exists in natural systems (Staudinger et al., 2013), the rapid and extreme shifts associated with anthropogenic climate change are challenging for most organisms to persist (Carlo et al., 2018; Levin & Poe, 2017).

If new climatic scenarios are no longer suitable for species to maintain their normal ecology and physiology, they either shift their geographical range or, in worst case scenarios, go extinct (Aitken et al., 2008; Bellard et al., 2012). Although, species may persist through enhanced physiological tolerance, phenotypic plasticity and/or genetic adaptation (Anderson et al., 2011; Christmas et al., 2016). Maintenance of standing genetic variation (within population allelic variation at a locus) is a key factor for adaptation to changing conditions in native habitats (Chhatre et al., 2019; Guzella et al., 2018) and for persistence through environmental stressors over generations (Kremer et al., 2012; Sexton et al., 2011). Genetic variation is critical for ecological adaptive capacity – the potential and ability to adjust to, and persist through, external factors – and consequently, the evolutionary potential of the species (Reed et al., 2011). Evolution to a specific environment through natural selection results in patterns of local adaptation, when a local population experiences higher fitness compared to nonlocal counterparts (Kawecki & Ebert, 2004).

Recent improvements in DNA sequencing and statistical methodology have made it possible to investigate genetic divergence and the effects of environmental factors on the process of local genetic adaptation (Gougherty et al., 2020; Honjo & Kudoh, 2019). Environmental association analyses (EAA) have been gaining traction in the last decade (Ahrens et al., 2018), allowing identification of possible candidate genes involved in environmental adaptation from tens of thousands of genome‐wide single‐nucleotide polymorphisms (SNPs) sourced from populations across environmental gradients. For example, EAAs have been used to explore adaptive genetic variation on diverse and widespread woody plant genera, like Quercus (Gugger et al., 2021; Martins et al., 2018), Populus (Gougherty et al., 2021; Ingvarsson & Bernhardsson, 2020) and Corymbia (Ahrens et al., 2019). These studies have identified functional genes involved in adaptation to climatic factors that can be interpreted as divergent selection linked to population‐specific environmental variables (i.e., local adaptation to climate). However, different climate factors identify different sets of adaptive candidates, and few studies have focused on how sets of putatively adaptive SNPs sort across the landscape. If adaptive SNPs independently sort across the landscape (here, we define independent assortment as contrasting distributions of adaptation across the landscape), then understanding these species’ adaptative patterns to climate may prove to be difficult.

Identifying the genetic basis of local adaptation and selective environmental factors is still challenging. Genetic patterns that confer climate adaptations are mostly polygenic (Savolainen et al., 2013) and complex to investigate (Lind et al., 2018), particularly for species with limited genomic resources (Capblancq et al., 2020; Mayol, 2019). Genomic resources allow investigators to identify mutations that are more likely to affect adaptation. For example, nonsynonymous mutations in coding genes result in amino acid changes, which can yield changes in gene functions (Kryazhimskiy & Plotkin, 2008), or mutations in cis‐regulatory regions, which can often result in quicker adaptive processes (Wittkopp & Kalay, 2012). Beneficial mutations can be under selection among populations spread across that environment. Groups of genes found to be significantly associated with environment can be categorised into broader functional groups using gene ontology (GO) enrichment analysis (The Gene Ontology Consortium, 2019). GO terms have been used to predict polygenic adaptive biological processes and molecular functions associated with putatively adaptive SNPs in tree species (Collevatti et al., 2019; Jordan et al., 2017). However, few studies investigate how genes of similar function develop patterns of adaptation across complex landscapes. If genes with related functions are found to be adaptive to the same climate variable, this might be indicative of additive genetic variation controlling adaptation to the environment.

This study investigated the putative patterns of local adaptation associated with climate gradients across complex landscapes. To test hypotheses associated with signals of adaptation, we focused on Eucalyptus marginata Donn ex. Sm. (jarrah) because of its high genetic diversity and low population differentiation (Wheeler et al., 2003), and its ecological importance in the biodiverse hotspot of South‐west Western Australia (SWWA). This region has prolonged periods of extensive drying, with an estimated reduction of 20% in rainfall, from the 1970s to the present (Water Corporation, 2020), documented impacts of drought and heatwave events (Matusick et al., 2013), and the future (2030) climate is projected to show increased frequency and intensity of extremes (BOM & CSIRO, 2020). Furthermore, jarrah provenance trials have demonstrated genetic variation in functional traits associated to precipitation factors (Koch & Samsa, 2007; O'Brien et al., 2007), indicating potential local adaptation to drought stress. Ecological studies have also confirmed that water availability is critical for jarrah seedling survival and persistence (McChesney et al., 1995; Standish et al., 2015; Stoneman et al., 1994). Considering these studies on jarrah, we hypothesize that (1) populations show strong genetic patterns of local adaptation to climate, (2) precipitation is a stronger determinant of genetic adaptation compared to temperature, (3) functionally related genes show similar signatures of adaptation to climate, and (4) adaptive variants are independently sorting across the landscape. Lastly, we use this information to map the biological turnover of loci across the landscape to facilitate informed strategies for forest management that incorporate current patterns of genetic variation. We discuss how active management strategies, such as assisted gene migration (Aitken & Bemmels, 2016; Hoffmann et al., 2015; Prober et al., 2015) may incorporate the results from this study and be employed to build adaptive capacity to climate change.

MATERIALS AND METHODS

Sample collection and DNA extraction

Leaf samples from a total of 280 individual mature trees from 28 natural jarrah populations across the geographic range of the species (Figure 1), including one outlier population (JIL; located outside jarrah's natural range), were collected during 2019 (Table 1). Geographic coordinates were recorded for all sampled individual from each population using a handheld GPS device (Magellan eXplorist 310) (Table 1 shows a central point for each population). The sampling, which covered a total area of approximately 80,000 km2, included independent (>50 km separation) and replicate (across similar climate of origin) populations over both temperature and precipitation gradients to ensure adequate partitioning of the adaptive and neutral genetic variation. For each population, mature leaves were collected from 10 trees at least 100 m apart from each other. Leaves were immediately stored in silica gel until freeze‐dried using FreeZone 6 Liter Benchtop Freeze Dryer (Labconco Corporation), and kept in silica gel at room temperature until DNA extraction could be performed. For each sample, genomic DNA was extracted from 40 mg of freeze‐dried leaf material. Each leaf sample was independently ground into fine powder and a modified CTAB‐DNA extraction protocol was employed (Doyle & Doyle, 1990), with 0.1 M sodium sulphite (Byrne et al., 2001) and 2% w/v polyvinylpyrrolidone (MW 40,000) added to the extraction buffer. Quality of extracted DNA was estimated using gel electrophoresis and quantified using the Qubit dsDNA BR assay kit on a Qubit fluorometer (Invitrogen).

FIGURE 1

Sampling locations of jarrah in South‐west Western Australia (black squares). Two climate gradients are shown for the species distribution area: (a) maximum temperature of the warmest month, (°C; T MAX) and (b) mean annual precipitation (mm; P MA). Bioclimatic layers from worldclim.org (Fick & Hijmans, 2017). Insert shows distribution of jarrah in Australia

TABLE 1

Locations and climatic variables for the 28 sampled populations of jarrah in SWWA

Population	Code	Lat	Long	T SEAS	TMAX	TMIN	PMA	PWQ
Mt Lesueur	LES	‒30.1644	115.1991	41.1	32.2	8.2	578	35
Julimar	JUL	‒31.3491	116.2470	49.0	33.1	6.1	635	44
Jilakin Rock a	JIL	‒31.6647	118.3261	52.8	33.2	5.0	326	46
Chidlow	CHI	‒31.8622	116.2266	47.4	32.3	6.1	876	54
Perry Lakes	PER	‒31.9436	115.7838	37.6	30.4	9.4	765	38
Dale	DAL	‒32.1017	116.1900	45.9	31.5	6.2	1053	58
Serpentine	SER	‒32.3451	116.072	43.9	30.6	6.4	1151	57
Lupton	LUP	‒32.5292	116.5003	48.3	31.4	4.3	705	45
Whittaker	WHI	‒32.5499	116.0100	43.1	29.9	5.8	1190	62
Peel	PEE	‒32.6920	115.7103	37.5	30.4	8.3	888	42
Saddleback	SAD	‒32.9967	116.535	46.1	30.8	4.3	681	44
Godfrey	GOD	‒33.2142	116.5712	45.0	30.2	4.1	661	45
Yourdaming	YOU	‒33.3035	116.2407	43.9	30.4	4.1	851	46
Eaton	EAT	‒33.3177	115.7482	39.2	30.5	6.7	853	47
Meelup	MEE	‒33.5939	115.088	30.1	27.4	9.1	839	43
Grimwade	GRI	‒33.7612	115.9988	40.2	29.6	5.3	881	53
Katanning	KAT	‒33.8294	117.5731	41.9	29.2	5.2	457	50
Bramley	BRA	‒33.9035	115.0871	28.8	26.1	8.8	1072	54
Mowen	MOW	‒33.9133	115.5434	34.7	27.8	6.9	965	54
Nannup	NAN	‒33.9852	115.7778	36.1	28.3	6.6	928	56
Kingston	KIN	‒34.0825	116.3374	38.8	28	5.1	785	61
Milylannup	MIL	‒34.1928	115.6654	32.3	26.6	7.4	1027	64
Stirling Range	STI	‒34.3850	117.9927	35.4	26.9	5.8	493	67
Carey	CAR	‒34.4257	115.8223	30.6	26	7.6	1112	72
Boorara	BOO	‒34.6126	116.2060	31.4	25.9	6.9	1126	79
Plantagenet	PLA	‒34.6402	117.4987	33.7	26.7	6.5	738	79
Beadmore	BEA	‒34.8171	116.4834	31.3	25.8	7.0	1088	83
Denmark	DEN	201334.9535	117.3805	30.3	25.8	7.6	976	88

Temperature (T) and precipitation (P) variables are expressed in degrees Celsius (°C) and millimetres (mm), respectively.

Abbreviations: Lat, latitude; Long, longitude; P MA, mean annual precipitation; P WQ, mean precipitation of the warmest quarter; T MAX, mean maximum temperature of the warmest month; T MIN, mean minimum temperature of the coldest month; T SEAS, temperature seasonality.

Outlier population.

Genotyping by DArTseq platform

Sequencing of the 280 jarrah individuals was undertaken using DArT‐Seq technology at Diversity Arrays Technology Pty Ltd. This technology uses a double digestion complexity reduction method for next generation sequencing (Kilian et al., 2012). The reduction of the genome is accomplished by using a combination of PstI and HpaII enzymes in digestion/ligation reactions with different adapters corresponding to two different restriction‐enzyme overhangs. The PstI‐compatible adapter is designed to include flowcell attachment sequence, sequencing primer sequence and varying length barcode region. Diversity Arrays Technology's proprietary bioinformatic pipeline was used to demultiplex and align the raw fastq files. Identical sequences were then collapsed into fastqcall files. These files were used in the secondary pipeline for DArT P/Ls proprietary SNP calling algorithm (dartsoft14). Minimum read depth for each individual was set to six and average read depth was 30.93 across all SNPs, guaranteeing call quality for all SNPs and individuals. For the SNP calling algorithm, only nucleotide substitutions were considered a SNP. Only one random SNP was retained on each 75 base pair (bp) sequence to avoid linkage disequilibrium bias. All SNPs were mapped to the Eucalyptus grandis genome to obtain chromosome number and bp position to support the linkage disequilibrium analysis. The full data set was then filtered in r (R Core Development Team, 2020) using custom scripts. We applied a minor allele frequency filter of 2% (i.e., an allele frequency of 0.02), which equates to a minor allele count of at least 11 calls, minimising inclusion of sequencing errors. Missing data was set to 6% across individuals (SNPs were kept if they were called in at least 263 individuals). These thresholds were chosen because this translates to, on average, the presence of genetic information from nine individuals per population, which is adequate for EAA type of method and identifying SNPs under selection (Ahrens, Jordan, et al., 2021). Linkage disequilibrium (LD) was calculated within each of the chromosomes using the function LD. Measures in ldcorsv (Mangin et al., 2012). To guarantee adequate independence between SNPs and prevent potential linkage bias, the data set was filtered by the within chromosome pairwise LD r 2 coefficient (only one of the SNPs was randomly retained for analysis if the r 2 is >.5).

Environmental variables

Temperature and precipitation variables have been commonly assessed as predictors for environmental adaptation in eucalypts for phenotypic and genotypic variants (Aspinwall et al., 2019; Correia et al., 2018; Pritzkow et al., 2020). Climatic data for all populations was downloaded from the 19 variables in the worldclim v2 database (Fick & Hijmans, 2017) at a spatial resolution of 30 arcsec. Climate data for each population was extracted using the r package raster (Hijmans, 2021) from the geolocated GPS coordinates of the sampled populations. Principal component analysis (PCA) of environmental variables was performed with r package ade4 (Chessel et al., 2004) and a Pearson's correlation coefficient matrix was calculated between all 19 climate variables using the cor function. For later environmental association analysis, we tested a total of five climate variables. Two of these represent extreme temperature and precipitation variables and we predicted they would drive patterns of adaptation in MTC regions: maximum temperature of the warmest month (T MAX) and precipitation of the warmest quarter (P WQ). Three other temperature and precipitation variables were selected as independent climatic factors, based on PCA and Pearson's correlation coefficients (If the r > .7 of magnitude between two variables, these were considered strongly correlated and thereafter not included in further analysis; Table S2, tabs 16–17) and are known to be important for local adaptation in eucalypts (Queirós et al., 2020; Rocha et al., 2020): minimum temperature of the coldest month (T MIN), mean annual precipitation (P MA), and temperature seasonality (T SEAS).

Environmental association analysis

To understand how genetic structure of jarrah populations might affect EAA, genetic structure was estimated by measure of genetic differentiation (F ST) (Weir & Cockerham, 1984) using the hierfstat package (Goudet, 2005) in r. We also estimated individual ancestry coefficients for input for the EAA in LFMM. For this, we used the sparse nonnegative matrix factorization (SNMF) method in the r package LEA (Frichot & François, 2015). SNMF was run for each k‐value between 1 and 10, with each k‐value ran 10 times (200 iterations each). The optimum k‐value across all 10 runs was estimated using the software clumpp (Jakobsson & Rosenberg, 2007), and the graphical parameters were drawn in the program distruct (Rosenberg, 2004). The ideal k‐value was selected by visualising the cross entropies as defined in the SNMF manual (Frichot & François, 2015) and choosing the k‐value with the lowest cross entropy score for the LFMM analysis.

To elucidate the association between climate and genetic variation, three approaches were applied: a redundancy analysis (RDA), latent factor mixed models (LFMM) and BAYPASS. The three distinct approaches apply different statistical frameworks to identify population structure, which should reflect on the performance for each method (Lotterhos & Whitlock, 2015), but simultaneously allow a greater confidence in the identified associations (Ahrens et al., 2018). The detailed description for each method, including advantages and limitations, are broadly documented elsewhere (Ahrens et al., 2018; Hoban et al., 2016; Lotterhos & Whitlock, 2015; Rellstab et al., 2015). Therefore, we applied a multiapproach methodology as suggested by De Mita et al. (2013), providing a powerful detection of potential adaptive loci, regardless of limitations from each method (Rellstab et al., 2015).

Briefly, RDA is a multivariate method that assumes linear relationships from explanatory variables on response variables, thus allowing the estimation of genetic variance related to each distinct environmental factor simultaneously (Forester et al., 2018). RDA and LFMM require complete data sets, therefore we imputed missing data as the most common allele in the locus from the optimal ancestral cluster (k) as defined in the SNMF output. The explanatory variables (i.e., climate) were then constrained by the dependent variables (i.e., individuals), using the rda function in the vegan package 2.5‐1 in r (Oksanen et al., 2018). The anova.cca function was used to test for RDA significance using 999 permutations (randomised environmental variables). We did not explicitly control for population structure because RDA without explicit population structure inputs improves the output (Forester et al., 2018).

We also used LFMM to test for climate associations (Frichot et al., 2013), which applies a univariate regression model to assess genotype‐environment associations while using the optimal k‐value estimated in SNMF to control for ancestral population structure. This method is described as highly efficient to identify polygenic associations, even across diverse demographic sampling (De Mita et al., 2013; Lotterhos & Whitlock, 2015). The LFMM analyses were independently performed for each of the climate variables, consisting of 30,000 iterations each (15,000 discarded as initial burnin). Median z‐scores were combined from a total of five runs for each variable and recalibrated by manually adjusting the genomic inflation factor, λ, and then dividing the scores by λ. Adjusted p‐values were computed by flattening the histogram (false discoveries were controlled with the Benjamin‐Hochberg algorithm using q = 0.01), which ideally should display a peak close to zero. We used λ = 0.45 (this optimally flattened the histogram after testing other λ values as recommended in LFMM manual) in the adjustment function to flatten the histogram and followed the steps and R script available from the LFMM manual. To account for multiple comparisons, we applied a false discovery rate (FDR) threshold of 0.05 to all runs.

Lastly, we used a hierarchical clustering model implemented in baypass (Gautier, 2015), based on the model from bayenv (Coop et al., 2010). A population covariance matrix (Ω) was generated by running the core model. Each run had 100,000 iterations (50,000 discarded as initial burnin), repeated five times and averaged. The covariance matrix was then used in the AUX covariate mode (100,000 iterations; 50,000 as burnin), repeated five times and averaged for final results. Significant SNPs were identified if they had a Bayes Factor (BF) >3 (Kass & Raftery, 1995). Like LFMM, BAYPASS is based on a mixed linear model to account for potentially confounding allele frequency variances due to population structure. However, the difference between the two approaches may provide a means of identifying any influence of population structure (Ahrens, Jordan, et al., 2021; Forester et al., 2018).

Annotation and gene ontology analysis

To investigate the potential role of adaptive SNPs, identified by the three EAAs, in coding regions of genes, genomic annotation was run using the blastn function (Altschul et al., 1997) from BLAST (https://blast.ncbi.nlm.nih.gov/). The 75 bp sequences associated with each putatively adaptive SNP were annotated using the E. grandis genome (Myburg et al., 2014) and gene annotation considered significant if their blast output values met two related thresholds: an E‐value <1 × 10−8 and a blast‐score of at least 60.0. Chromosome number and location of SNPs were recorded, as well as annotated gene functions. Annotated genes were used in a comparative analysis with the co‐occurring tree marri (Corymbia calophylla) from Ahrens, Byrne, et al. (2019) to identify shared adaptive genes between the two species.

The putatively adaptive SNPs in genes were also used to predict broader biological functions using GO enrichment analysis through the web interface PlantRegMap (Tian et al., 2020). GO terms are organized within three categories: molecular function, cellular component and biological process. We explored the biological process aspect from the GO analysis, which refers to a category of broad processes accomplished by multiple genes or gene products. For each climate variable, Fisher's exact test was used to test for significantly overrepresented GO terms, with a threshold of p‐value <.01. Significant GO terms are defined by a set of genes associated with a specific climate variable. We use this output to explore how functionally related candidate SNPs associated to each climate variable are linked to GO terms, which might be additive for environmental responses (e.g., abiotic stress response) (Kulbaba et al., 2019). We only developed a subset of significant GO terms for landscape scale patterns to illustrate our point of variation among similarly functioning genes. We chose GO terms with the highest number of SNPs and/or with biological processes directly related to the climate variable (specifically response and/or tolerance to heat, cold and drought) for further landscape genomic analysis.

Landscape genomics

We used generalized dissimilarity modelling (GDM) to visualise the relationship between allele frequency of putatively adaptive SNPs and climate (Fitzpatrick & Keller, 2015). GDM is a statistical method that predicts spatial patterns of allelic turnover across geographic regions due to climate by generating an I‐spline turnover plot for each tested predictor and uses percent deviance explained as a measure of model fit (Ferrier et al., 2007). Specifically, the GDM spline plots show the association between predicted ecological distances and genetic dissimilarities (F ST matrix) while partitioning out variance explained by geographic distance; the y‐axis on the spline plots is therefore labelled as partial genetic distance, as it describes a portion of genetic distance, and the height of each spline indicates the magnitude of genomic turnover of a SNP along the climate gradient. GDM analyses was run using the gdm package v 1.3.7 in r (Manion et al., 2018), considering an input genotypic matrix (pairwise F ST for single putatively adaptive SNPs or SNP groups from GO terms) and a pairwise climate matrix that includes geographic coordinates. The hierfstat package in r was used to create population pairwise F ST matrices with the putatively adaptive SNPs associated with each climate variable. GDM was independently applied to all the putatively adaptive SNPs identified by the EAA as significant. For each climate variable, SNP GDM models with the highest value of deviance explained was selected for plotting and mapping of predicted allelic turnover to test the hypothesis of adaptive variants being independently sorted across the landscape.

Following the GDM transformation of the climate variables for each SNP, we performed PCA on the GDM transformed environmental layer using r. The PCA was predicted across geographic space using the predict function in the raster package and visualised using the plot function, revealing a 0–1 allelic turnover map. To test the hypothesis of additive variation, we ran GDM analyses on groups of SNPs related to specific GO terms for each of the five climate variables and visualised how the allelic turnover within the GO term was related to that climate. By running GDM analysis on groups of SNPs, we were able to quantify an “additive score” of deviance explained for each set of SNPs to compare the importance of GO terms. To test the independent assortment hypothesis among SNPs, we directly compared maps derived from the GDM models that explained the most deviance for each of the five bioclimatic factors using a Spearman's correlation coefficient test in the spatialeco package (Evans, 2021).

RESULTS

Sequencing and population structure

A total of 78,198 SNPs were generated by DArTseq technologies and filtered down to 13,534 independent SNPs, with 8824 SNPs mapped to the 11 Eucalyptus grandis chromosomes. The number of SNPs per chromosome varied from 599 to 1083, with a mean of 802 SNPs per chromosome. Of the remaining SNPs, 477 fell on unspecified scaffolds and 4233 on regions that could not be aligned to the E. grandis genome (unknown location). Population differentiation was low (global F ST = 0.04), and similar to that identified in a previous RFLP analysis of variation (F ST = 0.034; Wheeler et al., 2003), with population pairwise F ST values ranging from 0.011 to 0.18 (Table S2, table 18). The cross‐entropy analysis estimated that the optimal number of clusters (k‐value with the lowest entropy score) was 6 (Figure S1). SNMF analysis with six clusters revealed substantial admixture in populations. Five of the clusters could be geographically described (Figures 2 and S2), one cluster was primarily located in the southern area, one in the central area and two in the northern area, where one cluster was dominant in populations along the coast. A fifth cluster occurred in the outlier population (JIL; blue colour), and the sixth cluster was present in four individuals (two individuals from BRA and BOO). The LES population, which is the northmost population, displays mixed affinity, being similar to both southern (green) and northern (yellow and red) populations. This is consistent with the northern areas harbouring ancestral variation in other co‐occurring species (e.g., marri; Sampson et al., 2018), or could indicate possible historic human influence through Aboriginal movement of plants (Lullfitz et al., 2020; Lullfitz, Dabb, et al., 2020).

FIGURE 2

Distribution of sampled jarrah displaying population membership proportion for K = 6 genetic clusters, depicted as pie charts. Table 1 provides more details on each population

The full set of 13,534 independent SNPs were considered for environmental association analysis. All three EAA approaches found putatively adaptive SNPs for each of the five climate variables (Table S2, tabs 1–5). The RDA approach identified fewer putatively adaptive SNPs than BAYPASS and LFMM that identified similar numbers (Figure 3; Table S1). The proportion of overlapped SNPs is different for each variable (Figure 3). Overall, 2336 unique SNPs were flagged to be associated with at least one of the tested climate variables across the three EAA approaches. RDA analysis (Figure S3) identified between 16 (T SEAS) and 57 (T MAX) SNPs significantly associated with each of the climate variables, for a total of 168 SNPs. All five climate variables were shown to be significantly associated with variation in the RDA (T SEAS: F = 3.98, p = .001; T MAX: F = 2.26, p = .001; T MIN: F = 2.11, p = .001; P MA: F = 1.80 p = .001; P WQ: F = 1.49, p = .001). LFMM identified between 263 (P MA) and 411 (T MAX) SNPs with significant correlations, with a total of 1753 putatively adaptive SNPs. BAYPASS identified between 284 (T MIN) and 888 (T MAX) SNPs with significant correlations, with a total of 2327 putatively adaptive SNPs. Putatively adaptive SNPs found for all environmental variables from each EAA method were used in further analyses.

FIGURE 3

Summary of environmental association analysis in jarrah. Venn diagrams show the intersections between the three approaches of environmental association analyses (RDA, red; LFMM, blue; BAYPASS, yellow) considering adaptive SNPs associated with each of the climate variables: P MA, mean annual precipitation; P WQ, mean precipitation of the warmest quarter; T MAX, mean maximum temperature of the warmest month; T MIN, mean minimum temperature of the coldest month; T SEAS, temperature seasonality

Annotation and gene ontology

Of the 2336 unique putatively adaptive SNPs associated with the climate variables, 1440 SNPs were linked to functionally annotated genes, which represents 10.6% of the total independent SNPs tested for EAA (13,534). Full annotation results for SNPs associated with each variable are given in Table S2, tabs 6–10. T MAX delivered the highest amount of linked functionally annotated genes (474), followed by P WQ (312), T SEAS (237), P MA (214) and T MIN (203). Annotated SNPs (Table 2) included, for example, JAR00198, associated with both T SEAS and T MIN, located in a transcinnamate 4‐monooxygenase (TCMO) gene; two SNPS associated with T MAX, JAR00038 and JAR00207 were found on transcription repressor MYB6 and transcription factor MYB44 genes respectively; P MA, JAR02395, located in a peroxidase 72 gene.

TABLE 2

Gene annotation showing the top five SNPs (Blast score, 125) for jarrah, with NCBI blast e‐value score, relative ranks based on levels of significance for each EAA and chromosome number (chr) from Eucalyptus grandis genomic mapping (un, unknown) for each environmental variable

Climate	SNP	RDA (p‐value)	LFMM (p‐value)	BAYPASS (BF)	chr	Blast e‐value	Gene annotation from the Eucalyptus grandis genome
T SEAS	JAR00166	–	.00034	4.708	un	1.0E‐28	Mitochondrion
	JAR00198	–	.00064	6.788	10	1.0E‐28	Trans‐cinnamate 4‐monooxygenase
	JAR00273	–	3.11E‐05	21.663	11	1.0E‐28	Mitochondrion
	JAR00499	–	.00071	–	8	1.0E‐28	Probable LRR receptor‐like serine/threonine‐protein kinase
	JAR00662	–	4.78E‐06	–	6	1.0E‐28	UPF0496 protein
T MAX	JAR00038	–	–	3.270	6	8.0E‐29	Transcription repressor MYB6
	JAR00207	–	–	3.054	6	8.0E‐29	Transcription factor MYB44
	JAR00209	–	–	9.466	11	8.0E‐29	AT‐hook motif nuclear‐localized protein 16
	JAR00214	–	–	11.262	6	8.0E‐29	Protein indeterminate‐domain 1
	JAR00262	–	–	6.154	4	8.0E‐29	Uncharacterized
T MIN	JAR00013	–	–	9.801	10	8.0E‐29	Mitochondrion
	JAR00166	–	–	18.303	un	1E‐28	Mitochondrion
	JAR00198	–	.00026	6.788	10	1E‐28	Trans‐cinnamate 4‐monooxygenase
	JAR00273	–	1.32E‐06	–	11	1E‐28	Mitochondrion
	JAR00620	0.242	–	–	11	8.0E‐29	Uncharacterized
P MA	JAR00027	–	.00014	10.316	7	1.0E‐28	Mitochondrion
	JAR00500	–	–	6.788	4	1.0E‐28	Putative yippee‐like protein Os10g0369500
	JAR01426	–	.0004	–	11	1.0E‐28	Tyrosine decarboxylase 1
	JAR01512	–	.0001	–	5	1.0E‐28	Uncharacterized
	JAR02395	–	.00092	–	9	1.0E‐28	Peroxidase 72
P WQ	JAR00214	0.454	.00053	–	6	8E‐29	Protein indeterminate‐domain 1
	JAR00273	–	–	11.889	11	8E‐29	10 kDa chaperonin
	JAR00499	–	.00021	–	8	1E‐28	Probable LRR receptor‐like serine/threonine‐protein kinase A
	JAR00690	–	.00031	6.266	1	8E‐29	Zinc finger protein ZAT5
	JAR01091	–	5.34E‐07	4.388	7	1E‐28	LOB domain‐containing protein 1‐like

SNPs that were also found associated with GO terms (Table 3) are in bold.

Abbreviations: P MA, mean annual precipitation; P WQ, mean precipitation of the warmest quarter; T MAX, mean maximum temperature of the warmest month; T MIN, mean minimum temperature of the coldest month; T SEAS, temperature seasonality.

TABLE 3

Overrepresented gene ontology (GO) terms for SNPs identified in jarrah for each environmental variable

Climate	GO id	GO term	GDM % Deviance explained	p‐Value	SNPs	Count
T SEAS	GO:0009314	Response to light stimulus	3.66	.00261	JAR02551, JAR04603, JAR00284, JAR02659, JAR06621, JAR04257, JAR07363, JAR00198, JAR01133, JAR07395	10
T MAX	GO:0000271	Polysaccharide biosynthetic process	4.61	.0067	JAR05227, JAR06489, JAR06314, JAR08046, JAR12549, JAR11847, JAR08134, JAR12439	8
	GO:0010104	Regulation of ethylene‐activated signalling pathway	2.80	.0097	JAR09402, JAR12137	2
T MIN	GO:0071840	Cellular component organization or biogenesis	2.80	.00859	JAR05151, JAR02381, JAR00166, JAR02528, JAR04603, JAR03088, JAR05858, JAR01284, JAR06869, JAR05668, JAR04700, JAR07286	12
	GO:0080167	Response to karrikin	14.27	.00391	JAR06869, JAR00198, JAR03623	3
P MA	GO:0044763	Single‐organism cellular process	3.01	.0061	JAR07368, JAR04995, JAR05607, JAR07695, JAR05954, JAR12984, JAR11156, JAR00284, JAR11454, JAR08984, JAR13223, JAR06091, JAR07363, JAR08184, JAR12280	15
	GO:0009411	Response to UV	4.62	.0059	JAR00284, JAR07363	2
P WQ	GO:1901566	Organonitrogen compound biosynthetic process	21.22	.0064	JAR12369, JAR00189, JAR00543, JAR11122, JAR05879, JAR12666, JAR02347, JAR11253, JAR11737, JAR06747, JAR07829, JAR06097, JAR10308, JAR12789, JAR12316, JAR00476, JAR13196, JAR11797, JAR11414, JAR11170, JAR10452	21

The GDM % deviance explained is expressed as additive score for groups of SNPs.

Abbreviations: P MA, mean annual precipitation; P WQ, mean precipitation of the warmest quarter by count of SNPs and/or relevant biological function; T MAX, mean maximum temperature of the warmest month; T MIN, mean minimum temperature of the coldest month; T SEAS, temperature seasonality.

Gene ontology enrichment analysis explored how groups of annotated SNPs relate to similar functions (Table S2, tabs 11–15). Enriched GO terms in the biological process category are highlighted for each of the five bioclimatic variables (Table 3). A GO term associated with response to light stimulus (GO:0009416) was found with the SNPs related to T SEAS. Genes associated with this GO term are linked to cellular response processes (in terms of components movement, enzyme production, and secretion and protein expression) from abiotic stimulus, specifically electromagnetic radiation and light. A GO term related to karrikin stimulus was found associated with T MIN (GO:0080167). As for P MA and P WQ, GO terms with high counts of SNPs were found for each variable (GO:0044763 and GO:1901566, respectively) as well as a term related to UV response (GO:0009411) associated with P MA.

Landscape modelling

We independently applied GDM analysis to all putatively adaptive SNPs associated with the five climatic variables (Figure S4), and the models that explained the highest deviance for each variable were selected to display spatial patterns of allelic turnover (Figure 4): T SEAS – JAR00269 (39.2%); T MAX – JAR11943 (25.5%); T MIN – JAR01172 (16.8%); P MA – JAR10596 (21.9%) and P WQ – JAR06621 (36.9%). The GDM model for the SNP associated with T SEAS (JAR00269) explained more deviance than any other putatively adaptive SNP across the five climate variables, followed by the model for P WQ (JAR06621). There was rapid turnover noticeable for the three temperature variables from the coastal to eastern populations in the north of the range, and more gradual turnover from the northern populations to the southern populations (Figure 4a–c). But even among the three temperature variables, there were major differences in adaptive patterns. For instance, while T SEAS and TMAX displayed a similar rapid turnover from the coastal to eastern populations in the north of the range, and fairly gradual turnover from the northern populations to the southern populations, T MIN followed the same trend in the northern region, but a rapid turnover is present between the coastal and inland populations in the south region. In contrast, the precipitation variables showed rapid turnover in the southern or central parts of the distribution, and more gradual turnover in the northern distribution (Figure 4d,e). In southern areas, P WQ showed a rapid turnover between coastal and inland southern populations, while P MA showed a more gradual pattern in this region. Correlation coefficients between allelic turnover maps showed clear differences between some of the adaptive landscapes (Figure 4; bottom right table, below the diagonal), as r 2 values ranged between the absolute values of .21 to .74, and often different than correlation between their respective climate variables (above the diagonal). These differences are indicative of contrasting spatial patterns.

FIGURE 4

Predicted spatial variation of allelic turnover based on the outputs from the GDM models that explained the most deviance for each climate variable (between climate and SNP) for jarrah. (a) T SEAS, JAR00269; (b) T MAX, JAR11943; (c) T MIN, JAR01172; (d) P MA, JAR10596 and (e) P WQ, JAR06621. Insets are spline plots of partial genetic distance (y‐axis) by climatic distance (x‐axis) for the individual SNP (dimensions of the plot are the same as in Figure 4). Table shows the pairwise Spearman's correlation coefficient (r 2) between the two allelic turnover maps (below diagonal) and between the climate variables (above the diagonal)

The groups of SNPs associated with selected GO terms (Table 3) were also used in a combined GDM analysis to measure allelic turnover across climatic gradients and interpreted as an additive pattern of adaptation (Figure 5). The patterns of allelic turnover varied by climatic variable: overall, GDM showed small to moderate response, in terms of deviance explained. The GDM model explained more deviance for the group of SNPs linked to the GO term associated with P WQ (n = 21; 21.22%, Figure 5e) than any other climate variable association using GO terms, followed by GO terms associated with T MIN (n = 3; 14.27%, Figure 5c). Models for T SEAS, T MAX and P MA explained a similar deviance for allelic turnover composition (<5% for each group of SNPs).

FIGURE 5

Geographic generalized dissimilarity modelling (GDM) in jarrah showing SNPs allelic turnover for gene functions (GO terms) across each environmental variable. P MA, mean annual precipitation; P WQ, mean precipitation of the warmest quarter; T MAX, mean maximum temperature of the warmest month; T MIN, mean minimum temperature of the coldest month; T SEAS, temperature seasonality. GO terms with different SNP sets in the same plot are represented with different colours (black or orange)

DISCUSSION

Our study identified putative patterns of climate adaptation in jarrah, with several strong associations between candidate SNPs and climatic gradients. The results provide support for our hypothesis of strong patterns of local adaptation to climate across the distribution of jarrah, although, contrary to our second hypothesis, we found adaptation to both temperature and precipitation variables rather than primarily with precipitation. As expected, annotation highlighted functional genes associated with biological processes, some of which relate to abiotic stress factors and provide good candidates for adaptations. Furthermore, the landscape genomics modelling assessed the magnitude of allelic turnover for putatively adaptive SNPs and highlighted temperature seasonality, mean maximum temperature of the warmest month and precipitation of the warmest quarter as explaining significantly more variation than other climate drivers. These patterns indicate that adaptive variants are independently sorting across the landscape, which is consistent with our fourth hypothesis. We discuss the mechanisms for adaptation to climate across complex landscape for forest trees, including a direct comparison with a codominant co‐occurring foundation species, before providing the scientific basis for implementation of management and conservation strategies to promote the resilience of foundation tree species. These associations, indicating potential local adaptation, were found despite high levels of gene flow among populations across the distribution, a common characteristic among eucalypt species (Ahrens, Byrne, et al., 2019; Jones et al., 2002; Murray et al., 2019; Supple et al., 2018). Low differentiation among populations indicates that application of EAA in jarrah is appropriate to identify alleles putatively under selection. As expected, the three EAA approaches identified different sets of putatively adaptive SNPs for each climate variable. One limitation of EEAs is the identification of SNPs that are found to be under selection but are in fact not (false positives). While false positives are an inherent limitation in EAA studies, EAAs can consistently identify adaptive SNPs, even if the selection coefficient is small (Ahrens, Jordan, et al., 2021). We focus our interpretation on SNPs that are within gene space to lessen the impact of false positives, despite that candidate SNPs identified outside of gene space could be true positives. For instance, SNPs could be in promoter regions, which are known to have high proportion of adaptive variants (Wittkopp & Kalay, 2012), SNPs could form a large haplotype block with genes that are under selection (Todesco et al., 2020), or SNPs could be in linkage disequilibrium with adaptive SNPs. Future work should focus on improving the genomic resources of the species to elucidate these complex issues that are beyond the scope of this work.

Adaptation to temperature and precipitation

Generalized dissimilarity modelling models on all putatively adaptive SNPs found the highest deviance explained for a SNP associated with T SEAS (39.2%), closely followed by P WQ (36.9%), with overall results showing low to moderate deviance across the five climate variables. Furthermore, P WQ was linked to GO:1901566, with the highest number of putatively adaptive SNPs (21) and also showed the highest deviance explained by the GDM analysis (Figure 5e). Overall, both temperature and precipitation variables are linked to adaptive genetic variants; although, GO and GDM analyses highlighted the specific precipitation variable (P WQ) as a stronger driver of adaptation.

The annotations of putatively adaptive genes were made based on the reference genome of E. grandis, a distant relative, so we provide a pertinent but cautious preliminary interpretation of functional results until a full jarrah reference genome becomes available. Gene functions associated with the temperature and precipitation variables show biological functions associated with response and adaptation to these abiotic factors. For example, the KCS gene family (JAR02659), that was associated with T SEAS, has been linked to cold and light responses (Joubès et al., 2008) in Arabidopsis, being involved in the biosynthesis of waxes that cover the leaves surface. Two SNPs, also associated to T SEAS, (JAR13256 and JAR08936) are linked to the ABC transporter gene families, which have been shown to be associated with heat response and abiotic stress tolerance during seed germination (Hwang et al., 2016; Zhang et al., 2012). Control of seed germination during periods of thermal stress could be a crucial mechanism for selecting phenotypes that are more adapted to Mediterranean type of climates with hot and dry summers. In that sense, we can envision how genes that control seed germination response could be selected for or against in such climates (Rix et al., 2015). A SNP associated with P WQ (JAR13490) was found in the chromatin‐remodelling factor PKL gene that has been consistently linked to multiple plant development processes, particularly to the abscisic acid (ABA) pathway regulation (Perruc et al., 2007). ABA is a phytohormone that is well known for controlling stomatal closure (Maheshwari et al., 2020; Rajab et al., 2019), thus being crucial for efficient drought response (Yu et al., 2019; Zhang et al., 2020). Drought is a known selective force in the region, and the ability to control stomates would also allow for the reduction of transpiration during dry periods, a physiological ability that has been shown to be crucial in the co‐occurring marri (Challis et al., 2020). These are just a sample of the many compelling gene functions associated to either temperature or precipitation found across the five tested climatic variables, identifying these variables as potential drivers of local adaptation.

Functionally related genes have similar adaptive patterns

In our GO enrichment analysis, we focused on biological processes related to abiotic stress responses such as drought, cold and heat. Generally, we found biological GO terms with gene overrepresentations, consistent with expectations under our third hypothesis: functionally related genes have similar patterns of correlation with climate. For instance, a GO term related to karrikin stimulus was found associated with T MIN (GO:0080167). Karrikins are a group of phytohormones that control several aspects of plant germination and growth (Nelson et al., 2012). A study with Arabidopsis showed that karrikin signalling can inhibit seed germination under heat stress (Wang et al., 2018), possibly to avoid germination under conditions unfavourable to seedling establishment.

Many plant functional traits are polygenic, involving complex interactions controlled by multiple genes, so it is expected that patterns of climate adaptation are also the result of combined effects from several alleles of small‐effect (Wadgymar et al., 2017). Indeed, climatic variables are expected to not be the main driver for variation in some putatively adaptive SNPs, as the genes associated can be pleiotropic and may be under selection from other biotic or abiotic factors. For example, although precipitation and temperature are consistently highlighted as key factors influencing plant distribution and ecology, soil properties greatly affect these settings, as water availability depends on the interaction between climatic variables and soil characteristics (Piedallu et al., 2013). The identification and understanding of adaptive genetic variations might then be improved by including other relevant abiotic factors such as soil characteristics. Nevertheless, by hierarchically categorising gene functions through GO enrichment, we were able to find adaptive patterns across the distribution, highlighting likely polygenic adaptations to climate variables in this species.

Adaptive variants are independently sorted

Across the species geographic distribution, climatic heterogeneity explains significant genomic variation. In particular, two distinct climate variables, T SEAS and P WQ, showed strong associations with genomic variants. The patterns of genomic turnover associated with the studied climatic variables are aligned with the climatic gradients of the region (Figure 4). These associations are indicative of the multidimensional patterns of adaptation resulting in uncorrelated intraspecific selection among SNPs (White & Butlin, 2021). Here, we define dimensionality as the interaction between uncorrelated climate variables to independently describe each habitat. Our dimensionality is driven by climate, and the independent sorting of putatively adaptive SNPs is indicative of this complex pattern. It has been modelled that local adaptation increases with dimensionality (MacPherson et al., 2015), and it probably leads to dimensionality of phenotypic traits (Kirkpatrick & Meyer, 2004; McGuigan et al., 2005). Indeed, there is evidence of intraspecific variation among growth traits (e.g., height and diameter) that are locally adapted in jarrah (Koch & Samsa, 2007; O’Brien & Krauss, 2010; O’Brien et al., 2007).

In some ways, increased dimensionality is ubiquitous with increased habitat heterogeneity, and habitat heterogeneity has been shown to drive signatures of adaptation to temperature and precipitation in tree species (Shryock et al., 2020; von Takach et al., 2021; Walters et al., 2021). While these studies did not explore dimensionality explicitly, their results nevertheless show that tree species are able to independently adapt to multiple types of environments. While such patterns of differential adaptation make management of the species more complex and nuanced in the future, our results provide a level of understanding that will allow for targeted responses to changing climatic conditions in different regions.

Landscape adaptations of forests

Comparative analysis can provide broader patterns for forest management, where concurrent genetic and spatial patterns of local adaptation within co‐occurring tree species provides strong evidence for environmental fitness and evolution (Bragg et al., 2015). Our analysis here identified SNPs associated with both temperature and precipitation in jarrah; while a similar study on a co‐occurring species, marri, found SNPs associated with temperature to explain more deviance than precipitation (Ahrens, Byrne, et al., 2019), thereby suggesting that temperature is a stronger driver of local adaptation for marri. It is interesting that there were similarities in functional genes associated with several adaptive variants between jarrah and marri (e.g., ABC transporters and CBL gene families). Comparison across both species identified a set of 26 genes that were also found to be associated with at least one of the five variables analysed (Table S3). Most of these shared genes are associated with either T MAX (16) or P WQ (12) in jarrah; while for marri, the majority of the shared genes are associated with T MAX (24), which is consistent with adaptation to both temperature and precipitation in jarrah and with temperature in marri. While this comparison shows that some functional genes share adaptive patterns, there were more genes that were different, indicating that the same adaptive management plan may not be effective for both species.

Management perspectives

Our analysis of standing genetic variation across the distribution of jarrah found potential links between putatively adaptive SNPs and climate factors, which may provide a source of adaptation to future climate conditions. The evidence that genetic variants are involved with climate adaptation occurred as either associations with specific annotated gene functions or biological processes associated to climate factors. Our analysis here, and that of the co‐dominant species marri (Ahrens, Byrne, et al., 2019), are also consistent with results from recent genomic studies on other eucalypt species in other regions of Australia (Jordan et al., 2017, 2020; Steane et al., 2017), providing evidence of adaptation to climate in natural populations and stressing the role of temperature (particularly T SEAS and T MAX) and precipitation (P WQ) variables. The presence of climate adaptation provides a basis for implementation of assisted gene migration for forest management strategies (Aitken & Bemmels, 2016) and climate adjusted provenance (i.e., sourcing of seed from populations in the direction of climate change for use in restoration sites to enhance adaptation to future climate) in restoration practices (Prober et al., 2015). As a foundation tree, jarrah is a vital component in the ecosystem and has a significant role in regulating local hydrological systems and carbon storage (Bradshaw, 2015; CCWA, 2013). Additionally, it offers abundant habitats for a wide variety of groups, from vascular flora and lichens to terrestrial vertebrates and birds (Whitford & Williams, 2002; Whitford et al., 2015), as well as unique food sources for fauna, especially birds (Lee et al., 2013; Wrigley & Fagg, 2012). The Forest Management Plan 2014–2023 (CCWA, 2013) for SWWA forests has provision for implementation of assisted gene migration in management strategies for response to climate change. Our findings of standing variation harbouring putative adaptations to climate associated with temperature and precipitation factors provides an evidence base for design and implementation of such strategies. In addition, phenotypic approaches on other eucalypt species have also highlighted the role of local climate in the development of adaptive traits (Ahrens, Andrew, et al., 2019; Ahrens, Rymer, et al., 2021; Costa e Silva et al., 2019). Expanding this work to a phenotypic approach in jarrah (such as O’Brien et al., 2007) for identifying patterns of plasticity and adaptation associated with climate would contribute to further understanding the association of genomic and phenotypic diversity across environmental gradients. While it appears that genetic variants associated with similarly functioning genes are adapting to the environment in similar ways, we also found that putative adaptations among climate variables are sorted through the landscape in contrasting ways. This makes implementation of assisted gene migration strategies more complex and multidimensional. In fact, our findings support recommendations for sourcing germplasm from multiple sources to bolster the adaptability in adaptively depauperate populations and provide a basis for more active selection of functionally related genes, potentially increasing the diversity and adaptability through new combinations of genetic variation.

CONFLICT OF INTERESTS

The authors declare that they have no conflict of interests.

AUTHOR CONTRIBUTIONS

Collin W. Ahrens, Paul D. Rymer, and Margaret Byrne conceptually designed the experiment. João Carlos Filipe, Collin W. Ahrens, Paul D. Rymer, and Margaret Byrne designed the sampling strategy. João Carlos Filipe and Richard Mazanec collected and prepared all samples for sequencing. João Carlos Filipe and Collin W. Ahrens analysed and interpreted the data, with outputs interpreted by João Carlos Filipe, Collin W. Ahrens, Paul D. Rymer and Margaret Byrne. João Carlos Filipe and Collin W. Ahrens wrote the manuscript with intellectual contributions from Paul D. Rymer and Margaret Byrne, and all authors edited the manuscript.

Supplementary Material

Click here for additional data file.

Table S2

Click here for additional data file.