Genomic diversity guides conservation strategies among rare terrestrial orchid species when taxonomy remains uncertain

Abstract

Background and Aims Species are often used as the unit for conservation, but may not be suitable for species complexes where taxa are difficult to distinguish. Under such circumstances, it may be more appropriate to consider species groups or populations as evolutionarily significant units (ESUs). A population genomic approach was employed to investigate the diversity within and among closely related species to create a more robust, lineage-specific conservation strategy for a nationally endangered terrestrial orchid and its relatives from south-eastern Australia.

Methods Four putative species were sampled from a total of 16 populations in the Victorian Volcanic Plain (VVP) bioregion and one population of a sub-alpine outgroup in south-eastern Australia. Morphological measurements were taken in situ along with leaf material for genotyping by sequencing (GBS) and microsatellite analyses.

Key Results Species could not be differentiated using morphological measurements. Microsatellite and GBS markers confirmed the outgroup as distinct, but only GBS markers provided resolution of population genetic structure. The nationally endangered Diuris basaltica was indistinguishable from two related species (D. chryseopsis and D. behrii), while the state-protected D. gregaria showed genomic differentiation.

Conclusions Genomic diversity identified among the four Diuris species suggests that conservation of this taxonomically complex group will be best served by considering them as one ESU rather than separately aligned with species as currently recognized. This approach will maximize evolutionary potential among all species during increased isolation and environmental change. The methods used here can be applied generally to conserve evolutionary processes for groups where taxonomic uncertainty hinders the use of species as conservation units.

Diuris, genotyping by sequencing (GBS), microsatellite markers, next-generation sequencing, Orchidaceae, population genomics

INTRODUCTION

Conservation of species relies on evidence-based, accurate and meaningful data collected with consistency and scientific rigour (Tear et al., 2005; Frankham et al., 2012). Conservation strategies can be improved greatly with knowledge of genetic variation identified through population genetic methodologies (e.g. Xiao et al., 2006; Martins et al., 2015; Zhang et al., 2015; Ahrens and James, 2016; Turchetto et al., 2016). Indeed, population genetics has recently been used as the basis for the development of a conservation prioritization framework for rare plants (Ottewell et al., 2016). The utility of understanding the genetic variation within a species can be crucial for prioritizing certain populations over others (Ahrens and James, 2015; Dillon et al., 2015). Particularly with the advent of newer genomic technologies over the past 5 years, we are able to clarify and resolve patterns of population structure that were previously unattainable. Reduced-representation molecular techniques are particularly useful in that they allow low-cost genome typing of hundreds or thousands of individuals [e.g. genotyping by sequencing (GBS), Elshire et al., 2011; double digest restriction site-associated DNA sequencing (ddRADseq) Peterson et al., 2012]. Having thousands of markers gives researchers and practitioners the power to uncover fine-scale structure of a conservation unit regardless of current taxonomic classification.

Defining species boundaries can be a difficult task (Fazekas et al., 2009; Hausdorf, 2011). Taxonomic uncertainty can lead to negative outcomes for conservation strategies (Frankham et al., 2012). Indeed, erroneous classifications can lead to potential financial, legal, biological and conservation repercussions (Hey et al., 2003). Biological reality suggests that a species can be difficult to confirm when morphology is influenced by phenotypic plasticity, hybridization or intraspecific ploidy variability, among other characteristics (Fazekas et al., 2009; Fitzpatrick et al., 2015). If maintaining evolutionary potential is a part of the final conservation management goal, categorization that is not constrained by species boundaries and includes genetic diversity metrics may lead to improved conservation outcomes. In such cases, next-generation sequencing techniques have the potential to differentiate conservation units by identifying genetic lineages that are part of dynamic evolutionary processes (i.e. contemporary gene flow, mating dynamics and environmental resilience; Allendorf et al., 2010; Funket al., 2012; Hoffmann et al., 2015). The use of evolutionarily significant units (ESUs) provides an alternative categorization of biological entities (e.g. ecotypes, subspecies or populations) and may be suitable as conservation units for managing closely related species. ESUs represent distinct lineages and take into account genetic interactions among species that influence gene pools but are not constrained by species delimitations (Moritz, 1994). Classifying populations that have substantial reproductive isolation is important for maximizing their potential for adapting to future environmental change (Funk et al., 2012). The use of ESUs instead of species classifications has the power to improve conservation strategies by incorporating genetic diversity in addition to biological and morphological information. As currently constituted, existing international conservation policy tends to overlook genetic diversity as an important component of the species (Laikre, 2010).

Consider the family Orchidaceae, which exhibits highly diverse floral morphology and is under pervasive threat from global change (Swarts and Dixon, 2009). This floral diversity is the result of unprecedented genetic diversity that originated early in orchid evolution (Mondragon-Palomino and Theissen, 2009). The breadth of floral variation within genera may create the illusion of greater species diversity than is actually present, and has led to a history of taxonomic exaggeration (Pillon and Chase, 2007). It has been found that taxonomically difficult groups tend to have more endangered species within them (Pilgrim et al., 2004). Further, species delimitations in the Orchidaceae remain tenuous in some groups (Flanagan et al., 2006; Devey et al., 2007; Pillon and Chase, 2007; Bateman et al., 2010; Hedrén et al., 2012), suggesting that using an approach other than species boundaries may sometimes be warranted to guide conservation. Under circumstances where orchid species are hard to define, understanding patterns of genetic diversity has been successful for developing conservation strategies (Swarts et al., 2009; Chung et al., 2012; Chen et al., 2013; Qian et al., 2014; Yang et al., 2014; Pandey et al., 2015; Pedersen et al., 2016; Ross and Travers, 2016). In many cases, genetic differentiation (F_ST) within rare orchid species remains very low (0·092) compared with more common species (0·279), suggesting that there is generally high gene flow and little population structure in rare orchids (Phillips et al., 2012) but no studies were found that recorded F_ST among closely related orchid species. Despite the power of reduced-representation next-generation sequencing protocols to resolve fine-scale population structure, there have been few, if any, studies that use these new techniques to guide the conservation of orchids with ‘fuzzy’ taxonomic margins. By using these new techniques, we can instead focus on protecting evolutionary processes in orchid populations, which shifts the conservation question from ‘which species should be conserved?’ to ‘which populations are evolutionarily important?’

We used five putative orchid species from south-eastern Australia to test the reliability of using current species delineations, made with morphological characteristics, as the units of conservation. We hypothesized that morphologically diagnosed species boundaries in the study group within Diuris (Orchidaceae) do not reflect the evolutionary lineages identified by population genomics. Our goal was to look beyond taxonomic classifications and investigate genomic diversity to create a robust conservation strategy. We used morphometrics, GBS and microsatellite markers to test our hypothesis by exploring genetic divergence and interspecific population structure at fine scales. The genetic results were compared with morphological characteristics to guide the identification of appropriate conservation targets, and we suggest new strategies for conservation and recovery of this Diuris species complex.

MATERIALS AND METHODS

Study species

Four co-occurring, closely related terrestrial orchid species on the Victorian Volcanic Plain (VVP; Fig. 1) were used to examine the relationship between current species delineations, morphological differentiation and genetic divergence. The VVP (Fig. 1) is a basalt-derived bioregion in south-eastern Australia that contains <1 % of its original vegetation due to fragmentation and urban development (Williams et al., 2015). Diuris basaltica is a nationally endangered terrestrial orchid that occurs in only three known locations on the urban outskirts of the greater Melbourne district. A fourth historical population is now extinct. At the time of sampling, two locations contained only 2–5 individuals. The third site supposedly had hundreds of individuals, but none flowered during the 2 years of the study. However, an ex situ population, maintained at the Royal Botanic Gardens Victoria glasshouses, contained several hundred individuals cultivated for conservation from plants sourced originally in the late 1990s and early 2000s from both the third population and the fourth now-extinct population. Of the remaining species, D. gregaria is considered endangered in the state of Victoria and is mostly confined there, with some occurrences in South Australia (Fig. 1). It exhibits a clumping habit not found in the other species. Diuris chryseopsis and D. behrii have large distributions that extend beyond Victoria (Fig. 1). These two species exhibit morphological floral and leaf characters that tend to be larger in size than the preceding two species. Neither one is considered to be threatened or endangered in the state of Victoria or nationally in Australia. This Diuris group mimics nectar-bearing flowers of other plants in the area and lures the same generalist pollinators via simple food deception (Backhouse and Lester, 2010). Species descriptions focus only on morphological characters and geographic locations as diagnostic features (Jones, 1998, 2006). Diuris monticola was used as an outgroup because it has similar floral morphology to the other four species and occurs outside of the VVP in an Alpine region, but within the distribution of D. chryseopsis and D. behrii (Fig. 1).

Fig. 1.

Diuris populations of four species sampled in the Victorian Volcanic Plain (VVP) and an outgroup sampled from the alpine region. Historical distributions for the four VVP species were generated from the Atlas of Living Australia (ala.org.au).

Sample collection took place in Spring of 2012 and 2013. A total of 298 individuals were collected from 17 populations across Victoria (Table 1; Fig. 1). Six other populations, recorded within the past 10 years, were visited but no flowering individuals were present during the two field seasons. At least 20 samples were collected from each population, except where populations were too small [4–12 individuals (Table 1)]. One population (PAR) exhibited high morphological variation so sampling was increased to cover the variation observed. At least one representative voucher specimen (3–4 in most cases) was taken from each population for post-hoc identification at the Royal Botanic Gardens Victoria by botanical experts and lodged at the National Herbarium of Victoria. While we used post-hoc identifications as putative population classifications, the presence of more than one ‘species’ was possible. However, if two or more species were present within a population, we would expect morphological and/or genomic data to distinguish them. Of the 17 populations sampled, two were identified as D. basaltica, six as D. gregaria, five as D. chryseopsis, three as D. behrii and one as D. monticola (Table 1). Only individuals whose flower morphology could be measured were sampled. Morphology was measured and recorded in situ, with the exception of the glasshouse- (GH) grown D. basaltica population. Only one D. basaltica population (eight plants) was measured in situ. The remaining 24 D. basaltica plants included in the analysis were measured in the ex situ GH population. Leaf material was collected and stored for DNA extraction.

Table 1.

Population identity, assigned species category and estimates of population diversity estimates from microsatellite markers are given for all populations and number of samples, reads and reads per sample for the genotyping by sequencing (GBS) analysis

Population	Putative species	Microsatellite markers					GBS
		N	Ne	Ho	He	FIS	N+reps	N reads	Reads/Ind
GH	D. basaltica	24	2·86	0·45	0·53	0·14	23	5 770 048	250 871·7
SK	D. gregaria	4	2·11	0·41	0·47	0·12	4	1 372 623	343 155·8
ML	D. gregaria	12	2·92	0·44	0·53	0·17	12	3 543 504	295 292·0
CG	D. gregaria	11	2·37	0·30	0·44	0·31	11	2 442 792	222 072·0
Woo	D. gregaria	10	2·32	0·31	0·44	0·29	9	2 898 598	322 066·4
Par	D. behrii	36	3·83	0·37	0·57	0·34	44	14 427 416	327 895·8
YP	D. gregaria	8	2·39	0·31	0·51	0·39	5	1 054 672	210 934·4
CW	D. chryseopsis	20	2·53	0·42	0·53	0·22	19	4 530 488	238 446·7
IV	D. chryseopsis	20	2·59	0·34	0·53	0·36	25	7 087 302	283 492·1
DK	D. chryseopsis	20	2·01	0·41	0·38	−0·08	27	8 244 430	305 349·3
WG	D. basaltica	7	1·85	0·34	0·33	−0·03	6	4 303 107	717 184·5
RC	D. chryseopsis	20	2·77	0·44	0·50	0·12	27	16 276 802	602 844·5
PF	D. chryseopsis	20	2·41	0·32	0·49	0·35	28	14 751 125	526 825·9
CR	D. behrii	20	2·46	0·39	0·56	0·30	15	7 661 890	510 792·7
CM	D. behrii	24	3·17	0·32	0·53	0·40	24	18 737 356	780 723·2
AC	D. gregaria	20	2·64	0·37	0·42	0·13	23	16 184 928	703 692·5
ALP	D. monticola	25					15	11 246 681	749 778·7
Overall		298	2·28	0·37	0·48	0·24	317	140 533 762	443 324·2

N, number of individuals; N_e, number of effective alleles; H_o, observed heterozygosity; H_e, expected heterozygosity; F_is, inbreeding coefficient, N_+reps, number of individuals + the number of replicates that made it through the GBS pipeline; N reads, number of reads for all individuals; Reads/Ind, average number of reads per individual.

Morphology

For each individual sampled, we measured 11 vegetative and floral characteristics considered diagnostic for the species and compared them with species descriptions. As diagnostic features of this Diuris group were based on floral morphology, eight of the measurements were taken from the flower. These included the number of flowers per flowering stalk, width (W) and length (L) of the dorsal sepal, W and L of the lateral petal, W and L of the labellum, and flower diameter. The non-floral measurements included the W and L of the leaves and height of the plant. Length and width ratios (L/W) for all traits where W and L were measured were used for morphological analyses. Using the ratio of some traits allowed us to account for petal shape (e.g. short and wide or long and thin) and variation in size due to growing environment, as some individuals were measured while growing in the glasshouse. Size for other traits (height and flower width) informed us if size was a deterministic character. A principal co-ordinates analysis (PCoA) was performed on the morphological data set using the Vegan package (Oksanen et al., 2013) and the ‘dist’ function in R v3.1.0 (R Core Team, 2014) to calculate the distance matrix for the PCoA with the ‘dudi.pco’ function. In addition, we analysed the 11 measurements as a principal components analysis (PCA) using the ‘prcomp’ function in R and a PCoA using the ‘dudi.pco’ function in R.

Microsatellite markers

Genomic DNA was isolated for both genetic protocols from 20 mg of silica-dried leaf tissue from each individual using the DNeasy Plant Mini Kit (QIAGEN, Hilden, Germany) by the Australian Genome Research Facility (AGRF; Adelaide, South Australia). The 13 microsatellite markers used in this study were developed previously for this Diuris group (Ahrens and James, 2014). Loci were amplified in three multiplex PCRs using a three-primer approach (Blacket et al., 2012). PCRs were carried out under the following conditions: 5 min denaturation at 95 °C, followed by 32 cycles of 95 °C for 30 s, annealing at 60 °C for 1 min 30 s and 72 °C for 30s, with a final extension at 60 °C for 30 min. Fragment analysis was performed via capillary electrophoresis on an ABI 3130×l (Applied Biosystems, Foster City, CA, USA) by Macrogen (Seoul, South Korea). Alleles were visualized and scored in Geneious Software version 7.0 (Biomatters Ltd, Auckland, New Zealand).

Allele size data were imported into the Polysat (Clark and Jasieniuk, 2011) package in R for PCoA on Bruvo’s genetic distance. This distance metric was developed to measure genetic distance similar to dominant data while taking into account mutational distances between alleles. Observed heterozygosity (H_o), expected heterozygosity (H_e) and effective number of alleles (N_e) were calculated on allele size data in Genodive (Meirmans and Van Tienderen, 2012). Inbreeding (F_IS) was calculated in Genodive using Weir and Cockerham’s (1984) methodology for f.

Genotyping by sequencing

The GBS library was prepared following a method adapted from Elshire et al. (2011) using 300 ng of genomic DNA per individual and the HpaII restriction enzyme. We barcoded, amplified and normalized individual libraries before pooling and sequencing, which reduced the range in per sample coverage (Nicotra et al., 2016). Two pooled GBS libraries of 154 and 178 multiplexed samples, respectively, were sequenced on an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) using 100 bp paired-end reads. Technical replicates (n = 34) of the first lane were added to the second lane.

AXE (v0.2.6, https://github.com/kdmurray91/axe) was used to demultiplex the reads. For each demultiplexed sample, a custom script was used to remove adaptor contamination due to fragment sizes <100 bp by aligning each read to its pair and removing read through if present. Overlapping reads were merged with FLASh v1.2.11 (Magoc and Salzberg, 2011), with maximum overlap = 70. Reads were filtered by quality, with 5′ trimming turned off, and trimmed to 64 bp with sickle se (v1.33, https://github.com/najoshi/sickle). For compatibility, pseudo barcodes were added to each read and those reads were imported into TASSEL v3.0.172 (Glaubitz et al., 2014). Candidate single nucleotide polymorphisms (SNPs) were identified in TASSEL using the UNEAK pipeline with default settings, including a minimum coverage per tag of 5, error tolerance of 0·03, maximum frequency of 0·5 and a minimum frequency of 0·05 for minor alleles. The output genotype matrix was further filtered in R using custom script (https://github.com/borevitzlab/GBSFilteR/blob/master/Diuris_hmn_clean ing Fig.R) first by removing failed samples with <4000 SNPs called (273 unique samples remained) and secondly by removing SNP loci with <200 (‘stringent’) or 70 (‘relaxed’) samples called. Potential paralogous loci with >45 samples being called as ‘het’ were also removed. Technical replicates were used to confirm reproducibility. Only one sample of each replicate was retained for downstream population analyses. The ‘relaxed’ data set was used whenever a genetic distance matrix was calculated, as genetic distances can be calculated with good results for genomic data with up to 90 % missing data (Fu, 2014). The ‘stringent’ data set was used for analyses with complex modelling components (i.e. STRUCTURE and SNAPP) where large amounts of missing data can give spurious results or make them difficult to interpret.

Both the ‘relaxed’ and ‘stringent’ data sets were used for genetic distance calculations. The ‘dist’ function in R was used to calculate the genetic distance (Euclidean) of both data sets with and without the outgroup, D. monticola. PCoAs were created for the ‘stringent’ and ‘relaxed’ data sets with adegenet v1.4-2 (Jombart, 2008). An unrooted dendrogram was created using all populations. A second unrooted dendrogram was built using paired F_ST values between populations calculated in adegenet. We used the coalescent model in SNAPP, an add-on to BEAST v2.3 (Bryant et al., 2012; Bouckaert et al., 2014), to investigate genetic divergence. The data set used for SNAPP included six individuals selected randomly per population from the ‘stringent’ data set to increase computational speed. For the priors, we calculated the mutation rates based on the data set as described in the SNAPP manual (Bouckaert and Bryant, 2012). We ran the MCMC (Markov chain Monte Carlo) algorithm for 300 000 generations, recording every 1000 tree and log of the output for visualization and convergence, respectively. Tree topology was visualized in DensiTree. Branch lengths are relative because we only use polymorphic SNPs (Bouckaert and Bryant, 2012; Bryant et al., 2012). The ‘stringent’ data set was used in STRUCTURE v2.3, a Bayesian clustering MCMC method (Pritchard et al., 2000) to determine the number of genetic clusters (K) in the data. A burn-in length of 100 000 iterations and a run length of 150 000 iterations were used in an admixture model with correlated allele frequencies among populations for each K-value between 1 and 10. Each K-value was run 12 times, and ten runs with the highest lnP(D) scores were used for later consensus. The optimal K-value was determined using the Evanno ΔK method (Evanno et al., 2005) in STRUCTURE Harvester (Earl and vonHoldt, 2012). A consensus of the optimum K-value was created in CLUMPP (Jakobsson and Rosenberg, 2007), and the graphical parameters were drawn in the program DISTRUCT (Rosenberg, 2004).

To understand the relationship between morphological and genetic differentiation, the natural log of the genetic distance matrix computed from the ‘relaxed’ data set was plotted in R against the natural log of the morphological pairwise distance matrix. A generalized linear model was fit to test the relationship between morphological distance and genetic distance.

Results

Morphology

During field collections, morphological variation was apparent within populations (Supplementary Data Fig. S1). Most characters measured in situ were in agreement with those given in published descriptions, except for D. basaltica, which were consistently different (Fig. S1). We controlled for size differences between glasshouse-grown and in situ plants as described in the Materials and Methods. The PCoA showed that there was no clear distinction between VVP species, but D. gregaria and D. behrii did show a small degree of separation (Fig. 2A). The D. monticola samples grouped together but showed no clear separation from other species. The x-axis explained 94·95 % of the variation, and the y-axis explained 4·83 % of the variation. The PCA and PCoA with all 11 measurements show the same pattern as the PCoA with ratios (Supplementary Data Figs S2 and S3). The ellipses in the PCA indicated that all species had overlapping morphological characteristics, but D. basaltica, D. gregaria and D. monticola tended to be slightly smaller.

Fig. 2.

Principal co-ordinates analyses for (A) morphology, (B) microsatellites, (C) genotyping by sequencing (GBS) ‘stringent’ data set and (D) GBS ‘relaxed’ data set. Colours represent putative species, and the shapes within colour groups represent different populations. Inset in (D) PCoA of the ‘relaxed’ data set with all 17 populations (including D. monticola).

Microsatellite analysis

Genotypes based on 13 microsatellite markers were obtained for 270 samples. None of the markers amplified in the D. monticola samples. Some individuals had more than two alleles per locus in three of the 13 loci, which could be due to locus duplication or polyploidy. Putative polyploidy was locus specific and not population or species specific. Global F_ST was estimated at 0·03 within the VVP group (four species on the VVP), indicating high gene flow and little genetic isolation. F_IS fell between –0·08 and 0·40, with an overall F_IS of 0·24, which indicated that inbreeding was moderate, and two populations had a slightly negative coefficient. The PCoA did not show significant species or population differentiation (Fig. 2B). Three species (D. basaltica, D. chryseopsis and D. behrii) form one overlapping cluster although the D. gregaria individuals showed off-centre aggregation. The two-axis PCoA accounted for 7·27 % of the variation on axis 1 and 6·40 % on axis 2.

Genotyping by sequencing

In total, 141 million reads were distributed amongst 273 unique samples, 25 samples removed prior to analysis and 34 replicated samples for 332 total samples (Table 1). We identified 118 764 SNPs across the 332 samples which is within range for GBS studies (Poland and Rife, 2012; Grabowski et al., 2014). Two replicates were removed due to low read count. In a dendrogram of the genetic distance matrix generated with the ‘relaxed’ data set (Supplementary Data Fig. S4), all 32 replicates were closest to their replicate counterpart. We removed each replicate counterpart for all of the following analyses. From the 273 unique genotypes, we isolated 9793 SNPs for the ‘relaxed’ data set and 874 SNPs for the ‘stringent’ data set.

The PCoA with the ‘relaxed’ data set showed similar patterns to the PCoA calculated from the ‘stringent’ data set (Fig. 2C, D). The two-dimensional PCoA from the ‘stringent’ data set explained 6·44 % of the variation (Fig. 2C), whereas the PCoA from the ‘relaxed’ data set explained 7·54 % of the variation (Fig. 2D). The PCoA from the ‘relaxed’ data set showed greater separation between groups, especially between the D. gregaria and the main cluster of D. chryseopsis, D. behrii and D. basaltica. In both GBS PCoAs, one D. behrii population (DK) separated from the main cluster. This was due to low levels of population structure with a small amount of genetic differentiation (y-axes explain 2·46 and 3·38 %, respectively).

The unrooted dendrogram showed the genetic relationship among individuals (Fig. 3). Individuals of the D. monticola population were the most distinct group. However, there were less obvious clusters present. A cluster consisting of mostly D. gregaria individuals with one D. behrii individual and one D. chryseopsis individual was isolated from the other three species. The remaining four D. gregaria individuals were intermingled with the other three species. The remaining clusters all showed a mixed representation of the three described species (D. basaltica, D. chryseopsis and D. behrii). The unrooted dendrogram calculated from population pairwise F_ST values showed a similar pattern. For some populations (ALP, WG, SK, CG and YP), F_ST values could not be calculated due to missing data and/or small population size, but the pattern observed in the unrooted dendrogram and the PCoA showed differentiation between D. gregaria and the other populations.

Fig. 3.

Phylogenetic representation of the ‘relaxed’ data set using genetic distance. Unrooted dendrogram calculated using PCoA on the ‘stringent’ data set. The inset shows the F_ST dendrogram from the ‘stringent’ data set. Some populations were removed from the F_ST dendrogram because the populations were either too small or had too much missing data.

The SNAPP analysis (Fig. 4; coalescent model) showed a similar pattern to those found with the genetic distance analyses. The D. monticola population, supported by 100 % of the trees, was separated from the main group with the greatest number of mutations (tree height) and was in agreement with the previously built Diuris phylogeny (Indsto et al., 2009). The only other structure detected by the coalescent model was the bifurcation between the D. gregaria populations and the other three species from the VVP, which was supported with 100 % of the trees. We considered this to be evidence of a separate genetic lineage. There was no structure between the remaining ten populations of D. basaltica, D. chryseopsis and D. behrii. The tree topology showed that the VVP populations were all connected within the 0·01 grid, meaning that there were few mutational differences between populations or the described species. Bayesian clustering analysis of the ‘stringent’ data set indicated that an optimal K-value of 2 was sufficient to explain the results via a strict Evanno interpretation (Fig. 5; Supplementary Data Fig. S5). We also showed a K-value of 3 because there was a secondary peak at three in the ΔK graph. Biologically, the K = 3 analysis showed finer structure within the VVP populations by separating D. gregaria populations from the other three species and showing admixture between D. gregaria and three of the D. chryseopsis populations. The DK population showed a nearly perfect admixture between the two VVP genetic clusters, which is in agreement with the separation of this population in the PCoAs.

Fig. 4.

A species complex tree obtained from a coalescent MCMC analysis of 97 individuals using SNAPP and BEAST on the ‘stringent’ data set. Support values are given for clades present in > 80 % of the iterations. Branch lengths are interpreted as relative time. Branches are identified by population abbreviations given in Table 1. Species are represented by colours.

Fig. 5.

Bayesian clustering analysis of the ‘stringent’ data set using (A) two genetic clusters (k = 2) and three genetic clusters (k = 3) in STRUCTURE. Each column (individual) is assigned a specific suite of colours corresponding to the membership coefficient. Populations are labelled according to Table 1. The line above STRUCTURE figures represents the five species (from left to right D. basaltica, D. gregaria, D. behrii, D. chryseopsis and D. monticola). (B) The posterior probability of the STRUCTURE analysis in geographical space. Genetic clusters are averaged across populations, and genetic cluster are analogous to Fig. 4A.

The log of genetic distance was plotted against the log of morphological distance with (Fig. 6A) and without (Fig. 6B) D.monticola. Two distinct clusters were present when D. monticola was included (Fig. 6A). The cluster to the right comprised D. monticola individuals whereas the cluster to the left comprised the four VVP species. The generalized linear model accounted for a small percentage of the variation in the data (r² = 0·023) but was significant (P = 2·1 × 10^–16). Here, genetic distance did have a weak relationship with morphological distance. When D. monticola was removed (Fig. 6B), even this weak relationship was lost, indicating that morphology was varying within a larger well-connected population.

Fig. 6.

Results from correlation analyses between the ‘relaxed’ data set [genotyping by sequencing (GBS)] and morphology for the Diuris species. Genetic distance vs. morphological distance along with a best-fit line from a GLM model, with D. monticola (A) where the grouping on the right contains the individuals from the outgroup and (B) without.

DISCUSSION

The five species we investigated accentuate the potential conservation risk and management shortcomings associated with reliance on species delineations based solely on non-exclusive morphological characters. Morphological variation within the group had a large amplitude and the characters used for species identification overlapped in the descriptions of currently recognized taxa (Fig. 2A;y-axis). On the basis of our results, morphological features alone were poor indicators of species divergence and not suitable to define conservation units. This is not surprising as Diuris is described as having few reliable morphological characters for classification (Indsto et al., 2009).

Our genetic data have provided evidence of some isolation between D. gregaria and the other three Diuris species found on the VVP. However, our findings did not confirm previously described groupings at the species level, nor did the findings correspond to population structure. There are important differences when comparing the F_ST between the microsatellite data sets provided here and in other studies. The F_ST among all four VVP species (0·03) is much lower than the average within-species F_ST (0·092) for 22 studies of rare orchid species (Phillips et al., 2012). This is indicative of a well-connected meta-population that should be conserved as such. Indeed, the genetic data indicated a high level of gene flow and a low level of inbreeding in the Diuris complex, including the nationally endangered D. basaltica populations. Current genetic structure suggests that gene flow has not been constrained significantly by reproductive, geographic or environmental isolation and is consistent with recent habitat fragmentation. Using the current species as our pre-defined populations, we would expect to find evidence of correlation between genetic differentiation and species delineations. Instead, our genetic results suggest that current species delimitations are not based on consistently different morphological characters, and are unsupported by genetic structure. Therefore, current species delimitations are not suitable as conservation units.

The measured divergence of D. gregaria implied some level of reproductive isolation despite sympatry with other members of the VVP group. However, the degree of differentiation between D. gregaria and the other VVP species is minimal because it represents only a tiny fraction of the total variation within the GBS data sets (4·16 %). This leads us to an important decision by which we must determine what demarcates an ESU and where the ESU threshold lies. In the case of this Diuris group, we are taking into account genetic diversity and current taxonomy. Therefore, for conservation purposes, we propose one ESU that includes D. gregaria and the other VVP species. We prioritize one ESU because we consider that the risk of negative consequences that can occur through isolation (unknown) is likely to be higher than through genetic connectedness (already established). In contrast to species growing on the VVP, D. monticola comprises a distinct genetic lineage despite its morphological similarity. This divergence may be a result of reproductive and environmental isolation between D. monticola and Diuris populations on the VVP. In its alpine environment, D. monticola flowers later than the VVP populations (November–January vs. September–October). As D. monticola occurs within the geographic range of D. behrii and D. chryseopsis in eastern Victoria, phenological data are required to determine the possibility of gene flow between the three species.

Evolutionary implications

When conservation focuses on one entity of a species complex, such as the nationally endangered D. basaltica, the downstream effects can include long-term negative impacts on evolution within the group (Ellstrand, 2014). By using taxonomically uncertain species as the conservation unit, land managers and conservation practitioners risk inadvertently prioritizing only a portion of the broader interactive gene pool. The Diuris populations to be conserved tend to be small, and thus are prone to negative impacts from stochastic processes such as genetic drift and non-random mating (Whiteley et al., 2015). As such, future speciation, species persistence and the rate of genetic change could be compromised by artificial genetic selection (Dynesius and Jansson, 2013). The outcome would be a skewed reduction in genetic diversity, particularly affecting small populations. A reduced gene pool is likely to be less resilient to the effects of rapid environmental change (Hoffman et al., 2015; Bragg et al., 2015). More generally, a meta-analysis of fragmented, inbred populations shows that outcrossing of inbred populations had benefits for 92·3 % of the cases included (Frankham, 2015). Frankham (2015) found that increased gene flow reduced inbreeding while increasing fitness effects and evolutionary potential. While not perfectly analogous to this situation due to low inbreeding of the Diuris group, the principles of genetic rescue remain germane because ongoing gene flow among populations within the group is critical for the maintenance of genetic diversity and evolutionary processes. In essence, any reduction in genetic diversity could have major ramifications on survival and fitness within the group. If the outcome of conservation strategies artificially restricts genetic connectivity by protecting few populations or those that contain only a sub-set of a larger well-connected metapopulation, overall evolutionary potential and long-term viability could be severely limited for the Diuris species studied here.

Conservation and management

When taxonomic uncertainty exists, extinction risks could be mitigated by managing gene pools more broadly to retain evolutionary processes, regardless of species boundaries. This inclusive strategy has the potential to simulate historical patterns of gene flow as well as minimizing risks associated with inbreeding, low genetic diversity and genetic drift. The benefit of maintaining allelic diversity for the group as a whole facilitates evolutionary processes (Sgrò et al., 2010) that could buffer the Diuris group from negative impacts from environmental change. High allelic diversity, as found in this study, facilitates genetic responses to an uncertain future and can aid conservationists in predicting evolutionary potential (Vilas et al., 2015). Given the rapid demographic decline but lack of distinct genetic lineages in the Diuris populations across the VVP, we propose that treating the group as one ESU is a less risky strategy than using taxonomically questionable species as conservation units.

With increased focus on active conservation and restoration (Barnosky et al., 2011), it is critical to incorporate genetic diversity and evolutionary potential into management practices. At the same time, funds for on-ground conservation work in the state of Victoria and within Australia more broadly are often predicated on the rarity of a species so conservation strategies must be created within that framework. We are unable to revise the taxonomic treatment of this Diuris group based on the data we have presented, but the challenging taxonomy of this group should not dissuade funding for the protection of the few remaining populations. A comprehensive revision of Diuris that includes the study taxa and considers the results presented here is needed. With this, we encourage that decisions made for one species should be made in the context of the others. Future reintroductions are planned for D. basaltica, and we believe that it would be in the best interest for the species to not be isolated from other species. Actions from this type of forward thinking could mitigate against the detrimental effects of genetic drift and other stochastic influences from environmental change. Our data suggest that anthropogenically assisted migration from any species on the VVP to the WG population is likely to be beneficial for longevity and genetic diversity.

Conclusions

When categorical species assignment cannot fully explain the variability of an organism, rigid adherence has the potential to result in weakened management outcomes. All organisms are exposed to changing environmental conditions and must migrate or adapt if they are to persist. High levels of genetic diversity that are subject to selection, gene flow and evolutionary processes provide the best strategy to optimize environmental resilience. This work illustrates a robust method that can be used for taxonomically difficult organisms to prevent uncertainty from impeding conservation. We envisage that this work will help to bridge the gap between theory and practice by identifying unique attributes of organisms or populations to maintain biological diversity and evolutionary trajectories. Ultimately, managing those attributes in the context of broader circumstances (e.g. climate change, fragmentation or higher levels of taxonomy) will lead to fruitful strategies that benefit biological diversity.

SUPPLEMENTARY DATA

Supplementary data are available online at https://academic.oup.com/aob and consist of the following. Figure S1: (A–F) Variation in floral morphology of six different Diuris individuals from the Parwan population. (G) Morphological variation measured during the study vs. that provided in descriptions. Figure S2: principal co-ordinate analysis of the 13 morphological traits. Figure S3: principal component analysis of 13 morphological measurements for all individuals in the study for five species. Figure S4: dendrogram from the ‘relaxed’ data set of the Diuris species complex on the Victorian Volcanic Plain in south-eastern Victoria with replicates annotated with arrows. Figure S5: line figure of Delta K from STRUCTURE Harvester.

ACKNOWLEDGMENTS

We thank N. Walsh, J. Jeanes and K. Brown for their valuable input. This work was supported by by Hansons Construction Materials Pty Ltd.