Introduction

Incorporating genetic, environmental, and ecological data can enhance our understanding of how ecology, evolution, and geography contribute to the processes that drive spatial distribution (Balkenhol et al. 2017). Exploring how environmental factors influence population genetic differentiation can also improve our knowledge of the elements responsible for speciation, which are central to the diversification of species (Barraclough and Nee 2001; Manel and Holderegger 2013). Species with narrow geographic distributions within heterogeneous environments represent ideal examples with which to examine the function of environmental factors in shaping genotypic variation. In this context, the rainforests of central and eastern Australia provide the opportunity to explore these elements in one of the most complex and diverse subtropical biomes in the world (Chester and Bushnell 2005).

Fontainea Heckel (Euphorbiaceae) is a genus of dioecious rainforest shrubs and small trees with links to the Gondwana landmass (Brunton et al. 2023) that is patchily distributed among rainforest fragments of Australia, Papua New Guinea, New Caledonia, and Vanuatu. Several Fontainea species have low to moderate genetic diversity and occupy restricted environments as a result of historic rainforest contraction and expansion cycles, coupled with more contemporary effects of habitat loss, modification, and fragmentation from human disturbance (Lamont et al. 2016; Conroy et al. 2019; Brunton et al. 2022). The critically endangered F. oraria (Coastal Fontainea) and vulnerable F. australis (Southern Fontainea) belong to a closely related subclade that occurs across the heterogeneous landscape of the Wollumbin Caldera in central, eastern Australia (Act 1999; Brunton et al. 2022). Similar to most Fontainea species, Southern Fontainea occur in small, clumped, and isolated populations (Forster, 1997), with recent surveys estimating there are < 200 mature, wild individuals of this taxon from 13 known populations (Brunton et al. 2023). Coastal Fontainea is presently only found on the outer, southeastern margin of the Wollumbin caldera within several small scattered littoral (coastal) rainforest fragments on basaltic soils (Brown et al. 2016) derived from the Wollumbin volcano lava flows in a 600 m radius. With just 10 mature individuals, it is considered at-risk of extinction, which has prompted in-situ and ex-situ conservation management to enhance the long-term survival of the species (Rossetto et al. 2000b). Rossetto et al. (2000b) and, more recently, Brunton et al. (2022) showed that Southern Fontainea and Coastal Fontainea have shared plastid genotypes and that morphological differences between species may represent the influence of distribution across the highly variable environment of the Wollumbin Caldera. Another hypothesis is that the two species are at various points along a speciation continuum as a result of geographic isolation and differing environments (Brunton et al. 2022).

Local adaptation to specialised habitats can be a key driver of genetic differentiation between species (Holderegger et al. 2006; Li et al. 2017) and can develop in populations in response to environmental factors such as climate and topography (Kooyman et al. 2012; Lajoie et al. 2018). Such responses can result in changes to the genetic structure of populations over multiple generations, leading to Isolation by Environment (IBE; Wang and Bradburd 2014). IBE is the relationship between genetic divergence among populations and environmental differentiation, outside of geographic distance. Genetic differentiation can also be related to geographic distance due to limited gene flow between distant populations in the phenomenon known as Isolation by Distance (IBD; Wright 1943). For example, Elaeocarpus sedentarius, a species endemic to the Wollumbin caldera, has a restricted range as a result of dispersal limitations rather than highly specialised habitat requirements, although fine-scale environmental factors do affect different populations of the closely related E. grandis (Rossetto et al. 2008). By contrast, using reduced-representation genome sequencing, Rutherford et al. (2018) showed that gene flow was not a strict constraint to population differentiation for six closely related Eucalypt species. Overall, the study showed that the various Eucalypt lineages were in different stages along a continuum of speciation, indicating a number of mechanisms such as vicariance, divergent selection between environments (ecological speciation), and origins from partially merged lineages (reticulate evolution) acting in combination were likely responsible.

Recent evidence suggests that Southern Fontainea and Coastal Fontainea are poorly conserved, with most populations located on private land, and are at risk of decline from limited habitat suitability under future climate projections (Brunton et al. 2023). Therefore, investigating the influence of spatial drivers of population genetic differentiation is crucial to understanding how the effects of geographical and environmental isolation have shaped the current distribution and genetic composition of these threatened Fontainea species. This will aid in identifying geographical barriers to gene flow and the role that the local environment has played in the genetic divergence of their populations. Furthermore, exploring genome-wide data for evidence of loci correlated with local adaptation, and identifying the processes of selection, can ultimately be used to inform conservation management decisions to preserve populations of Southern Fontainea and Coastal Fontainea and to guide potential genetic rescue efforts.

Using reduced-representation single-nucleotide polymorphism (SNP) markers, we utilised a multiple matrix regression approach to examine the contributions of geographic and environmental distance to spatial genetic variation of the sister species, Southern Fontainea and Coastal Fontainea. As such, we aimed to disentangle the individual and combined effects of IBD, IBE on the patterns of genetic differentiation both within and among the two Fontainea species. In addition, we explored genotype-environment interactions by applying latent factor mixed models (LFMM) and a multivariate partial redundancy analysis (pRDA) to identify putative loci under selection and identify key ecological drivers of population differentiation between these two closely related, threatened rainforest species.

Methods

Study system and field sampling

The Wollumbin caldera is the largest erosion caldera in the southern hemisphere, covering an area of approximately 2,400 km2, with its northern extent spanning the border of Queensland (QLD) and New South Wales (NSW), central eastern Australia (Solomon 1964). Intermittent lava flows from the eruption cycles of Mount Warning (Wollumbin) ~ 23 million years ago, have generated a banded geology that has undergone variable weathering to form a complex, geodiverse landscape with a unique range of landforms in a circular arrangement around the central mass (Willmott and Blight 2004). This system is thus characterised by distinct habitat shifts marked by steep slopes, cliffs, gullies, and ridges, which influence rainfall, temperature, and topographical gradients at a local and landscape scale. Essentially, the caldera represents an area with distinct physical boundaries and microclimates formed by natural processes and the more contemporary effects of land-clearing, which has created the isolated, patchy distribution of subtropical rainforest in the region (Laidlaw et al. 2016).

Here, we focussed our analyses on the two closely related sister species (Brunton et al. 2022) Southern Fontainea (Fig. 1a) and Coastal Fontainea (Fig. 1b), endemic to the rainforest fragments across the Wollumbin caldera in central, eastern Australia (Fig. 1c). Both species only occur in a small number of isolated and scattered populations (Rossetto et al. 2000b). Small population size can limit innovative conservation methods, yet this is the demographic reality for many rare, threatened species that require particular attention to ensure their long-term survival (Frankham 2005). Genetic evidence indicates limited genetic distance between some populations of these two geographically proximate and phenotypically similar Fontainea species. In addition, phenotypic distinctions between the species are expected to be driven by local adaptation to variable habitats (Rossetto et al. 2000b; Brunton et al. 2022). In order to further understand the genetic distinctions between Southern Fontainea and Coastal Fontainea, samples from across the broader species’ range are necessary for a more complete representation of the genetic profile for each respective taxon. Therefore, leaf samples were collected from 139 Southern Fontainea individuals from 12 sites and seven mature individuals of Coastal Fontainea throughout the natural distribution of the two species across the Wollumbin Caldera between the years 2020 and 2021 (Supplementary Table 1 and Fig. 1c). Leaf samples were dried with silica in zip-lock plastic bags in preparation for genotyping.

Fig. 1
figure 1

Distribution and images of two sister species of threatened Fontainea; a Fontainea australis (Southern Fontainea), with inset image of fruit and; b F. oraria (Coastal Fontainea), including male flowers and inset image of fruit and; c satellite image of the Wollumbin caldera region and sampling locations throughout the distribution of the two study species in central, eastern Australia. Inset maps indicate relative area (in orange) to larger map extents

Full size image

Reduced-representation genotyping and SNP data filtering

Silica-dried leaf samples were genotyped by Diversity Arrays Technology Pty Ltd (DArT, Canberra, Australia), who performed DNA extraction and generated reduced-representation SNP markers using a proprietary restriction enzyme and Next Generation Sequencing method, known as DArTseq (Kilian et al. 2012). DArTseq has been used for a variety of plant phylogenetic and population studies including commercial and natural populations (Steane et al. 2011; Iqbal et al. 2021; Brunton et al. 2022). Quality filtering steps on the raw DArTseq SNP markers was implemented by an in-house designed R package RRtools v.1.0 (Rossetto et al. 2019) to remove poor quality SNPs that have a reproducibility score of < 0.96, a call rate of ≥ 0.95, missingness of < 20%. Low-quality samples (high proportion of missing loci) were assessed with a 50% missing loci limit to check and remove samples above this value. To reduce the possible effect of linkage disequilibrium, SNPs were filtered to retain one SNP per locus for each sample. Thus, our filtering process retained > 17, 000 SNPs with a mean read depth of 47.27, and no samples presented high missingness values (> 50%) to be omitted from downstream analyses. Finally, due to the large number of SNPs recovered post-filtering, we retained a sub-sample of 10,000 using a down-sample approach for all downstream analyses to maintain computational efficiency. This method randomly samples the filtered SNP dataset so that all classes have the same frequency as the minority class.

Phylogenetic relationships of Southern Fontainea and Coastal Fontainea

Sequences of the filtered DArTseq SNPs were converted into nexus format for alignment and phylogenetic analysis. We used the R packages, decipher v.2.28.0 (Wright 2016) to generate FASTA alignments of the SNP dataset, followed by phangorn v.2.11.1 (Schliep 2011) and phytools v.1.9–16 (Revell 2012) to construct a maximum-likelihood tree with the best-fit GTR + I model and nearest-neighbour branch arrangements with 100 permutations to generate bootstrap support, representing the phylogenetic relationships for all 146 samples across both Fontainea species. For the final phylogeny, we pruned tips to the node which represented the position of divergence for each population.

Genetic structure and differentiation

To explore patterns of genetic structure we conducted a discriminant analysis of principal components (DAPC) using the R package adegenet v.2.1.10 (Jombart et al. 2020) to identify posterior probabilities of individual assignment to K-means groups. DAPC analyses were conducted twice to examine the influence of a priori groupings on the results. In the first analysis, sampling locations were used as a priori groups. In the second analysis, the find.clusters function was used to determine the number of groups (K) de novo, with optimal K selected as that with the lowest BIC value. For both analyses, the optimal number of principal components (PCs) to use in the DAPC was determined using the optim.a.score function.

For comparison to the DAPC, we also investigated admixture analysis using the R package LEA v.3.12.2 (Frichot and François 2015), which applies a sparse non-negative matrix-factorization (snmf) algorithm to generate ancestry coefficients. In this case, to estimate the best representation of K-means groups, we ran a cross-entropy analysis between K = 1–20 values, with 10 runs for each value of K. The values of K were then visualized in a cross-entropy plot with the lowest value of K prior to an increase used to indicate the closest estimate for the number of K-means groups.

We evaluated patterns of genetic differentiation between the populations using FST, which we calculated in the R package StAMPP v.1.6.3 (Pembleton et al. 2013) with 1,000 bootstraps to generate 95% confidence intervals and to evaluate whether the observed differentiation was statistically significant (p-value ≤ 0.05).

Environmental data

For the assessment of spatial genetic variation and to identify candidate SNPs for environmental adaptation we utilised eight explanatory variables: slope, aspect, hillshade (Hs), isothermality, mean temperature of the wettest quarter (TWQ), mean temperature of the driest quarter (TDQ), precipitation seasonality and precipitation of the warmest quarter (PWQ) downloaded from WorldClim database (Hijmans et al. 2005) at 30-arc second resolution (~ 1 km2 grids). For Hs we employed a multidirectional hillshade method that was a combination of slope, aspect, azimuth, and altitude to indicate the sun’s relative position and provide a representation of low-light areas across the study site. It is especially relevant to the study region in this research, where due to the unusual topographic relief, some areas only receive limited-no direct sunlight. Therefore, we might expect Fontainea populations that occur in low-light areas to be correlated to high hillshade in GEA analysis to indicate a driver of genetic differentiation at the local scale, population and/or species level.

Brunton et al (2023) identified this set of environmental factors were key predictors to the distribution of Southern Fontainea and Coastal Fontainea from a broader set of 23 variables using a stepwise procedure to exclude correlated factors and reduce overfitting for species distribution model predictions. In this way, the combination of landscape and climate factors may accurately represent the complex environmental characteristics across the study region.

Quantifying spatial genetic variation of Fontainea

To explore the influence of landscape heterogeneity on spatial genetic variation of Fontainea, we used a multiple matrix regression with randomization (MMRR) analysis on genetic, geographic, and environmental distance matrices. This method is similar to Mantel and partial Mantel tests, but is extended to incorporate multiple regression, which allows us to quantify how the dependent variable (i.e., genetic distance) changes in response to the independent variables (i.e., environmental distance, geographic distance; Wang 2013). Therefore, we can estimate and disentangle the relative influences of IBD and IBE on genetic variation (Wang 2013). Moreover, we can estimate the significance of individual variables as well as their combined effects on genetic variation (see Wang 2013 for details). We tested for IBD, IBE, and the effects of geographic distance against environmental distance by implementing MMRR on different matrices of pairwise genetic distance matrix of all Southern Fontainea and Coastal Fontainea and Euclidean distances for geographic and environmental factors between all sampling locations. A pairwise FST (genetic distance) matrix and geographic matrices were generated using the R package adegenet (Jombart et al. 2020) with the environmental matrix scaled and also generated in R using vegan v.2.6–4 (Oksanen et al. 2013).

MMRR was implemented in R using the function MMRR (Wang 2013) with 10,000 permutations to assess the additive effect of both independent factors (geographical distance and environmental). All matrices were scaled to a mean of zero and unit standard deviation using the scale function implemented in R before running the MMRR. With simulated data, this approach has produced robust results across a range of dispersal rates and proved accurate in situations of low-to-moderate gene flow (Wang 2013). This is especially relevant with the potential dispersal limitations suggested of extant Fontainea species (Lamont et al. 2016; Conroy et al. 2019).

Genome-environment associations and functional relationships of Fontainea

To detect genetic-environment relationships of the two Fontainea species, we tested the association between the set of climate and topographic variables described above, and the genetic variability across the sampled populations of each species. Although differentiation-based detections of statistical outliers have been frequently applied to identify signals of adaptive differentiation, GEA allow researchers to link signals of local adaptation to environmental variables. (De Mita et al. 2013; Forester et al. 2018).We identified putative signals of local adaptation using a combination of univariate and multivariate methods, following similar methods of Carvalho et al. (2021). We used latent factor mixed models (LFMM) and a partial redundancy analysis (pRDA) to represent the univariate and multivariate methods, respectively. Both methods have been shown to perform well in simulation studies and provide a compromise between power and error (Forester et al. 2018; Carvalho et al. 2021; Lotterhos 2023). All samples were included in GEA, that is, we analysed both Southern and Coastal Fontainea samples simultaneously to understand the relationships between the species and populations at candidate SNPs. We also ran these analyses excluding Coastal Fontainea samples to ensure that their inclusion did not influence our results (see Appendix S1 for results).

Latent factor mixed models analysis was implemented using the lfmm R package v.1.1 (Caye et al. 2019). We considered six uncorrelated environmental variables hypothesized to influence the distribution of Southern and Coastal Fontainea (Brunton et al. 2023). We scaled predictor variables and ran LFMM using the identified number of genetic clusters (K = 4) as latent factors to account for neutral genetic structure. Genotype information was represented by allele counts for each individual. We calculated and modified the genomic inflation factor (λ) to calibrate the p-value distribution (LFMM – Supplementary Fig S2; pRDA – Fig S3) according to the recommendations of François et al. (2016). Candidate SNPs under selection were identified after correcting for multiple tests using a False-Discovery Rate (FDR) of 10% (Benjamini and Hochberg 1995).

We performed pRDA in the R package vegan v.2.6–4 (Oksanen et al. 2013), using the same genotype matrix and environmental variables from the LFMM as input in pRDA. To control for population structure effects in the pRDA, we used the snmf assignments as a conditioning variable. Multicollinearity between predictors was determined using the variance inflation factor (VIF). We evaluated the overall pRDA significance and significance of each pRDA axis using permutation tests with 999 permutations and a p-value threshold of 0.05. Candidate SNPs were identified on all significant pRDA axes using a Mahalanobis distance approach (Capblancq et al., 2018) to compare pRDA and LFMM results (Carvalho et al. 2021). We adjusted p-values using λ and identified candidate loci using an FDR of 10% (Benjamini and Hochberg 1995). Selective gradients associated with each pRDA axis were calculated by correlating the environmental variables with each axis. Finally, we ran the pRDA excluding neutral SNPs using the candidate SNPs identified in either pRDA or LFMM analysis to visualize the relationships between the two species and populations of Southern Fontainea.

Additive polygenic score

We used additive polygenic scores (APS) to evaluate the cumulative effect of candidate SNPs identified from GEA and to assess if there were differences between Southern Fontainea and Coastal Fontainea. Additive polygenic scores were calculated as the total number of adaptive alleles within a particular environment. We identified the adaptive allele for each candidate SNP as the allele that exhibited a positive relationship with increasing values of a given environmental variable (Xuereb et al. 2018). Then we compared the average APS between F. australis and F. oraria to determine if there was a significant difference in APS between the two species using Wilcoxon Mann Whitney tests.

Results

Genetic relationships of Southern Fontainea and Coastal Fontainea

The maximum-likelihood tree we constructed was generally congruent with previously published assessments of these taxa (Rossetto et al. 2000b; Brunton et al. 2022), with (LW) weakly clustering separately from other Southern Fontainea populations and Coastal Fontainea. However, with our broader set of populations, LW was nested within a subclade of the Southern Fontainea populations; NB, CC, CRY, MTC and EHR (Fig. 2a) at the northern area of the Wollumbin caldera (Fig. 2b). Three Southern Fontainea populations; MNP, KRR and MRR formed a well-supported (100) subclade. Additionally, the population of NC, which is a notable geographic outlier for its position at a much higher elevation (> 750 m) compared to the remaining populations (~ 50–380 m), and WG were positioned on long branches, with moderate bootstrap support to indicate they are genetically distinct populations. Interestingly, HB, which is found in a gorge north of WG (Fig. 2b), was nested with Coastal Fontainea (FO).

Fig. 2
figure 2

Population differentiation in the two sister-species, Southern Fontainea and Coastal Fontainea from central, eastern Australia, including; a population-species relationships based on a GTR + I, maximum-likelihood tree generated with 10,000 reduced-representation markers with bootstrap values > 50 represented at branch nodes, with pie charts indicating K = 4 genetic clusters; b position of sampling locations in the Wollumbin caldera overlaid with current mean summer temperatures; c A priori discriminant analysis of principal components (DAPC) results from seven principal components identified using optimal a score and; d FST heatmap for comparing the defined Southern Fontainea populations and Coastal Fontainea. For figures a, c, and d FO represents Coastal Fontainea, with all other abbreviations representing Southern Fontainea populations

Full size image

Results obtained from the DAPC and SNMF algorithms presented conflicting results when for the number of genetic clusters for Southern Fontainea and Coastal Fontainea. The DAPC assignment without predefined groups (Supplementary Fig. 4a) the BIC value showed a similar value between 4 and 6 however, there was a slight increase in the value from four genetic clusters to support selecting K = 4. However, in contrast the snmf cross-entropy was less definitive, with separate plateaus at K = 6 and K = 8 clusters (Supplementary Fig. 4b). Due to the potential of the snmf algorithm to over-split, and therefore over-estimate the number of genetic clusters under a strong influence of IBD (Bradburd et al. 2018), we used a combined approach to explore outputs from the DAPC results and phylogenetic tree relationships to estimate a biologically appropriate K value for Southern Fontainea and Coastal Fontainea. The phylogram (Fig. 2a) also showed four distinct nodes of divergence among the various populations for both species to provided further support to using K = 4.

Investigating genetic structure and selecting an ‘optimal’ number of genetic clusters (K) remains an ongoing challenge that is frequently debated. More recently, a series of approaches is preferred to reach a consensus and justification to validate a suitable choice (Janes et al. 2017). As such, based on this collection of analyses, we determined K = 4 as the suitable representation of genetic clusters for the two Fontainea species. However, we acknowledge that K = 5 or possibly K = 6 could also be valid choices.

When sampling locations were used as a priori groups, DAPC results (Fig. 2c) supported a number of clear genetic groups that mostly reflected the geographic-genetic relationships indicated by the maximum-likelihood tree and FST values (Fig. 2d). Without predefined groups, DAPC analysis showed genetic assignments that also represented the FST results, which identified four distinct genetic clusters (based on optimal K value, see Supplementary Fig. 4a), to indicate that generally, genetic differentiation and structure of the two Fontainea species could be strongly associated with the geographic location of the sampled populations (Fig. 2b). A number of sampling locations were clustered together in the DAPC that generally reflected their geographical dynamics. For example, LW and NB formed a compact group on the plot in proximity to a dense cluster made up of the four other populations (CC, CRY, EHR, MTC) from the northern area of the study region. NC, HB, WG and FO formed a loosely grouped collection well differentiated from the northern populations. Notably, three populations (KR, MNP, MRR) were in a distinct group separated from all the other sampling locations.

Comparisons of pairwise genetic differentiation (FST) between all population pairs indicated significant differentiation among Southern Fontainea and Coastal Fontainea, highlighting a number of discrete populations (Fig. 2d). Although the difference in FST values were significant between all pairs of Southern Fontainea populations (p-values =  < 0.001), the effect size of differentiation between some populations was relatively small and typically reflected their geographical associations. For example, CRY had low genetic differentiation values compared to LW, NB, MTC, EHR and CC, which are all located around the northern area of the Wollumbin caldera. This pattern of geographic association and low relative magnitude of genetic differentiation was observed amongst several other Southern Fontainea populations, including between those in the Mooball area (MNP, MRR, and KR). Interestingly, perhaps the most geographically isolated Southern Fontainea population of Nightcap National Park showed strong levels of genetic differentiation only in comparison with the Mooball populations. Importantly, FST values of FO (Coastal Fontainea) showed relatively high genetic differentiation from all the Southern Fontainea populations (Fig. 2d).

A pairwise FST matrix based on the genetic clusters identified in the admixture analysis (K = 4) and phylogenetic relationships showed a significant difference (p-values =  < 0.001) in genetic distance between all four genetic groups (Supplementary Fig. 5).

Contributions of geography and ecology to the genomic divergence of Fontainea

Overall, the explanatory power of the MMRR was significant (R2 = 0.424, p < 0.001), in explaining the geographic and ecological variables that influence genetic divergence among Southern Fontainea and Coastal Fontainea. In particular, the MMRR showed that IBD (β = 0.424, p < 0.001) played a significant role in genetic differentiation between Southern Fontainea populations and among the Southern Fontainea and Coastal Fontainea (Fig. 3a). The regression coefficient for IBE was marginally associated (β = 0.330, p = 0.069) with genetic distance (Fig. 3b). There was no evidence of a statistically significant correlation between geographic distance and environmental distance (p =  > 0.05).

Fig. 3
figure 3

Scatterplots generated by multiple matrix regression with randomisation (MMRR) analysis showing relationships of isolation by distance a and b isolation by environment for Southern Fontainea (represented by blue points) and Coastal Fontainea (represented by red points)

Full size image

Evidence of divergent selection of Fontainea and functional associations

We identified 1,567 candidate SNPs using LFMM. The global pRDA model was significant (p-value < 0.001) six of the constrained pRDA axes were significant (Supplementary Table 2), with the first three axes found to explain 62% of the genetic variation in the SNP dataset (adjusted R2 = 0.11). Variance inflation factors ranged from 1.96–6.07 in pRDA, indicating that there is some collinearity between environmental predictors (Supplementary Table 3). GEA using pRDA identified 682 candidate SNPs (Fig. 4; Supplementary Table 4) that were correlated with six environmental variables. Mean temperature of the wettest quarter was associated with nearly a third (212) of candidates (Supplementary Table 2). Additionally, the topographic variables Hillshade (22% of total candidates) and Aspect (15% of total candidates) also showed significant association with candidate SNPs to indicate putative signals of local adaptation across Southern Fontainea and Coastal Fontainea populations. Combining LFMM and pRDA results, we identified 1,626 candidate SNPs.

Fig. 4
figure 4

Redundancy analysis (RDA) for the first two constrained ordination axes showing the candidate loci (colour-coded to environmental predictor) and neutral loci (grey circles) with six environmental variables (arrows) detected from 10,000 reduced representation SNPs for Southern Fontainea and Coastal Fontainea

Full size image

Additional pRDA analysis using only candidate SNPs explained 30% of the variance in the genetic dataset and revealed divergence between Coastal Fontainea and Southern Fontainea. The first constrained axis separated Coastal Fontainea individuals from Southern Fontainea populations and showed a relationship between Coastal Fontainea and slope and the temperature of the driest quarter (TDQ; Fig. 5a). Candidate only RDA split Southern Fontainea individuals into groups that generally correspond to their sampling location. The second constrained axis placed Southern Fontainea groups close together with the exception of individuals from the NC locality which show a relationship with temperature of the wettest quarter (TWQ.; Fig. 5a). The third constrained axis reveals a greater degree of separation among Southern Fontainea individuals. NB and KR individuals for example, show a relationship with TWQ while CRY individuals exhibit a relationship with Aspect and PWQ (Fig. 5b).

Fig. 5
figure 5

Redundancy analysis (RDA) for the first two constrained ordination axes showing: a relationship between individuals from Southern Fontainea and Coastal Fontainea populations (colour coded circles) and directionality of association with six correlated environmental variables (arrows) and; b for the first three constrained axes detected from 10,000 reduced representation SNPs

Full size image

Additive polygenic score

Estimates of APS did not identify a significant difference between the average number of adaptive alleles between Southern Fontainea and Coastal Fontainea (p-value = 0.46; Table S4).

Discussion

The capacity for local adaptation among small, threatened plant populations and the migration of advantageous alleles to track suitable habitats is critical to their long-term survival, especially in the context of rapidly changing and extreme climate projections (Leimu and Fischer 2008; Weeks et al. 2011; Kooyman et al. 2012). Our study indicated that geographic isolation and ecological differentiation are key drivers of genetic divergence of Southern Fontainea and Coastal Fontainea across the environmentally heterogeneous region of the Wollumbin caldera. Additionally, we found putative signals of local adaptation for the two Fontainea species related to a range of topo-climatic gradients that could be linked to important biological functions. These findings suggested that the number of putatively adaptive alleles varies between populations, which may indicate differences in the degree of adaptive potential of individual populations to certain factors such as climate gradients between inland and coastal populations. We also showed that Southern Fontainea and Coastal Fontainea together, could form part of a species continuum, that have become genetically differentiated as a result of geographic isolation and adaptation to local habitats. Based on genetic associations we identified four distinct groups that require separate conservation management to preserve the unique genetic resources and adaptive potential within each group. Consequently, to preserve these two threatened species into the future, a strategy that considers the adaptation potential of the distinct genetic groups, their local habitats and assisted migration is suggested.

Contributions of geography and environment to genomic divergence

We expected the diverse environments of the Wollumbin caldera to impose selective filtering across both Fontainea species. As such, we explored the influence of geographical features and environmental gradients on genomic divergence among Southern Fontainea and Coastal Fontainea. The observed genetic differentiation and the presence of an IBD pattern among some Southern Fontainea and Coastal Fontainea populations suggests that geographic isolation and limited gene flow may be contributing to a pattern of genetic divergence between these species. This is indicative of the numerous geographic barriers or isolated habitats typical of the study region, that may have restricted gene flow and promoted the genetic distinctiveness observed between some Southern Fontainea and Coastal Fontainea populations. In addition, the limited support of IBE suggests there may be some ecological divergence between species. Our genetic analysis also revealed patterns that align with the geographic positioning of both Fontainea species. While our study did not explicitly investigate the historical factors responsible for these patterns, such as the presence of refugia, these genetic findings suggest that the study region may have served as an area of refuge or genetic differentiation during past climatic fluctuations. Further research is needed to confirm and explore the nature of these patterns and their implications for the evolutionary history of these species. Microrefugia are a characteristic of complex, mountainous environments (González-Martínez et al. 2007; Holderegger and Thiel-Egenter 2009), which are a prominent feature of the Wollumbin caldera that contains numerous isolated valleys, gorges, steep slopes, cliffs and mountain ranges. This region is a noted climate refugium for a variety of species endemic to the area (Weber et al. 2014) including Southern Fontainea (Brunton et al. 2023), and the orographic features have had a major influence on the phylogeographic processes within the system (Kooyman et al. 2013). For example, some Southern Fontainea populations occur in largely continuous areas of complex, subtropical rainforest (LW, NB, EHR, CRY, CC), while others (MTC, KR) are found in regrowth which can be dominated by sclerophyll species such as Brushbox (Lophostemon confertus). Certain populations (MNP and MRR) survive in forest margins exposed to extreme edge effects with the latter found in a small, narrow strip of road reserve vegetation. Other populations are comprised of isolated individuals within regenerating rainforest patches, or mid-slope in deep gorges with creeks that often flood during significant rain events. One of the most isolated Southern Fontainea populations (NC) is located atop a steep, narrow escarpment at an elevation > 750 m in a once heavily logged forest. Finally, Coastal Fontainea occurs in littoral rainforest fragments exposed to salt-bearing winds. With these observations, we can begin to understand how local adaptation, geographic isolation and dispersal limitations could influence genetic heterogeneity along the range of the distribution for Southern Fontainea and Coastal Fontainea.

Phylogeny and genetic structure

Phylogenetic topology and structure analyses suggested a strong genetic structure related to geographical and ecological barriers in Southern Fontainea and Coastal Fontainea. There were limited signs to indicate gene flow between most of the Southern Fontainea populations and Coastal Fontainea to suggest some may have been reproductively isolated for some time. However, we identified a similar genetic structure among Coastal Fontainea and the Southern Fontainea populations of HB and WG, with only a weakly supported phylogenetic distinction from the additional Mooball populations (MNP, MRR, KR). Dispersal dynamics of Fontainea remain a critical knowledge gap however, we may assume our evidence of clear population differentiation, strong signal of geographic divergence and highly clumped distribution, indicates both species are dispersal limited. Generally, the population dynamics of other Fontainea taxa are similar to our findings, with geographically isolated populations displaying strong genetic differentiation and minimal evidence of gene flow among some adjacent populations (Conroy et al. 2019; Grant et al. 2019).

However, admixture was present in some geographically proximate Southern Fontainea such as between NB, CC, CRY, LW, MTC, and EHR to indicate genetic exchange has occurred among these neighbouring populations. This finding was somewhat unexpected given the NB, EHR and MTC populations are found at relatively higher elevations at the northern edge of the caldera, whilst CC and CRY are located south of the other three populations, in different sub-catchments at the foot of steep cliffs, which was an anticipated physical barrier to gene flow. Although directionality of gene flow has not been tested here, this implies there are two possible genetic interchange scenarios: 1) propagules have been transported by gravity or hydrochorous dispersal from the populations at higher elevations (NB, EHR, MTC), in a southerly trajectory down the steep escarpment to where the present LW, CC and CRY populations are found or, 2) a series of dispersal events in either direction (north or south) over the cliffs via animal dispersal. It is widely acknowledged multiple dispersal vectors are involved in seed transport (Hämäläinen et al. 2017). As such, it is possible a combination of both scenarios or an alternative version is involved in this case.

Three Southern Fontainea populations within a ~ 3 km radius of each other, MRR, MNP, and KR also had a level of population admixture and formed another genetically distinct group. Notably, these populations occur on the southern side of the Tweed valley and Tweed River, in a habitat corridor with limited physical barriers to gene flow. The remaining Southern Fontainea populations (NC, HB, WG) and Coastal Fontainea, are all generally geographically discrete and separated by prominent physical barriers (mountain ranges and gorges), which likely impedes gene flow between them. Even though Coastal Fontainea is well differentiated from most Southern Fontainea populations, phylogeny and structure analyses showed this species shared genotypes with certain Southern Fontainea populations (HB and WG). This finding was also highlighted in previous research (Rossetto et al. 2000a; Brunton et al. 2022).

Evidence of adaptive potential

We found numerous candidate loci that showed clear signatures of local adaptation that indicated ecological divergence resulting from geoclimate gradients across the species’ range may play a significant role in driving continued divergence among Southern Fontainea and Coastal Fontainea. Thermal tolerance is considered a major factor in shaping species diversity and determining plants species’ range (Reyer et al. 2013). Loci that are potentially under selection exhibited the strongest associations with temperature gradients in the case of Coastal Fontainea. This was not unexpected given where Coastal Fontainea is situated within 1 km of the ocean which experiences higher mean annual temperatures than the Southern Fontainea populations (Brunton et al. 2023). Southern Fontainea showed a strong relationship with both temperature and topographic factors that generally represent the distinct and heterogenous nature of the Wollumbin Caldera. For example, the most elevated population of NC showed a strong association with temperature at the wettest quarter (TWQ). TWQ, which corresponds to the hottest period in this region, implies that NC might possess adaptations to cope with the cooler temperatures anticipated at significantly higher elevations (~ 700 m asl) in contrast to all other Southern Fontainea and Coastal Fontainea populations. Despite NB and KR populations being located at lower elevations, they also showed a relationship with temperature during TWQ. It is possible that these populations experience unique climatic conditions or localised temperature patterns that influence the genetic traits of Southern Fontainea individuals in a way that correlates with TWQ, even at lower elevations. Further research and analysis would be needed to pinpoint the exact drivers behind this relationship. In the case of CRY individuals exhibiting a relationship with Aspect, this could suggest that the topographic aspect of their habitat, such as the direction in which slopes face or the angle of sunlight exposure, may play a role in shaping their genetic traits. The orientation of slopes and the resulting variations in sunlight and temperature can create microclimates that influence the adaptation of the population (Pincebourde et al. 2016; Bilir et al. 2021), demonstrating the importance of considering local environmental factors, including topography, when studying the genetic variation and adaptation of populations.

Consequences for Southern Fontainea and Coastal Fontainea taxonomy

Whilst Coastal Fontainea shows greater adaptation to the local conditions, it does not imply that it is a distinct lineage altogether. More plausible, is that Coastal Fontainea is an example of a population that has undergone differential environmental filtering through time, that is currently geographically separated and has undergone pronounced genetic drift, while potentially still being part of a broader species continuum that includes Southern Fontainea. Although our overall findings seem to challenge the current species boundaries between Southern Fontainea and Coastal Fontainea, a conclusive judgement would benefit from an analysis of the phenotypic variation associations, ideally in a common-garden experiment where the local conditions are unlikely to induce plastic responses. Consequently, we recommend genotype–phenotype sampling is necessary and should include representatives of the two species populations as well as an ecological outlier such as the Fontainea congener F. rostrata, which is found in drier rainforest habitats (Forster 1997). This would bridge taxonomic and geographical gaps in our dataset and could also provide greater information on the relative strength of selection and drift in the two southern Fontainea species, which is of importance given that the current inter-population and inter-species isolation and fragmentation is likely to remain at current levels or most likely be further exacerbated.

We also highlight that even with the several approaches of GEA analysis to indicate SNPs under putative selection and minimise false-positive rates it is possible the individual SNP loci we identified are putatively involved in local adaptation, and therefore, must be investigated further. Also, the evidence of local adaptation we uncovered does not necessarily support a role of ecological opportunity in initiating divergence between the two Fontainea species as adaptation may have occurred sometime after taxa became isolated from one another. Importantly, to assess the relative strength of divergent selection among the two species is complex, and any number of selective pressures may be operating independently, or acting together. This multidimensional aspect of selection was not specifically tested, and requires further research to understand the mechanisms involved between Southern Fontainea and Coastal Fontainea. Nevertheless, the overall pattern presented clear evidence of spatial climatic and topographic variation that appear to be important factors shaping local adaptation and possibly endemism between the two Fontainea species.

Implications for conservation

Threatened species found to have a narrow range of tolerance to projected shifting biomes, with limited ability to evolve in a rapidly changing environment may place species from southern rainforest lineages at significant risk of extinction (Crisp et al., 2009) and active conservation strategies could be necessary. Therefore, identifying genomic regions of local adaptation is essential to conserving genetic resources of Fontainea in the context of habitat loss and climate change. Several geographically and genetically distinct population clusters were found to support the genetic divergence between the northern locations (EHR, MTC, NB, CC, CRY, LW) of Southern Fontainea and the other Southern Fontainea and Coastal Fontainea populations. This enhances the importance of the northern margin of the Wollumbin caldera providing critical climate refugia for a host of ancient plant groups, which were spared the extreme conditions associated with paleoclimate fluctuations (Weber et al. 2014).

At the population and genetic group level, we uncovered evidence supporting the existence of four distinct genetic clusters, which included differentiation between certain Southern Fontainea populations and Coastal Fontainea. Interestingly, Coastal Fontainea was found to share some genetic structure with several Southern Fontainea populations. More detailed ecological and phenotypic studies are needed to determine whether these taxa should each be considered their own species, or as varieties or ecotypes of a single species; other genetic studies of Fontainea suggest that either scenario may be appropriate based on observed values of genetic variation (Rossetto et al. 2000a) and phylogenetic topology (Brunton et al. 2022). The recent discovery of a proposed new Fontainea species, F. sp. Coffs Harbour A.S. Benwell (PlantNET), also suggests a taxonomic revision of the genus could be timely to resolve uncertainties presented from our genetic data.

Irrespective of taxonomic classification, conservation and management efforts could be enhanced by conserving the genetic and adaptive traits of the four genetic groups highlighted in this study. Irrespective of taxonomic classification, in light of our pRDA analysis, prioritising the conservation of genetic variants identified through LFMM and pRDA, particularly the candidate SNPs (> 1600), can significantly inform and enhance conservation and management strategies. These SNPs reflect substantial genetic differentiation among Coastal Fontainea and Southern Fontainea populations and their associations with specific environmental factors, offering valuable insights for targeted conservation efforts. Both Southern Fontainea and Coastal Fontainea face an elevated risk of population decline due to genetic factors (low genetic diversity, strong genetic differentiation) as well as ongoing threats from continued habitat decline and fragmentation, under-representation of sites within protected areas, land use modification and climate change (Rossetto et al. 2000a; Brunton et al. 2023). These threats have led to the establishment of ex-situ and in-situ Coastal Fontainea populations for conservation purposes (Brown et al. 2016). There is an opportunity to adopt similar management objectives that include adaptive potential outlined in this study for Southern Fontainea. Equipped with the information from this research and previous studies (Rossetto et al. 2000a, b; Brunton et al. 2022), we advocate for immediate action to preserve Southern Fontainea populations, also using a similar strategy proposed for F. rostrata or, Deep Creek Fontainea (Conroy et al. 2019) to ensure the long-term survival of both species. Specifically, for longstanding conservation of both Fontainea species, assisted migration of populations to suitable habitats may be necessary but should consider the genetic composition, adaptation potential and local environment conditions of each population. Considering an adaptation potential perspective will provide environmental managers an improved understanding of the broader ecology of the two species, which can be used to review and identify significant conservation units. This approach also aids in developing conservation plans that consider the capacity to survive with the rapid onset of changing environments, and consequently enhance the conservation of Southern Fontainea and Coastal Fontainea. Collectively, our study has provided critical insights into the combined contributions of geo-climate and local adaptation on the diversification of these two narrow-range, endemic species, which could also have consequences for the conservation of other rare, threatened rainforest species.

Conclusion

While it is well-recognised that local adaptation is crucial for species’ survival in a complex, changing environment and plays a central role in diversification (Kawecki and Ebert 2004; Hereford 2009), there are few studies investigating the evolutionary dynamics of narrow-range, rainforest endemics. Closely related species such as Southern Fontainea and Coastal Fontainea provide ideal models to explore this subject, owing to their geographically discrete distribution across a fine, but complex landscape, while their threatened status adds to the significance of this study. We presented molecular evidence that the spatial isolation of Southern Fontainea and Coastal Fontainea has led to local adaptation by some populations to their distinct environments. Survival of rainforest species through paleoclimate oscillations relied upon tracking suitable conditions, however this necessitates adequate dispersal capacity and accessibility to appropriate habitats. All Fontainea species produce a fleshy-fruit suitable for animal-mediated dispersal, however the key mode of seed transport and genetic movement is unknown for these cryptic species. Crucially, we showed that Southern Fontainea and Coastal Fontainea may be part of a larger, albeit spatially segregated species continuum, but suggest that distinct genetic clusters should be managed as separate units due to putative signals of local adaptation within some populations. Small, disjunct populations are expected to continue to decline from increased genetic differentiation as a result of genetic drift, bottlenecks and limited adaptive potential to climate projections (Lesica and Allendorf 1992). Threatened species such as the two focal study taxa face an elevated risk of extirpation from a variety of stochastic events including wild-fires, droughts or demographic fluctuations that can lead to a build-up of deleterious alleles. Since Southern Fontainea and Coastal Fontainea are confined to narrow, isolated areas of suitable habitat, mostly lacking formal protection, they Qare at a high-risk of decline caused by these stochastic processes, as such immediate conservation action is required to ensure their long-term survival.