Genomic, Habitat, and Leaf Shape Analyses Reveal a Possible Cryptic Species and Vulnerability to Climate Change in a Threatened Daisy

Olearia pannosa is a plant species listed as vulnerable in Australia. Two subspecies are currently recognised (O. pannosa subsp. pannosa (silver daisy) and O. pannosa subsp. cardiophylla (velvet daisy)), which have overlapping ranges but distinct leaf shape. Remnant populations face threats from habitat fragmentation and climate change. We analysed range-wide genomic data and leaf shape variation to assess population diversity and divergence and to inform conservation management strategies. We detected three distinct genetic groupings and a likely cryptic species. Samples identified as O. pannosa subsp. cardiophylla from the Flinders Ranges in South Australia were genetically distinct from all other samples and likely form a separate, range-restricted species. Remaining samples formed two genetic clusters, which aligned with leaf shape differences but not fully with current subspecies classifications. Levels of genetic diversity and inbreeding differed between the three genetic groups, suggesting each requires a separate management strategy. Additionally, we tested for associations between genetic and environmental variation and carried out habitat suitability modelling for O. pannosa subsp. pannosa populations. We found mean annual maximum temperature explained a significant proportion of genomic variance. Habitat suitability modelling identified mean summer maximum temperature, precipitation seasonality and mean annual rainfall as constraints on the distribution of O. pannosa subsp. pannosa, highlighting increasing aridity as a threat for populations located near suitability thresholds. Our results suggest maximum temperature is an important agent of selection on O. pannosa subsp. pannosa and should be considered in conservation strategies. We recommend taxonomic revision of O. pannosa and provide conservation management recommendations.


Introduction
Populations of threatened plant species are typically small and fragmented, suffering low genetic diversity associated with increased genetic drift and elevated inbreeding rates [1][2][3][4][5]. Additional to these factors, a changing climate places increasing pressure on threatened populations. Plants can respond in one of three ways in the face of climate change: adapt to the new environment, migrate to locations with suitable conditions, or perish [6,7]. In this era of accelerated climate change, plants may not have the capacity to respond through in-situ adaptation, particularly those already suffering loss of genetic associates with environment and used habitat suitability modelling to reveal environmental limitations to the species' distribution. We discuss the outcomes from these analyses in terms of their implications for the conservation management.

Study Species and Sampling
Olearia pannosa is a low shrub (~1.5 m) which spreads by producing root suckers [34]. The species can be found in a range of floristic communities and climatic conditions within southeastern Australia. There are currently two subspecies with similar ecology but can be differentiated by leaf shape: O. pannosa subsp. pannosa (silver daisy bush), and O. pannosa subsp. cardiophylla (velvet daisy bush). O. pannosa subsp. pannosa is listed as vulnerable under the Australian Government's Environment Protection and Biodiversity Conservation Act 1999 [36]; and vulnerable under the South Australian National Parks and Wildlife Act 1972 [37]. O. pannosa subsp. cardiophylla is also of conservation concern and has been listed as vulnerable under the Victorian Advisory List of Rare or Threatened Plants in Victoria-2014 [38].

Study Species and Sampling
Leaf samples were collected from a total of 191 individuals from 25 localities spanning the ranges of both subspecies ( Figure 1, Table A1). Each sample (individual plant) consisted of five leaves: one leaf collected from each aspect/side of the plant (North, South, East, West), and one from the central apex point of the plant-generally the highest/near highest growing point of the plant. Nearest neighbour individuals were avoided to minimise relatedness between genotyped individuals and obvious clones and root suckers were avoided. Samples were dried in a sealed container with silica beads. Collections were made during winter and spring of 2018. Voucher specimens were also collected at each locality (except the Victorian stands), and where there was evidence of morphological variation between individuals at a site, multiple specimens were collected. Voucher specimens consisted of a 20-30 cm small branch taken from a growing point of the plant, which was labelled and pressed at the earliest convenience ( Figure A1). Vouchers were submitted to the State Herbarium of South Australia, and an identification was provided for each collection locality (Table A1).

DNA Extraction, Sequencing, and Filtering
All DNA extraction, library preparation and sequencing were conducted at the Australian Genome Research Facility (Adelaide, Australia). The Machery-Nagel Nucleospin Plant II Kit was used to extract DNA prior to undergoing double digest restrictionassociated DNA sequencing (ddRADseq) [39]. The restriction enzymes Pstl and MseI were used for the digestion, then barcoded adaptors were ligated to the restriction site overhangs for each sample. Samples were then pooled, and the DNA size selected using Blue Pippin (Sage Science, Beverly, MA, USA) (Wide) to retain fragments 280-342 bp. The resulting library was amplified via Polymerase Chain Reaction (PCR) for 11 cycles using indexed primers. Libraries were assessed by gel electrophoresis (D1000 ScreenTape Assay, Agilent, CA, USA), quantified by qPCR (KAPA Library Quantification Kits for Illumina, Roche, Basel, Switzerland), and then single read sequenced (150 bp) on the NextSeq 500 system using NextSeq 500 high Output Kit v2 reagents (150 cycles).
A total of 637,395,304 sequences were generated across 192 O. pannosa samples and we processed these using the STACKS pipeline [40,41] at AGRF (Melbourne, Australia). Reads were de-multiplexed using the unique sample adapter barcodes and filtered for both read quality and the presence of a restriction site. The resulting FASTQ files were then trimmed to the size of the shortest read, minus 2 bp to account for differences in read length due to any variation in barcode length. Reads were then stacked according to similarity, forming tags of reads. Tags which appear across all samples were collated (catalogue tags), and genotypes were then allocated to the common polymorphic sites. The collated SNPs across all individuals were then filtered using the following settings: minimum number of reads required at a stack to call a homozygous genotype = 5; minor allele frequency, below which a stack is called a homozygote = 0.05; minimum minor allele frequency to call a heterozygote = 0.1 (above 0.05 but below 0.1, a stack is called 'unknown'). Computed genotypes were then output as a single VCF file (Supplementary File S1) for further filtering and analysis.
We used VCFtools [42] to filter SNPs with a minor allele frequency of less than 0.1 and with a depth of coverage of less than10 reads per individual. We removed SNPs with >25% missing data across all individuals and individuals with >25% missing SNP calls. We then used PLINK [43] to identify SNPs in linkage disequilibrium with the 'indep-pairwise' command, with one SNP in a pair of linked SNPs (r 2 > 0.8) being removed from the dataset.
Levels of ploidy were assessed by examining the distribution of allele dosage at heterozygous sites. For this, the counts of reads supporting each allele at a heterozygous site were extracted from the VCF file using a custom Perl script, and each individual's distribution was plotted. True allelic proportions for diploids, triploids and tetraploids at heterozygous sites are 1:1 (0.5), 1:2 or 2:1 (0.33 or 0.66), and 1:3, 2:2, or 3:1 (0.25, 0.5, or 0.75), respectively. No sample deviated from a normal distribution centred on a frequency of 0.5, suggesting that all samples are diploid.
For each SNP across all samples, inbreeding coefficient (F) values were estimated using the G-stats calculator in GENODIVE v. 3.2 [44,45]. As negative F indicates greater than expected heterozygosity under Hardy-Weinberg equilibrium, and could be indicative of paralogous reads being stacked together in the STACKS pipeline [46], we used VCFtools to remove SNPs with significantly negative F (significance assessed using permutation tests with 10,000 permutations and a significance threshold of p < 0.05).
After the filtering steps, 175 individuals and 9997 SNPs remained. We used this dataset to perform the PCA, DAPC, and ADMIXTURE. After filtering for highly related individuals, 169 individuals and 9997 SNPs remained in the dataset. We used this dataset to perform F ST , genetic diversity and inbreeding statistics, and AMOVA analyses.

Genetic Structure Analysis
We then analysed the population genetic structure across all samples in the filtered SNP. To identify the most likely number of genetic clusters (K), we used two different genetic clustering analyses. First, we used the non-model-based Discriminant Analysis of Principal Components (DAPC) of the R package adegenet [47]: genetic data were transformed into uncorrelated components using principal component analysis (PCA). The number of genetic clusters was then defined using k-means, a clustering algorithm that looks for the value of K that maximises the variation between groups. The Bayesian Information Criterion (BIC) was calculated for K = 1-20, and the K value with the lowest BIC suggesting the optimal number of clusters. A discriminant analysis was then performed on the first 80 principal components with the function 'dapc', implemented in R, in order to efficiently describe the genetic clusters and assign samples to each cluster. We also used ADMIXTURE [48], a model based genetic clustering algorithm, to estimate the most likely number of clusters (K) in our dataset. ADMIXTURE was also run for K values 1-10 and the most likely value of K was assessed by comparing the cross-validation errors (cv errors) between runs, with the lowest cv error indicating greatest support. We repeated each ADMIXTURE analysis ten times and used CLUMPAK to combine the results and construct bar plots of individual assignment to clusters for the most supported values of K using DISTRUCT [49,50]. We then calculated pairwise F ST between all collection localities using the R package StamPP [51]. This was then exported to Splitstree to create a neighbour-joining (NJ) tree [52].

Kinship, Genetic Diversity, and Inbreeding Analysis
The dataset was filtered for all first-degree relations (parent-siblings and siblingsiblings) following methods outlined in Dutheil [53] to avoid these influencing our statistical results. The king-cutoff option in PLINK 2 [54] was used to remove one in a pair of samples with a kinship coefficient greater than or equal to 0.177 (6 individuals were removed, Tables A2 and A4). A kinship coefficient of 0.25 corresponds to first-degree relations, so the 0.177 threshold should capture all of these. This step should also have removed any likely clones from the dataset, which was an important consideration given O. pannosa is known to reproduce clonally. We then used PLINK 2 [54] to calculate pairwise kinship between individuals. We then calculated average kinship within each locality in R v.4.0.4 As the genetic structure analysis identified such strong genetic differentiation between three of the clusters, we separated the dataset according to these groupings. Each dataset was refiltered for minor allele counts of 1 with VCFtools [39] to remove non-variants for subsequent analyses (See Table A2 for number of variants remaining at each collection locality).
For the collection localities in each of the three genetic groups, we separately estimated the inbreeding coefficient (F), and observed (H O ) and expected (H E ) heterozygosity in GENODIVE v3.2 [44,45]. For the collection locality COUL (COULTA), there were insufficient samples to calculate these statistics and so this collection locality was not included in this analysis. A one-way ANOVA and Tukey pairwise comparisons were used to assess group differences in kinship, F, H O and H E using R v.4.0.4 [55]. To assess the distribution of genetic diversity amongst and within the genetic groups, we conducted a nested Analysis of Molecular Variance (AMOVA) using GENODIVE v3.2 to explore genetic variation within individuals, among collection localities nested within groups, and among groups. Further AMOVA analyses were then run separately for O. pannosa subsp. pannosa and O. pannosa subsp. cardiophylla (there was not enough data to complete the analysis for Olearia (FR)) to compare how genetic variation was distributed within individuals and among collection localities between these two genetic groups.

Leaf Trait Analysis
In order to investigate the environmental influence on leaf shape and the variation that exists between collection localities, leaf traits were identified and analysed. In total 110 individuals were measured from across 22 collection localities (note that Victorian collection localities were excluded from this component of the study as we did not have suitable material for the measurements). Leaf traits were based on the methods of Lowe and Abbott [56] and measured using ImageJ software v 1.53e (https://imagej.nih.gov/ ij/download.html, accessed on 10 November 2020) [57] from photo images of voucher samples (See Figure A1 for examples). Three measurements were taken; leaf length (L) (measured from leaf base to apex), leaf width (W) (measured at the widest point of leaf), and the distance along the midrib to the widest point of the leaf from the base of the leaf, giving width location (WL). Leaf area measurements were square root transformed to satisfy the requirement of normality of residuals. Leaf length (L), width (W) and the distance at which the width measurement was taken from the base (WL) were used to estimate leaf area (LA), employing the formula We estimated leaf ovality or roundness by dividing the length along leaf mid rib to widest point of the leaf by total leaf length. A value around 0.5 indicates an oval leaf, <0.5 indicates a more heart-shaped leaf, and >0.5 a 'gingko' shaped leaf. The ratio of leaf length to leaf width (calculated as leaf length divided by width), was also determined. Where possible, five leaves were measured from each voucher specimen, with five individuals sampled from collection locality. A standardised leaf sampling protocol was conducted by sampling the first mature leaf from the base of the specimen, with subsequent leaves sampled working towards the apex. In some circumstances direct leaf measurement was not possible due to folding or damage, in which case width was calculated by measuring distance from leaf edge (at widest point) to leaf midrib and multiplying by two.
A principal component analysis (PCA) was performed using the 'prcomp' function in R v.4.0.4 [55] on the measurements for all measured leaf samples to assess differences in leaf shape among the three genetic clusters. Trait values were centred and scaled. Oneway ANOVA and Tukey pairwise comparisons were used to assess group differences in principal component axes 1 (PC1) and 2 (PC2) using R v.4.0.4 [55].

Redundancy Analyses-Both Genetic Data and Leaf Variation
To assess the degree to which the environment explains observed genetic and phenotypic variation, partial redundancy analyses (RDA) were performed for O. pannosa subsp. pannosa samples. We focused on this genetic group as they were distributed over a wide environmental gradient and by omitting O. pannosa subsp cardiophylla samples confounding genetic structure can be removed from the analysis. A set of 44 environmental (e.g., climate, soil, and landscape) variables were sourced from stacks of 9-s resolution (i.e., approximately 250 m pixels) rasters [58][59][60]. Environmental values were extracted for each raster at the coordinates of collection localities of O. pannosa subsp. pannosa using the RASTER analytical package in R [61]. A subset of these extracted environmental variables with low collinearity were assessed using the 'cor' function in R v.4.0.4 [55]. They were: mean annual maximum temperature ( • C), mean annual minimum temperature ( • C), annual precipitation (mm), total nitrogen (%), total phosphorus (%), organic carbon (%), and elevation (m).
Redundancy analyses (RDA) were performed on (a) the allele frequencies of all SNPs in collection localities of O. pannosa subsp. pannosa and (b) leaf trait measurements for all samples of O. pannosa subsp. pannosa in order to assess how much of the genetic and phenotypic variation among our collection localities can be explained by spatial (i.e., resulting from isolation by distance) and environmental (i.e., indicating selection) factors (Supplementary File S2) using the VEGAN analytical package v2.5-6 in R [62]. First, spatial coordinates for (a) each sampling stand or (b) each sample were centred and converted into third degree polynomials. An RDA of allele frequencies-all spatial polynomials was then performed using the 'rda' function. The 'ordistep' function was then used to carry out forward stepwise model building for the RDA, adding in one polynomial to the model at a time and assessing significance using permutation tests. This revealed x, y, and x 2 to be significant contributors to the model and as such only these variables were retained for the RDA. The variance between the selected spatial and all environmental variables was assessed using the 'varpart' function and significance of the partitioning was assessed using an analysis of variance (ANOVA)-like permutation test for RDA with the 'anova.cca' function. In each case, we tested whether any of the genetic or phenotypic variance was significantly explained by environmental variables by running an RDA on the models (a) 'allele frequencies-environmental variables conditioned on spatial variables' and (b) 'leaf traits-environmental variables conditioned on spatial variables'. To assess which of the environmental variables provided the most explanatory power in the models, we performed the same stepwise model building described above for the spatial variables on the environmental variables.

Habitat Suitability Modelling
For the largest genetic group, O. pannosa subsp. pannosa, MaxEnt and Random Forest were used to assess which environmental variables best explain habitat suitability given the species' current distribution. Occurrence data were sourced from collection locality sampling records, consisting of 105 locations. We modelled current distribution using the BCCVL online tool (http://www.bccvl.org.au/, accessed on 25 February 2021); Ref. [63] by statistically comparing confirmed presence sites against predictor spatial layers. Predictors included national soil data provided by the Australian Collaborative Land Evaluation Program ACLEP (www.clw.csiro.au/aclep, accessed on 25 February 2021), and a current Australian climate baseline of 1976 to 2005 provided as bioclim layers by BCCVL [64].
Distributions were predicted based on current climate and soil using MaxEnt and Random Forest, which are widely used machine-learning algorithms for mapping habitat suitability [65,66]. Models were run using 10,000 pseudo-absences confined to the land area within 250 km of a convex polygon surrounding presence locations. Starting with the full set of potential variables included in the edaphic and climatic layer sets, we sequentially removed both highly correlated (i.e., r > 0.8) and poorly performing predictors contributing to 1% or less of the model.

Genetic Structure Analysis
Strong genetic structure was identified across the range of O. pannosa and was supported by all genetic analyses (Figures 2, 3, A2 and A3, Table A3). While the BIC criterion from the DAPC analysis ( Figure A3) and CV error from the ADMIXTURE ( Figure A2) both indicate that K = 5-6 is the most likely value of K, there is strong evidence of three key genetic groups within O. pannosa which have genetic barriers that cannot be explained by geographic isolation (Figure 2). This was shown in samples from the Flinders Ranges in semi-arid South Australia identified as O. pannosa subsp. cardiophylla. These were genetically distinct from all other samples, even when collected from the same geographic location (see collection site Melrose (MEL) in Figures 2, 3, A1 and A2, Table A3). The second such genetic break is found on the Fleurieu Peninsula of South Australia. Despite evidence of admixture in these collection localities (Figure 2), the geographically proximate Newland Head (NEW) and Victor Harbor (VIC) (~12 km apart) fall into different clusters in the PCA and different genetic groups in the DAPC analysis (Figures 2, 3 and A3, Tables A1 and A3).
We therefore present three key genetic groups, acknowledging that there is further substructure within them. These groups are:

1.
The southern Flinders Ranges samples identified by the State Herbarium of South Australia as O. pannosa subsp. cardiophylla (Table A1, Peter Lang, pers. comm.). From here on, we will refer to this group as Olearia (FR).

2.
Samples from Coulta (COUL) on the Eyre Peninsula to Victor Harbor (VIC) on the Fleurieu Peninsula (excluding the Olearia (FR) specimens collected at Dutchman's Stern (DUT) and Melrose (MEL)) were mainly identified as O. pannosa subsp. pannosa (State Herbarium of South Australia, P. Lang 2021, personal communication; Table S1). There is evidence of some substructure within this group, with geographically isolated sites Coulta (COUL) and Cummins (CM) on the Eyre Peninsula forming a distinct cluster. Notably, there was some admixture between genetic groups, particularly in the samples collected at Victor Harbor (VIC), which were identified as O. pannosa subsp. cardiophylla (Table A1). From here on this genetic group will be referred to as O. pannosa subsp. pannosa.
All other samples ranging from the NEW (Newland Head) collection locality on the Fleurieu Peninsula across to the Victorian collection localities were identified as O. pannosa subsp. cardiophylla (State Herbarium of South Australia, P. Lang 2021, personal communication; Table A1). However, admixture was also detected in several collection localities geographically close to stands of O. pannosa subsp. pannosa (e.g., Newland Head (NEW)). From hereon we will call this genetic group O. pannosa subsp. cardiophylla.   Table A3 for full output and Figure A1 for images of full voucher specimens).
Across the whole dataset, our analysis of pairwise genetic differentiation (F ST ) showed a considerable range of differentiation between collection localities (F ST = 0.00-0.75; Figure 3, Table A3). Within the Olearia

Genetic Diversity, Kinship, and Inbreeding Analysis
Across the dataset, the majority of genetic variation was found within individuals, but that a significant proportion of variation was also distributed between the three genetic groups (AMOVA; 42.7% and 29.9% respectively; Table 1), however there were differences between subspecies. Within O. pannosa subsp. pannosa, most genetic variation was found within individuals (67.8%) and only 11.2% between collection localities (Table 1). However, within O. pannosa subsp. cardiophylla a high proportion of the genetic variation was distributed between collection localities (43.2%; Table 1).

Leaf Trait Analysis
Principal component analysis of all leaf measurements revealed clear clustering of samples, which broadly aligned with the three genetic groups, but with some overlap, indicating morphological intergrades ( Figure 5). Principal component 1 (PC1) explained 68.5% of the variance and mainly explained by variation in length, width, and surface area of leaves ( Figure A5)

Environmental Associations
The redundancy analysis of O. pannosa subsp. pannosa collection localities revealed that spatial distribution explained a significant proportion of the variance in allele frequencies (Table 2, Figure A4). The partial redundancy analysis, which looked at the proportion of variance explained by the seven environmental variables (mean annual maximum temperature ( • C), mean annual minimum temperature ( • C), annual precipitation (mm), total nitrogen (%), total phosphorus (%), organic carbon (%), and elevation (m)) after conditioning on spatial distribution, revealed that 11% of the variance in allele frequencies is significantly explained by the included environmental factors (Table 2, Figure A4). Forward selection of the environmental variables retained only mean maximum temperature as a significant contributor to the model, suggesting that this environmental variable associates with genetic variation.
The redundancy analysis of leaf shape within O. pannosa subsp. pannosa samples included the same seven environmental variables as for the allele frequency RDA. When running the RDA with only spatial distribution as the explanatory variable, 13% of the variance in leaf shape was significantly explained Table 2). With environment as the explanatory variable conditioned on spatial distribution only 3% of the variance is explained and was non-significant (Table 2). When partitioning the variance in leaf measures among the spatial and environmental variables, 17% of the variation is explained by the combination of both environment and spatial distribution. The environmental variables (particularly temperature and precipitation) do co-vary with spatial distribution, so it is difficult to tease apart the effects of the environment on leaf shape (i.e., environmentally determined) from differences in leaf shape that may be due to demography. However, as for the genotype RDA, the forward selection process revealed maximum annual temperature, as well as total nitrogen, as significant contributors to the model, once again suggesting that temperature associates with variation among these populations.

Habitat Suitability Modelling
Final models were fitted with two edaphic (% soil clay-clay30; % available water capacity-pawc1m), two temperature predictors (mean minimum temperature of the coldest month-bio6; mean maximum temperature of the warmest month-bio5), and three rainfall (mean annual rainfall-bio12; mean rainfall of the driest month-bio14; precipitation seasonality/coefficient of variation-bio15) predictors ( Figures A5 and A6). Distribution models showed O. pannosa subsp. pannosa occurs within a consistent climatic band intermediate between the mesic areas and the more arid inland. Peak suitability was associated with mean annual rainfall of 400-500 mm (occupied niche 330-670 mm) and mean summer temperature maxima <29 • C (MaxEnt; <30 • C for Random Forest; occupied niche 25-31 • C). For MaxEnt, precipitation seasonality had the highest variable relative contribution followed by mean maximum temperature of the warmest month. For random forest, mean maximum temperature of the warmest month had the highest variable importance followed by mean annual rainfall.

Discussion
Our analysis of genomic and leaf shape variation in the threatened daisy bush, Olearia pannosa, unearthed a complex picture that is not in agreement with the current subspecific classifications for the species. A key finding of our study was the identification of a likely cryptic species Olearia (FR) in the Flinders Ranges which exists only in small, fragmented collection localities. For the two subspecies, O. pannosa subsp. pannosa and O. pannosa subsp. cardiophylla, we identified a strong east/west genetic divide with signatures of admixture between the two subspecies and substructure within them. Genetic diversity was reduced in Olearia (FR) compared to O. pannosa subsp. pannosa, and was substantially lower in the small, isolated localities of O. pannosa subsp. cardiophylla. Our analyses of associations between environmental and genetic and phenotypic variation identified a strong association with maximum temperature. This finding was further supported by habitat suitability modelling, which identified maximum temperature and precipitation as key variables in determining species' distribution. Our results have implications for the conservation management for O. pannosa, particularly for the newly identified genetic group Olearia (FR) in the Flinders Ranges.

Genetic Structure, Diversity, Inbreeding, and Kinship
We found evidence of strong differentiation among collection localities, with clear support for three major genetic groups across all analyses and substructure with them. Previously, O. pannosa plants found in the southern Flinders Ranges were thought to be the same subspecies as O. pannosa subsp. cardiophylla found in southeastern South Australia and Victoria. However, our analysis reveals that these southern Flinders Ranges plants (Olearia (FR)), sampled at Dutchman's Stern (DUT) and Melrose (MEL), are not only strikingly distinct genetically from their neighbouring collection localities of O. pannosa subsp. pannosa but also from all other collection localities of O. pannosa subsp. cardiophylla in the southeastern South Australia and Victoria. This is in clear contrast to the current classification based on morphology, where O. pannosa in the Flinders Ranges is considered to be part of the same subspecies as O. pannosa subsp. cardiophylla found in Victoria. The apparent lack of gene flow among Olearia (FR) and neighbouring O. pannosa subsp. pannosa collection localities strongly suggests that there are reproductive barriers present between these two clusters and these collection localities likely represent two distinct species. Notably, both Olearia (FR) collection localities have high genetic differentiation between them, low genetic diversity, and high historical inbreeding. These findings are consistent with results from other recent genomic studies of threatened and endangered species [67][68][69]. As far as we know, this is a new finding, signalling the need for taxonomic description of this newly discovered genetic group and a conservation management strategy to be developed to ensure its survival.
The remaining samples could be assigned to either the O. pannosa subsp. pannosa or O. pannosa subsp. cardiophylla genetic groups, which are split by an east to west genetic divide. These two genetic groups differed significantly in levels of genetic diversity, inbreeding and kinship, with O. pannosa subsp. cardiophylla having lower levels of genetic diversity and inbreeding but higher levels of kinship. Interestingly, a self-compatibility test was undertaken at the Brisbane Ranges locality which found that individuals could self-pollinate and set seed, and inbreeding was flagged as a concern in O. pannosa subsp. cardiophylla [34]. Certainly O. pannosa subsp. cardiophylla shows signs of population genetic decline indicated by the low genetic diversity and high kinship, and our findings consistent with the fact that the Victorian populations are now small and demographically isolated. There is also potential for elevated levels of clonality in these collection localities [70,71], but as our sampling protocol was designed to minimise the likely number of clones in our dataset, it is not possible for us to calculate the proportion of clones within collection localities. We did, however, detect one likely clone in the Kangaroo Island (KIN) collection locality, which also had the lowest levels of genetic diversity. Nevertheless, due to the limitations in our study, we advise further research to test mating system dynamics.
We observed genetic admixture on the southern Fleurieu Peninsula and southeastern SA, with a shift from a collection locality with high assignment to the O. pannosa subsp. pannosa group at Victor Harbor (VIC) to a collection locality with high assignment to the O. pannosa subsp. cardiophylla group at the closely located (~12 km) collection locality at Newland Head (NEW). Interestingly, Newland Head (NEW) had the variation in pairwise kinship between individuals within O. pannosa subsp. cardiophylla. Genetic admixture patterns were reflected in the leaf shape analysis, with intermediate morphologies described at these sites. Distribution records for both subspecies do not align with these findings, likely due to the difficulties with identification in these admixed sites, highlighting the intrinsic difficulty of defining subspecies boundaries. There was further substructure within O. pannosa subsp. cardiophylla, which is represented by several isolated and geographically distant collection localities. Our findings in Victoria are in agreement with the results of Smith, James and Ladiges [34] who also found that the collection site Rushworth (RUS) was genetically differentiated from the Brisbane Ranges (BR) and Anglesea (ANG) sites which could be explained by the large geographic distance between these localities (~150 km and 220 km respectively).
Overall, the genetic structure results reflect the life history traits and disjunct range of O. pannosa as gene flow will be limited, resulting in greater genetic drift. Therefore, the strong differentiation in isolated collection localities is unsurprising. Similar patterns have been found in other threatened species [67,68]. For example, genomic analysis of the endangered Australian daisy (Rutidosis leptorrhynchoides) found strong genetic differentiation between geographically isolated collection localities [69].

Leaf Shape and Taxonomy
Our leaf shape analysis showed that the three main genetic clusters display relatively distinct leaf morphologies. Specimens measured from the O. pannosa subsp. pannosa genetic cluster had smaller, more oval shaped leaves compared to plants belonging to the Olearia (FR) and O. pannosa subsp. cardiophylla genetic groups, which we found to overlap more in size than shape. Previous taxonomic work on the species has classified subspecies by morphological traits and demonstrated clear differences between the leaf shapes of O. pannosa subsp. pannosa in South Australia and O. pannosa subsp. cardiophylla in Victoria [34], similar to what we find here. However, this previous analysis showed that Olearia (FR) samples (previously O. pannosa subsp. cardiophylla in South Australia) clustered with O. pannosa subsp. cardiophylla samples from Anglesea (ANG), Victoria. Here, we show that whilst the leaf shape of Olearia (FR) is more similar to O. pannosa subsp. cardiophylla than to O. pannosa subsp. pannosa, there are differences between them, with Olearia (FR) having more heart-shaped leaves. Geographic distance between Olearia (FR) and O. pannosa subsp. cardiophylla is extremely large (~315 km between Melrose (MEL) and the closest O. pannosa subsp. cardiophylla locality at Newland Head (NEW)). We also show that genetic divergence is high between Olearia (FR) and all other O. pannosa collection localities, suggesting they may be a distinct species, although reproductive isolation is yet to be tested.
We observed evidence for admixture among O. pannosa subsp. pannosa and O. pannosa subsp. cardiophylla on the Fleurieu Peninsula and in parts of the Southeast of SA, indicating that gene flow is possible between these groups. If Olearia (FR) and O. pannosa subsp. pannosa were not reproductively isolated, then the strong signal of genetic structure we observe between such geographically close collection localities inhabiting very similar environments is difficult to explain. The advance in sequencing technologies in recent years has led to numerous studies that have unearthed previously unrecognised 'cryptic species' [71][72][73][74][75][76][77][78][79]. These are genetically distinct taxa that are classified as the same species due to morphological similarity [80]. This species complex has undergone several taxonomic revisions since its first description by Hooker in 1851 (see Smith, James and Ladiges [34] and references therein), reflecting the large amount of phenotypic variation within the species. These results are also consistent with the ongoing difficulties with taxonomy and phylogeny in the Olearia genus [29][30][31][32][33]. Our findings strongly suggest a cryptic species complex in O. pannosa. For effective threatened species management plans, a firm understanding of likely species or subspecies boundaries is essential, and our results highlight the importance of combining both genetic and morphological results methods for this approach.

Environment Likely Shapes Genotype, Phenotype, and Distribution
We conducted further analyses to identify the effect of the environment across the genotype, phenotype, and distribution of O. pannosa subsp. pannosa. These analyses identified maximum temperature as a key environmental variable that is a likely selective pressure on this group. Our analyses of genetic and leaf shape variation both identified maximum temperature as explaining a proportion of the measured variation. These findings are consistent with a number of studies that have found genetic signals [81][82][83][84][85] and morphological variation [86][87][88][89][90] associated with temperature being an agent of selection on plant species. Furthermore, the current distribution of O. pannosa subsp. pannosa appears to be closely related to climate. Our climate suitability models indicated specific constraints with respect to mean annual rainfall and summer maximum temperature gradients, with the species occurring within a transition zone between the more mesic and arid zones of the surrounding region [91]. Many collection localities of O. pannosa subsp. pannosa sit at the warm/arid thresholds of these key climatic gradients beyond which modelled probability of occurrence drops off markedly.
Taken together these results suggest that the current climatic trends towards higher summer maxima and increasing aridity are likely to have a detrimental impact on many collection localities. However, the few localities in cooler, wetter habitats may, in fact, benefit in the short term from climate change, leading to shifts in distribution and abundance. Combined, these results highlight the potential vulnerability of O. pannosa subsp. pannosa to the imminent effects of climate change and should be considered when planning conservation strategies for the species.

Olearia pannosa (FR)
Olearia pannosa (FR) had low levels of genetic diversity across both collection localities and the collection locality MEL (Melrose) had the highest level of inbreeding across the whole dataset. We also found moderately high genetic differentiation between the two collection localities Dutchman's Stern (DUT) and Melrose (MEL). Concerningly, as there are only a few known stands of these individuals, all with small census estimates, this suggests an endangered species, or at least a subspecies which contains distinct genetic variation that will be lost if these they are not conserved. An immediate taxonomic review is essential to address these issues. We suggest a multi-level management strategy centred for this group. Firstly, we recommend in-situ threat abatement (i.e., protection from grazing pressure) and natural recruitment (i.e., creating disturbance to encourage recovery from the seedbank) [92]. We also recommend developing a strategy to aid the recovery of genetic diversity through a translocation plan which facilitates gene flow between the two isolated stands [92]. This could be developed through the ex-situ establishment of a managed seed production orchard to provide high quality seed [93]. Furthermore, further research could be informative, such as whole genome sequencing and breeding experiments, to ascertain (a) whether Olearia (FR) truly is a distinct, isolated species from O. pannosa and (b) what the functional genomic differences are between these divergent taxa.

Olearia pannosa subsp. pannosa
Within O. pannosa subsp. pannosa, there is evidence of some degree of historical gene flow across collection localities and low population substructure. This genetic group had higher levels of genetic diversity compared to O. pannosa subsp. cardiophylla. Observed heterozygosity was consistently lower than expected heterozygosity and was reflected in the significantly higher levels of historical inbreeding in O. pannosa subsp. pannosa than the other two genetic groups. Inbreeding has previously been flagged as a concern for O. pannosa [34] and historical selfing in this group may be one explanation for the results.
Due to the substructure that we found between isolated collection localities on the Eyre Peninsula (Coulta (COUL) Cummins (CM) and Cleeves (CLE)) and the rest of the group, we recommend that management efforts be focused here. Facilitating gene flow between these stands through translocations or increasing connectivity between may be advantageous.
As our results identify maximum temperature as potential agent of selection for this group, stands in more arid regions may experience detrimental effects from climate change. We recommend close monitoring, particularly of the more northern and arid localities. Apart from collection localities on the Eyre Peninsula, there is evidence of gene flow and higher levels of genetic diversity compared to other genetic groups. For these stands, we recommend a management approach which maintains connectivity and gene flow among localities. However, the potential future effect of aridification on this genetic group should be considered in future management strategies, and a climate-adjusted seed sourcing strategy [18] could be implemented (i.e., introducing seed from more arid areas to southern stands to bolster their adaptive capacity to arid conditions).

Olearia pannosa subsp. cardiophylla
We found that O. pannosa subsp. cardiophylla had the lowest levels of both historical inbreeding and genetic diversity, yet the highest levels of kinship. We also found that three of the four sites which had lowest genetic diversity and inbreeding in this group were the only collection sites that had to have highly related individuals removed across the whole dataset. There was considerable substructure in this genetic group, particularly between sites in Southeast SA and Victoria, and between sites within Victoria. It is likely that the large range disjunction of this group is driving genetic drift and genetic divergence. To develop an appropriate conservation management plan for O. pannosa subsp. cardiophylla, there are several knowledge gaps which still need to be addressed, particularly around the genetic divergence we observed between O. pannosa subsp. cardiophylla and the two other genetic groups, and the consequences that the high levels of substructure we found within this genetic group may have on mixing seed between localities (e.g., outbreeding depression). We suggest the establishment of a common garden trial which explores the effects on fitness and genetic diversity of genetic mixing between localities within this genetic group and amongst the three genetic groups. Until this information is known, in a scenario with high structure but low genetic diversity and inbreeding, Ottewell, Bickerton, Byrne, and Lowe [88] suggests a cautious management approach to minimise the risk of outbreeding depression with a focus on in-situ recovery of genetic diversity through the seed bank and a composite translocation strategy [19].

Conclusions
In summary, combining range-wide genomic and leaf shape analysis of the threatened Australian daisy Olearia pannosa revealed a likely cryptic species (Olearia (FR)) and identified the need to redefine the geographic intergrade between existing subspecific boundaries (O. pannosa subsp. cardiophylla and O. pannosa subsp. pannosa). Levels of genetic diversity and inbreeding differed between the three genetic groups, suggesting each requires a separate management strategy which we define below.

•
Olearia (FR)-We suggest a management strategy centred around in-situ threat abatement and disturbance to encourage seedbank recovery, the recovery of genetic diversity through a translocation plan to encourage gene flow between the two isolated stands, and the development of a seed production orchard to supply the seed. • O. pannosa subsp. pannosa-The main priority is to facilitate gene flow between existing stands or to increase connectivity between them, especially in an arid-tomesic direction since it appears that maximum temperature is an important agent of selection.

•
O. pannosa subsp. cardiophylla-It is likely that the large range disjunction of this group is driving genetic divergence. To develop an appropriate conservation management plan, several knowledge gaps still need to be addressed, particularly to assess the potential for outbreeding depression. Until then, we recommend priority is given to in-situ recovery of genetic diversity.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/life11060553/s1, Supplementary File S1-The variant call format (VCF) file will be made available as supplementary (Olearia_pannosa.vcf). Reads and mapping files will be archived at the NCBI SRA (accession number TBA). Supplementary File S2-the covariates from the redundancy analysis (genetic_RDA_covariates.xlsx). Data Availability Statement: Reads and mapping files will be archived at the NCBI SRA (accession number TBA).

Acknowledgments:
The authors wish to thank Amelia Hurren, Rohan Cleeves, Randall Bates, and Vicki-Jo Rusell from Trees for Life (SA), Rani Hunt from the Department of Environment, Land, Water and Planning (Victoria) and Chris Lindorff from Trust for Nature for providing their assistance and expertise throughout this and for assisting us with sample collection. We'd like to thank Peter Lang and Helen Vonow from the Sate Herbarium of South Australia for providing their expertise in identification of vouchers specimens and advise on taxonomy. Thanks to Dan Duval at the South Australian Seed Conservation Centre for sharing advice and expertise. Thanks also go to Sam Turner-Blyth and Kiri Marker for their assistance with drafting and editing the manuscript. Finally, a special thanks to Jeremy Austin for his assistance behind the scenes.

Conflicts of Interest:
The authors declare no conflict of interest.           Table A6. Results of the one-way analysis of variance (ANOVA) conducted to detect differences in leaf shape between the three genetic groups O. pannosa subsp. pannosa, Olearia (FR), and O. pannosa subsp. cardiophylla. We have compared the differences between PC1 (width, length, and leaf surface area), and PC2 (ovality and location of widest point of the leaf).