Morphological and molecular evidence for major re-circumscriptions in and eight new species of Melichrus R.Br. (Ericaceae subfam. Epacridoideae) in eastern Australia

Melichrus urceolatus has been maintained as a species through various treatments, albeit with differing limits. Paterson circumscribed M. urceolatus to include a vast geographical area; from Herberton, northern Queensland to Saint Arnaud, Victoria, acknowledging that ‘within this range a large variety of form is shown…’ (Paterson 1958, p. 313). Paterson (1958) provided a brief description of five groups that differed morphologically from the ‘typical condition’ and the approximate geographic range but did not attempt formal subdivision (Paterson 1958, p. 313). A great deal of this variation was noted to occur in Queensland, where three phrase-named putative species have been accepted by the APC (CHAH, accessed 4 May 2023), namely M. sp. Inglewood (A.R.Bean 1652), M. sp. Isla Gorge (P.Sharpe + 601) and M. sp. Tara (D.Halford Q2259).

Several putative species have been segregated from Melichrus urceolatus and informally phrase-named at the N.C.W. Beadle Herbarium (NE) but are not yet listed on the APC, including M. sp. Mareeba (H.T. Kennedy 92) NE Herbarium, M. sp. Galambary (H.T. Kennedy 7) NE Herbarium, M. sp. Gilgandra (H.T. Kennedy 128) NE Herbarium, M. sp. Gardens of Stone (H.T. Kennedy 46) NE Herbarium and M. sp. Salvator Rosa (M.B. Thomas 244) NE Herbarium that were identified during exploratory herbarium and field observations made in the early stages of this study.

Paterson (1958) retained the synonymy of Melichrus medius under M. urceolatus established by Bentham, explaining that no difference was observed in the protologue descriptions or subsequent collections made from the type localities. Improved access to the type specimen and other early collections allows the original definition of M. medius to be considered anew in this study.

Paterson reinstated Melichrus erubescens to species status, explaining that ‘the taxonomic differences are such as to justify specific rank, although the two species [M. urceolatus and M. erubescens] are undoubtedly closely related’ (Paterson 1958, p. 310). As currently delimited, the boundary between M. urceolatus and M. erubescens is unclear. Melichrus sp. Gardens of Stone, M. sp. Gilgandra and M. sp. Galambary, though treated as M. urceolatus by Paterson, all share morphological affinities with M. erubescens, exemplifying the unworkable boundary between the species as currently described.

Several populations of Melichrus that are morphologically similar to both M. erubescens and M. urceolatus have been collected for the first time since Paterson’s revision. These gatherings cannot be confidently assigned to either species and have been given the phrase names M. sp. Silent Grove (H.T. Kennedy 33) NE Herbarium and M. sp. Kaputar (H.T. Kennedy 25) NE Herbarium (Table 1). A key aim of this study was to decisively clarify the boundary between M. erubescens, M. urceolatus and morphologically similar putative species M. sp. Inglewood, M. sp. Isla Gorge, M. sp. Tara, M. sp. Gardens of Stone, M. sp. Gilgandra, M. sp. Galambary, M. sp. Silent Grove and M. sp. Kaputar.

The nomenclatural history of Melichrus adpressus is convoluted. The name M. adpressus was first applied to a collection made by Cunningham near Wellington Valley, New South Wales (De Candolle 1839). Later, Mueller listed M. adpressus as synonymous with Styphelia cunninghamii (von Mueller 1868, p. 39), which was described from specimens collected in northern New South Wales and southern Queensland. In doing so, an illegitimate name was erected and applied to two morphologically and geographically disjunct entities. Paterson (1958, p. 308) stated that:

the type locality of Melichrus adpressus is Wellington Valley but specimens subsequently collected from this area do not agree with the original description, and unfortunately Cunningham’s specimens were unavailable for comparison.

There is uncertainty regarding which specimens Paterson was referring to as ‘subsequently collected’ because no specimens from the Wellington Valley were cited in the revision. Paterson chose to adopt Mueller’s concept of M. adpressus and excluded the strict application of the name in the revision. Easy access to digital images of herbarium specimens (e.g. from JSTOR Global Plants) now greatly facilitates comparison of specimens that were largely inaccessible to earlier workers. See, for example, a specimen of M. adpressus collected by A. Cunningham near Wellington Valley that was seen by De Candolle, compared with a specimen cited by Paterson (Fig. 2). The morphological disparity between these specimens, including differences in leaf size, shape, margin texture, indumentum and density on stems, suggests that the name M. adpressus is misapplied in Paterson (1958). Herein, the name Melichrus adpressus sensu Paterson (or abbreviated to ‘s. Paterson’ in figures) indicates the use in the sense of von Mueller (1868) and Paterson (1958). When referred to in the original, strict sense of De Candolle (1839) this is denoted M. adpressus sens. str. (see Table 1).

Fig. 2.

Comparison of the application of the name Melichrus adpressus. (a) H. Beckler’s M. adpressus sensu Paterson specimen (MEL 2185645) from the northern rivers region of NSW, likely seen by F. von Mueller and cited as M. adpressus in Paterson (1958). Image reproduced with permission from Royal Botanic Gardens Victoria. (b) A specimen of M. adpressus sens. str. collected by Cunningham near Wellington Valley NSW that was seen by De Candolle (1839) and is a candidate for lectotypification (G 00454459). Image © Conservatoire et Jardin botaniques de la Ville de Genève.

Paterson (1958) noted only minor morphological variability across the range of Melichrus adpressus sensu Paterson. However, observations of herbarium and field specimens made early in this study suggested that this taxonomic consideration may neglect substantial morphological variation. We represent this variation by several phrase names. The northernmost populations of M. adpressus sensu Paterson are restricted to the Herberton area in Queensland and highly geographically disjunct from the southern populations in the state, and were designated as M. sp. Herberton (H.T. Kennedy 91) NE Herbarium. Specimens from the North Coast of New South Wales were morphologically distinctive and phrase named M. sp. Yuraygir (H.T. Kennedy 29) NE. Melichrus sp. Gurulmundi (H.T.Kennedy 69) NE Herbarium was applied to populations of M. adpressus sensu Paterson, found in southern Queensland that also share some morphological affinities with M. urceolatus. Isolated populations of Melichrus resembling M. adpressus sensu Paterson are known from slightly below the rocky outcrop at the summit of Mount Maroon, Queensland (Wahlmoorum of the First Nations’ Yuggera language), the upper slopes of Kroombit Tops National Park, Queensland and the sandy flats of Boonoo Boonoo National Park, New South Wales. The earliest collections of these populations were made in the latter half of the 20th Century, after the publication of Paterson’s (1958) revision. Each population is morphologically distinctive and to each we assigned the respective phrase names M. sp. Wahlmoorum (H.T. Kennedy 78) NE Herbarium, M. sp. Kroombit Tops (H.T. Kennedy 87) NE Herbarium and M. sp. Boonoo Boonoo (H.T. Kennedy 35) NE Herbarium for this study.

Melichrus procumbens was the first species of Melichrus described (as M. rotatus) and the taxonomic status of this morphologically distinctive species has never been questioned (Bentham 1868; von Mueller 1868; Paterson 1958). This species has a geographically disjunct distribution and has been recorded from coastal sands at Hervey Bay south to Sydney, where the type was collected and on usually granitic soils in the New England Tableland bioregion from Stanthorpe to Wollomombi. Paterson (1958, p. 307) described the plants across the disjunct range as ‘morphologically similar’ with the exception of a population near Glen Innes that was ‘in all respects larger … than typical specimens’ and hypothesised to be an autopolyploid. This putative new species was phrase-named Melichrus sp. Torrington (H.T. Kennedy 32) NE Herbarium. Since Paterson’s (1958) revision, a morphologically allied population was collected in Colo Gorge, New South Wales as M. procumbens and is here given the phrase name M. sp. Colo Gorge (I.R. Telford 8659) NE Herbarium. Sampling across the geographic range of M. procumbens is needed to test the relationship between morphologically dissimilar and geographically disjunct populations.

Melichrus hirsutus and M. gibberagee were each collected separately for the first time several years after Paterson’s (1958) revision was published. Despite the geographic proximity in north-eastern NSW, the taxonomic status was strongly supported by distinct morphological characteristics and these were published as new species (Kennedy et al. 2020). These two species are included in this study to test their status as distinct species based on molecular evidence.

A total of 90 individuals from 68 populations was included in the final morphometric and phenetic analyses; 97% of the sampled populations were also included in the DArTseq sample set, to maximise the complementarity of the datasets. A list of specimens and associated metadata is provided in Supplementary Table S1. Specimens were chosen to capture the full geographic and morphological range of each OTU but availability of suitable material was a limiting factor in some cases. Two OTUs, Melichrus sp. Colo Gorge and M. sp. Salvator Rosa, could not be located in the field. These were excluded from the analyses due to inadequate sampling. Wherever possible, a minimum of four samples per OTU was included. Several OTUs represented by fewer samples in the analyses were either geographically narrowly defined, such as M. sp. Kaputar, M. sp. Wahlmoorum, M. sp. Inglewood and M. sp. Tara, or appeared to be naturally rare and narrowly distributed; e.g. M. sp. Torrington.

During an initial exploratory phase, 205 characters were examined across 85 specimens, to represent a cross-section of the genus. This pilot yielded a total of 21 quantitative and 5 (binary) qualitative characters optimal for analysis (Table 2). The characters selected for statistical analysis were those that: (1) varied between samples, (2) could be scored reasonably easily and accurately, and (3) were selected with the aim of including as many characters as possible while prioritising quantitative characters and minimising missing data.

Morphological characters were scored either in the field or from photographic, dried and 70% ethanol-fixed specimens (Table 2). Macroscopic characters were measured using digital callipers or a ruler. Micromorphological characters were examined using a Nikon SMZ25 stereo microscope with a Nikon SHR Plan Apo 0.5× or 1× lens and measurements were made using Nikon NIS-Elements software. Descriptive erminology follows Radford et al. (1974) or if absent there, the ‘Flora of Australia’ glossary (McCusker 1999).

Analyses were conducted on a genus-wide dataset and on numerous subsets (Table 3). R (ver. 4.2.0, R Foundation for Statistical Computing, Vienna, Austria, see https://www.r-project.org/) and RStudio (ver. 2.3.4, Posit Software, PBC, Boston, MA, USA, see https://posit.co/products/open-source/rstudio/) were used for all analyses. To assess the suitability of the dataset for phenetic analysis, R package DataExplorer (ver. 0.8.2, B. Cui; see https://CRAN.R-project.org/package=DataExplorer) was used to plot the distribution of scores for each character and calculate the percentage of missing data by character, sample and for the entire dataset. The Pearson’s coefficient of correlation between each pair of characters was calculated. Where characters had a very high degree of correlation (<−0.8 or >0.8) and there was a biological rationale indicating that these may not have been mutually independent, one of the pair of characters was removed. Four ratio characters were included in the dataset. Where one or both components of the ratio were strongly correlated to the ratio character, one of the component characters was removed.

The character matrix was converted to a dissimilarity matrix using Gower’s distance coefficient (package cluster, function ‘daisy’, ver. 2.1.3, M. Maechler, P. Rousseeuw, A. Struyf, M.Hubert and K. Hornik, see https://cran.r-project.org/package=cluster). Gower’s distance coefficient is the most appropriate coefficient for use with datasets comprising a mixture of quantitative and qualitative characters, and a small amount of missing data that are weighted as zero in analyses (Gower 1971). Phenetic methods in which qualitative characters have been used have received criticism, as these can potentially have a disproportionate influence on the result. To better understand the influence of each qualitative character on the results, analyses were run numerous times, using quantitative characters only initially and subsequently adding the qualitative characters in one by one and interrogating the influence on sample distributions (Supplementary Fig. S1).

Sample similarity was analysed and visualised in two complementary ways. A hierarchical clustering phenogram was generated (R package stats, function ‘hclust’, ver. 4.4.2, see https://stat.ethz.ch/R-manual/R-devel/library/stats/html/00Index.html) using the unweighted pair group method with arithmetic mean (UPGMA). Non-metric multidimensional scaling (NMDS) was performed in vegan (function ‘metaMDS’, ver. 2.6-2, J. Oksanen, G. L. Simpson, F. G. Blanchet, et al., see https://CRAN.R-project.org/package=vegan) to visualise clustering patterns in the data using 2–4 dimensions (k). UPGMA and a non-metric ordination algorithm are appropriate metrics for this dataset, as most characters were not normally distributed (Supplementary Fig. S2). The NMDS was run with 30 random starts for k = 2–5. This protocol was determined experimentally to be sufficient to reach a convergent solution for all values of k. To better understand how many dimensions are needed to best visualise the dataset, the ordination distance was plotted against the observed dissimilarity (function ‘stressplot’). A stress value is considered ‘fair’ if <0.2 and ‘good’ if <0.1 (Kruskal 1964). In addition to considering the stress inherent in each ordination, ordinations were visually inspected and where adding an additional dimension did not change the pattern observed, the lower value of k was preferred.

Ordination was performed for subsets 1–7 (Table 3). The subsets were based on the clustering pattern observed in the genus-wide ordination and phenogram, and designed to examine fine-scale morphological patterns in the data.

The phenogram and ordinations were visually inspected for clustering patterns that were potentially representative of taxonomic groupings. Additionally, the mclust program (function ‘mclust’, ver. 5.4.9, see https://CRAN.R-project.org/package=mclust; Scrucca et al. 2016, 2023) was used to investigate clustering patterns that best represent the data under a series of models (Fraley and Raftery 2002; Meudt et al. 2020). mclust chooses among 14 different models describing the shape, volume and orientation of clusters. For a given number of clusters, G, ranging from 1 to 9, mclust will choose the best-fitting model based on the Bayesian information criterion (BIC). mclust also allocates individual observations to each cluster under the optimal model and assesses the uncertainty of the classification for each specimen. Only results from the best-fitting model are shown.

To be potentially recognised as a cluster in mclust analyses, an OTU must have at least one more specimen sampled than the number of dimensions used (Hausdorf and Hennig 2010; Edwards and Knowles 2014). In our dataset the OTUs that did not fulfil this requirement at k = 2 and k = 3 were M. sp. Torrington (n = 2), M. sp. Inglewood (n = 2), M. sp. Kaputar (n = 2), M. sp. Tara (n = 2), M. sp. Wahlmoorum (n = 1) and additionally at k = 3, M. adpressus sens. str. (n = 3).

The influence of each character on the distribution of points in the ordination was examined to help identify diagnostic characters for putative species. The ‘envfit’ function in vegan was used to fit the quantitative characters to the ordination and those with a significant influence (P < 0.05) were visualised in a biplot with arrow length proportional to the correlation between the data for that character and the ordination configuration. The direction of each arrow indicates the axis of the ordination space across which that character changes most rapidly. Examining the magnitude and orientation of each arrow provides insight into potentially diagnostic characters for OTUs. For qualitative characters, the average value for each state was plotted on the ordination to examine the influence of the character on the distribution of samples.

Ordination analysis was replicated for subsets 1–7 (Table 3) to search for potentially diagnostic characters between morphologically similar OTUs that were not detectable at the coarser, genus-wide scale. Subsets were based on the clusters apparent in the genus-wide phenogram and mclust analysis of the NMDS ordination (clusters 1–7 in Fig. 3; subsets summarised in Table 3). These subsets were not subjected to mclust analysis as sample numbers were too low for this analysis to be effective, therefore the NMDS ordinations were visually inspected for clustering patterns. Cluster 8, as defined in the k = 3 genus-wide ordination (Fig. 3c), was not considered suitable for ordination due to inadequate representation of M. sp. Wahlmoorum in the analysis. Cluster 6 also had a low number of samples that should be considered when interpreting NMDS results for this subset.

Fig. 3.

Morphometric analysis of 90 specimens of Melichrus scored for 26 characters showing NMDS ordination and mclust (grey ellipses) results. (a) k = 2, model ‘EII’, best number of clusters = 8. (b, c) k = 3, model ‘EEI’, best number of clusters = 9, axes 1 & 2 and axes 1 & 3. (d) A hierarchical clustering phenogram. Blue bars indicate clusters that correspond between the ordinations and phenogram, yellow bars indicate further clusters that are apparent only in the phenogram. Points represent individual plants. See Fig. 1 for a visual key and Table 1 for a textual description of symbols.

A table of differences was constructed to highlight diagnostic characters (Table 4a–c), using characters that varied between OTUs but that were not captured in the phenetic analyses. These characters were scored for the same specimens used in the phenetic analyses.

A total of 572 dried leaf samples were submitted to Diversity Arrays Technology Pty Ltd (DArT) for extraction of total genomic DNA (Supplementary Table S1). In addition to the 548 samples from 110 populations of eastern Australian Melichrus, this included 12 samples representing the three western Australian phrase-named entities, ten outgroup samples from other genera of Styphelieae (for use in phylogenetic analyses reported on in a future publication) and two technical duplicates (that were found to be surplus to the requirements of the analyses). A single sample, obtained from the only collection of Melichrus sp. Bruce Rock housed at the Western Australian Herbarium (PERTH), failed to produce sequence data.

Single nucleotide polymorphism (SNP) genotypes were determined from sequence reads by DArT using proprietary DArTseq pipelines (Kilian et al. 2012). These pipelines filter and remove poor quality sequence reads before de novo assembly of loci and calling of genotypes. Technical replicates of multiple samples are included and consistency of genotype calls between the replicates scored to provide a ‘reproducibility’ score for each locus. SNP genotypes were called on the whole sample set and subsequently called on a separate partition of samples that included only eastern Melichrus and a reduced selection of outgroup samples. Exploration of the two datasets revealed very little difference in the quality of SNPs between datasets and no difference in the results of analyses. We report only on the original, pre-partitioning sample set, as this returned a larger number of SNPs (95,319 SNPs) than the reduced partition (82,287 SNPs).

The statistical software R (ver. 4.2.0) and RStudio (ver. 2.3.492) were used for all statistical analyses. The package dartR (ver. 1.9.9.1, see https://CRAN.R-project.org/package=dartR; Gruber et al. 2018; Mijangos et al. 2022) was used to manipulate, filter and analyse the data. Only the eastern Australian Melichrus samples were required for this study, therefore the outgroup taxa and two technical duplicates were removed. Any resulting monomorphic loci were removed and the dataset metrics recalculated. The resultant dataset consisted of 548 individual samples and 78,474 SNPs, 76% of which had a reproducibility of 100% (93% of SNPs had a reproducibility of ≥98%). Approximately 24% of SNPs (18,940 SNPs) had <10% missing data and the average call rate for samples was 76%.

Initial exploratory filtering of DArTseq data at different levels of stringency confirmed only minor effects on the clustering of samples in PCA or SplitsTree4 (ver. 4.1, see https://github.com/husonlab/splitstree4; Huson and Bryant 2006). Refined sets of filtered SNPs were subsequently used in all analyses for improved calculation efficiency (Rutherford et al. 2018, Collins et al. 2022). For the entire dataset, all SNPs with a reproducibility of less than 98% were removed (gl.filter.RepAvg) and all remaining loci with more than 5% missing SNP calls were subsequently removed (gl.filter.callrate). To prevent linkage disequilibrium from influencing analyses, one SNP was randomly selected to retain from each of the loci that contained more than one SNP (gl.filter.secondaries). Any individual sample with more than 10% missing data was removed (gl.filter.callrate), except that for data subset F (Table 3), the filter was relaxed to allow up to 11% missing data to retain the single specimen of M. sp. Wahlmoorum in the analysis. SNPs where the frequency of the minor allele (MAF) was less than 0.05 were removed (gl.filter.maf). Finally, any monomorphic loci (gl.filter.monomorphs) were removed and the data metrics recalculated (gl.recalc.metrics). This filtration protocol was used on all sets of data for all analyses except for the population statistics. The filtered dataset used for population statistics was not filtered for minor allele frequency as this may affect estimates of observed heterozygosity (Linck and Battey 2019).

To improve resolution among OTUs and the associated populations in clustering analyses, several subsets of the dataset were made (summarised in Table 3). The dataset was divided into subsets based on the clusters observable in the results of a genus-wide PCA and the major clades in a Neighbour-Net analysis. Subsets were derived from the unfiltered in-group only dataset and subsequently filtered following the procedure outlined above.

PCA (dartR function ‘gl.pcoa’) was used to explore genetic structure between putative species and populations. PCA was performed on the entire dataset and all subsets of the data outlined in Table 3.

SplitsTree4 was used to conduct a Neighbour-Net analysis (Bryant and Moulton 2004) producing a distance-based network of relatedness for samples. The network can represent evolutionary histories with substantial reticulation or discordance, caused by incomplete lineage sorting and introgression (Huson and Bryant 2006). A relationship network of the total dataset was computed using the default software settings.

A Bayesian-model based procedure implemented in STRUCTURE (ver. 2.3.4, Pritchard et al. 2000) was used to further assess population structure and individual ancestry. STRUCTURE assigns individuals or a proportion of an individual’s ancestry, when using an admixture model, to K genetic groups based on allele frequencies for multi-locus genotypes without any prior assumptions of population, geographic or ancestral groups.

An admixture model with correlated allele frequencies and no location priors was used. Estimation of ancestry was performed on all the molecular subsets listed in Table 3. The ANCESTDIST option was used to return 90% probability intervals around each ancestry estimate. Each run included 400,000 MCMC steps and a burn-in period of 100,000 steps. Ten independent replicate runs were conducted for each value of K, for at least K = 1–10 and up to K = 13 for larger subsets.

Determining the best-supported number of genetically distinct groups (K) can be challenging. Several biases in a sample set, such as population divergence caused by isolation-by-distance, subtle subdivisions within diverged groups and recent genetic bottlenecks can lead to misinterpretation of results (Lawson et al. 2018). We took a cautious approach to identifying the best supported values of K, following the recommendations of Janes et al. (2017), Lawson et al. (2018) and Funk et al. (2020). This involved a two-pronged approach. We closely examined all bar plots produced by each STRUCTURE analysis at different values of K (from multiple runs and iterations), seeking consistency among results. We also plotted the natural logarithm of the likelihood (ln Pr(X/K)) for all values of K tested, using STRUCTURE HARVESTER (ver. 0.6.94, see http://taylor0.biology.ucla.edu/structureHarvester, accessed 30 September 2023). When interpreting the STRUCTURE results, we described the ancestral assignment for several values of K in relation to known natural history data (Funk et al. 2020). Admixture plots were generated using CLUMPAK (beta version, see http://clumpak.tau.ac.il/; Kopelman et al. 2015); all admixture plots are included in the Supplementary material (Supplementary Fig. S6–S9). CLUMPAK identifies distinct clustering modes from multiple independent runs at each value of K in the space of possible solutions for that value of K, and the runs are sorted into the predominant or ‘major’ clustering mode and any alternative or ‘minor’ clustering modes that occur less frequently. CLUMPAK also aligns inferred clusters across runs with different values of K for ease of plotting.

The R package conStruct (ver. 1.05, see https://cran.r-project.org/package=conStruct; Bradburd et al. 2018) was used to investigate the role of isolation by distance in producing some of the observed clustering patterns in the STRUCTURE and PCA results. conStruct is a model-based clustering method similar to STRUCTURE but with several important differences. When assigning ancestry, conStruct models ancestral populations as spatial layers incorporating isolation by distance to explain observed genetic variation, only ascribing discrete population structure to the dataset when the genetic differentiation significantly deviates from that explained by geographic distance. conStruct also provides a non-spatial model similar to STRUCTURE for comparison with the spatial model. When combined with STRUCTURE analyses conStruct provides a robust method to test whether sampled variation is discrete or continuous across a landscape. conStruct uses population-based allele frequencies, rather than individual genotypes as used by STRUCTURE and incorporates geographic information. Only populations with reliable latitude and longitude information available, and population level sampling of three or more individuals were used in this analysis, resulting in 16 populations being excluded on this basis. Only one OTU was entirely lost from the analysis through filtration, M. sp. Wahlmoorum, otherwise all OTUs retained adequate representation.

Basic runs (function ‘construct’) were performed on the specified data subsets. Initially the number of iterations, replicates and value of K, were explored and 20,000 iterations for 10 replicates of each value of K from 1 to 10 were found to produce results that converged largely on one solution and for which the most trace plots showed little autocorrelation and had settled into a solution space (see vignette ‘model-comparison’; Bradburd et al. 2018).

To help determine the best value of K for each subset analysed, cross validation (function ‘x.validation’) analysis was performed using the same run parameters as for basic runs. We plotted the mean predictive accuracy for each value of K to find the value(s) for which the data best fitted the model, and to compare the spatial and non-spatial models. To avoid over-parameterisation, we also calculated layer contributions (function ‘layer.contribution’) for each value of K. This analysis provides insight into how much each layer contributes to the total covariance; if the contribution of a layer is small to negligible it indicates that fewer layers are likely to be equally descriptive. Similar to STRUCTURE analysis, considering multiple values of K and placing the results within a biological framework when interpreting the implications of the results are important.

Genetic diversity parameters were estimated for each population, with three or more individuals sampled, to provide insight into gene flow within and between populations, and to look for evidence of hybridisation, polyploidy and apomixis. The expected heterozygosity (H_E), observed heterozygosity (H_O) and degree of inbreeding (F_IS) were estimated using the hierfstat package (function ‘basic stats’, ver. 0.5-11, J. Goudet and T. Jombart, see https://cran.r-project.org/package=hierfstat).

To investigate genetic differentiation between all populations, the pairwise F_ST was estimated using STAMPP (function ‘stamppFst’, ver. 1.6.3, see https://cran.r-project.org/package=StAMPP; Pembleton et al. 2013). Where pairwise F_ST values between taxa are high relative to intraspecific pairwise values within congenerics, while maintaining high levels of outcrossing (low F_IS), these can be considered genetically differentiated (Rutherford et al. 2018).

Two characters, ‘bud width’ and ‘corolla lobe width’ were removed from the analysis as these were highly correlated with several other characters and there was a biological rationale to suggest a lack of independence (Supplementary Table S2). An analysis of only the quantitative characters showed the basic structure of an emerging pattern in the data. However, several clusters in the phenogram and ordination were not clearly resolved, including OTUs from the M. erubescens sensu Paterson and M. urceolatus sensu Paterson groups. Qualitative characters were included one at a time and analyses rerun (Supplementary Fig. S1). The underlying pattern in the data did not change with the addition of each character but the distinction between clusters became clearer. After filtration, a dataset comprising 90 samples (Supplementary Table S1), 21 quantitative variables and 5 binary characters (Table 2) was deemed suitable for the final analysis.

The total amount of missing data in the genus-wide final dataset was 3.3%. Corolla characters had the most data missing (up to 8.9% per character; Supplementary Table S3a) due to the practical difficulty of aligning fieldwork with the period of flowering. All the missing data were contributed by 10 individual samples, from six putative taxa, for which incomplete fertile material was available (Supplementary Table S3b). The analyses were run both with and without these samples and the sample inclusion did not influence the general patterns observed. These were retained in the final analysis to boost representation for otherwise under-represented taxa.

Most quantitative characters were not normally distributed. Some characters were bimodal and many showed a skewed distribution (Supplementary Fig. S2). This supported the justification for use of a non-metric ordination approach and the UPGMA hierarchical clustering method.

The pattern of clustering in the phenogram (Fig. 3d) was consistent with that visualised in the NMDS ordinations (Fig. 3a–c). In Fig. 3d, clusters that correspond between the phenogram and NMDS ordinations are indicated by the blue bars, whereas further clusters apparent only in the phenogram are shown with yellow bars. The clustering pattern reflected many of the species boundaries hypothesised early in this study, expressed as phrase-named OTUs (Table 1). The NMDS ordination for both k = 2 and k = 3 returned fair stress values of 0.1920 and 0.1062 respectively, whereas k = 4 returned a good stress value of 0.0799. The clustering pattern produced by each analysis was similar and iterative, with the addition of more dimensions improving the resolution of clusters. Results for k = 2 and k = 3 NMDS ordinations are shown (Fig. 3a–c). The mclust model selected as the best model to represent the genus-wide data for k = 2 was ‘EII’ (spherical, equal volume) and the best number of clusters was eight. For k = 3, ‘EEI’ (diagonal, equal volume and shape), the optimal number of clusters was nine. For k = 3, eight of the nine clusters identified were also recovered in the k = 2 analysis (Fig. 3a–c).

The influence of each character on the clustering patterns in the ordination was examined to better understand what was driving the clustering pattern and identify potentially diagnostic characters for putative species. Results for the k = 3 ordination (Supplementary Fig. S3) were concordant with the k = 2 results (Fig. 4). All 26 characters had a statistically significant impact on the ordination configuration for k = 3, providing some further insight into diagnostic characters, especially for Melichrus gibberagee.

Fig. 4.

Character importance analysis of an NMDS ordination (k = 2) for a morphometric dataset comprising 26 characters scored for 90 specimens of Melichrus. Points represent individual plants. See Fig. 1 for a visual key and Table 1 for a textual description of symbols. Each vector represents the magnitude and direction of the influence of a quantitative character that had a statistically significant influence on the ordination space. The average value for each state (Y, affirmative; N, negative) for each binary character is plotted, the corresponding character for numbers 1–5 are listed in Table 2. Abbreviations: C., corolla; L., length; W., width.

The distribution of samples across axis 1 primarily reflected the fact that taxa in clusters 1 and 2 have a cup-shaped or rotate corolla (C. tube L./C. lobe L.; C. tube W.; 2N), with corolla lobes that are spreading from the base (1N; e.g. see Fig. 5), relatively narrow leaves (Leaf L./W.) and young leaves that are straight or slightly incurved (4N; e.g. see Fig. 5a, b). By contrast, taxa in clusters 6 and 7 have an urceolate corolla, with lobes that are held erect before spreading, leaves are generally ovate and young leaves are recurved (Fig. 4).

Fig. 5.

Selected morphological characters in Melichrus. (a) Young leaves recurved (M. urceolatus, H.T. Kennedy 130). (b) Young leaves straight to slightly incurved (M. sp. Galambary, H.T. Kennedy 136). (c) Fruit elliptical, sepal apex acuminate and exceeding fruit (M. gibberagee, H.T. Kennedy 28). (d) Fruit ovoid, sepal apex obtuse and not exceeding fruit (M. adpressus sens. str., H.T. Kennedy 140). (e) Fruit spheroid, sepal apex acute and exceeding fruit (M. urceolatus, H.T. Kennedy 61). (f) Corolla urceolate (M. urceolatus H.T. Kennedy 130). (g) Corolla cup-shaped (M. sp. Galambary, H.T. Kennedy 137). (h) Corolla rotate (M. procumbens, H.T. Kennedy 123). Images by J. J. Bruhl (c & e) and H. T. Kennedy (a, b, d, f, g & h).

Cluster 3 occupies a similar area of the ordination space to cluster 2, and these share several characters, including a broad sepal apex (Sepal apex angle), a domed habit (5Y) and a prominent central nerve on the abaxial leaf surface (3Y). OTUs in cluster 3 differ from those in cluster 2 by the broader leaves and some samples have an urceolate corolla shape (2Y).

Several characters contributed to the pronounced isolation of cluster 5 in the ordination space (Fig. 4). These characters included dense, imbricate leaves (Leaf density), robust branchlets (Branchlet diameter) and a wide lamina margin (Lamina margin W.; e.g. see Fig. 6a–c). The maximum length of hairs on the lamina margin (Lamina margin hair max. L.) strongly distinguished clusters 5, 8 and 1 across axis 2 from taxa with much shorter lamina hairs in clusters 2, 3, 6 and 7 (e.g. see Fig. 6a–c).

Fig. 6.

Selected micromorphological characters in Melichrus. (a) Upper leaf surface showing narrow coloured margin, minutely scabrous (M. adpressus sens. str., H.T. Kennedy 11). (b) Upper leaf surface showing wide, white, thin margin, scabrous to hispid (M. adpressus sensu Paterson, H.T. Kennedy 82). (c) Upper leaf surface margin not distinct from leaf surface, hairs pilose (M. hirsutus, H.T. Kennedy 112). (d) Lower leaf surface showing many veins with narrow furrows between veins that have minute barbs on the margins (M. sp. Herberton from QCO, H.T. Kennedy 91). (e) Lower leaf surface showing medium number of veins with wide furrows between veins that have short papillae on the margins (M. sp. Mareeba, H.T. Kennedy 94). (f) Lower leaf surface showing few veins with very wide furrows between veins that have long papillae on the margins (M. sp. Isla Gorge, H.T. Kennedy 84). Scale bars: 1 mm. Images: H. T. Kennedy.

Clusters 6 and 7 shared a similar area of the ordination space reflecting shared floral characters and that the young leaf growth is often recurved (4Y). Taxa in cluster 6 are distinguished from cluster 7 in the analyses by the especially small, sparse leaves (Fig. 4).

Cluster 4 is the most broadly spread cluster identified in the analysis, suggesting that this has the most intracluster differences (Fig. 4). This cluster also occupies the middle of the ordination, indicating a lack of unique characters captured in the analysis for taxa in this group. Table 4b highlights some additional characters that distinguish OTUs in this cluster from the proximal, cluster 5. Particularly, OTUs in cluster 4 all share a character in which the furrows between the veins on the undersides of leaves are particularly wide and have long hairs on the margins of the furrows (e.g. see Fig. 6d–f). Melichrus hirsutus (cluster 1) and M. gibberagee (cluster 9) are the only other taxa that share this character.

The k = 3 ordination split cluster 8 (from the k = 2 ordination) into two groups separated by significant morphological differences (Supplementary Fig. S3). Taxa in both clusters 8 and 9 have long, urceolate flowers, long buds and sepals with a narrow apex. However Melichrus gibberagee has far narrower leaves (Leaf L./W.) and flowers (Corolla lobe L. & Corolla lobe L./corolla lobe W.), far fewer veins on the abaxial leaf surface (Lamina veins count) and longer lamina margin hairs (Lamina margin hair max. L.) than in M. sp. Wahlmoorum and M. adpressus sensu Paterson. Further differences between the three OTUs are included in the table of differences (Table 4c).

NMDS ordination of the subsets of the data (see Table 3 for a summary of data subsets analysed) recovered the same additional clusters that were observable in the genus-wide phenogram, indicated in Fig. 3d by yellow bars. Supplementary Fig. S5 includes the most informative NMDS ordinations and character importance analyses for each data subset, the results of which are summarised here.

The NMDS ordination of a subset of the data including only OTUs from cluster 1 (subset 1), distinguished strongly between Melichrus hirsutus and M. procumbens plus M. sp. Torrington (Supplementary Fig. S4a). Several characters had a significant influence on the ordination pattern, particularly the lack of a white lamina margin, taller stature, cup-shaped corolla, more slender buds, acuminate sepal apex and wider leaves in M. hirsutus. Table 4a highlights several more characters that morphologically differentiate M. hirsutus from other OTUs in cluster 1, including an ascending habit, lamina margin hairs are all a similar length or longer at the base of the lamina and a pink corolla.

Ordination of subset 2 showed Melichrus erubescens, M. sp. Silent Grove and M. sp. Gardens of Stone samples clustering in three discrete groups corresponding with OTU designation (Supplementary Fig. S4b, c). Samples of M. sp. Gardens of Stone were distinguished by a relatively long, narrow, cup-shaped corolla with long corolla lobes. Melichrus sp. Silent Grove was distinguished in the ordination by a tall, erect stature and narrow corolla tube with relatively short corolla lobes. Melichrus sp. Gardens of Stone and M. sp. Silent Grove are further distinguishable from M. erubescens by acuminate sepals at fruiting that are longer than the fruit. Melichrus sp. Silent Grove also has a different fruit texture, similar to that seen in M. urceolatus where there is a distinctly smooth, slightly raised distal section (Table 4a).

Morphological analyses did not return any reliable diagnostic characters differentiating the three morphologically similar OTUs Melichrus sp. Kaputar, M. adpressus sens. str. and M. sp. Galambary (Supplementary Fig. S4d, e). The results did suggest that some of the southern populations of M. sp. Galambary (Gal3 & 4) were morphologically different from the more northern populations. In particular, southern populations had a longer corolla tube that was subtly urceolate with lobes that were held erect before spreading, the leaves were generally narrower and the dome-shaped habit was lacking.

Melichrus sp. Mareeba specimens in the ordination of subset 4 were differentiated primarily by the distinctive leaves that were larger and had a broader white margin, similar to the OTUs in cluster 5 (Supplementary Fig. S4f, g). Melichrus sp. Isla Gorge specimens clustered loosely in the ordination, demonstrating the variability for the chosen characters. The sample Isl5 of M. sp. Isla Gorge is nested within the M. sp. Mareeba samples in the genus-wide phenogram and M. sp. Isla Gorge and M. sp. Mareeba samples were not well separated in the ordination. Melichrus sp. Mareeba can be distinguished from M. sp. Isla Gorge and M. sp. Tara by several characters not included in phenetic analyses, such as a decumbent plant habit, the absence of pilose hairs on leaves and glabrous sepals (Table 4b). The two M. sp. Tara specimens included in the analysis were distinguished in the ordination by a broad sepal apex angle and short buds and flowers.

Samples in subset 5 separated into three distinct groups across the ordination space (Supplementary Fig. S4h). Samples of Melichrus sp. Herberton clustered together, defined by broader leaves and sepal apex angle and as the only OTU in this cluster with glabrous sepals (Table 4c). Melichrus sp. Yuraygir samples showed a high degree of variability in some corolla and indumentum characters but clustered discreetly based on the wide white margin, high density of short, wide leaves and a short pungent tip at the leaf apex. Samples of M. sp. Boonoo Boonoo and M. sp. Kroombit Tops occupied a similar part of the ordination, though these were not tightly clustered. This part of the ordination space was defined by longer, slightly narrower leaves and generally longer hairs across various plant parts.

Several vegetative and fertile characters separated Melichrus sp. Gilgandra and M. sp. Inglewood in the ordination of subset 6 (Supplementary Fig. S4i). Melichrus sp. Gilgandra samples generally had a wider lamina with a longer pungent tip, and flowers with a broad sepal apex and relatively short corolla lobes. Additionally, M. sp. Gilgandra leaves are strongly glaucous, the texture of the fruit is variably wrinkly all over and the calyx does not exceed the apex of the fruit (Table 4c).

All samples of Melichrus sp. Gurulmundi shared a similar area of the ordination for subset 7 (Supplementary Fig. S4j, k). This area was defined by a high density of leaves on flowering branches, a broad branchlet diameter, broad sepal apex angle and more veins on the abaxial leaf surface. Some samples of M. urceolatus, particularly Urc4 and Urc6 also fell within this area of the ordination, in evidence of the similarity that M. sp. Gurulmundi has with some of the more robust forms of M. urceolatus. Although M. sp. Gurulmundi is morphologically similar to M. urceolatus, these can be distinguished by a combination of characters. Melichrus sp. Gurulmundi rarely has more than one stem at the base, is robust, erect, with a candelabra-like branching pattern, has densely imbricate leaves, and sepals are pilose and have a broad apex angle (Supplementary Fig. S4j, k and Table 4c). Samples of M. medius were not distinguished from M. urceolatus in the phenetic analyses. The table of differences highlights the particularly densely pilose sepals but some other populations of M. urceolatus also exhibit this character, particularly populations on the Northern Tablelands of New South Wales (Table 4c).

Filtration of the dataset for PCA and Neighbour-Net analysis produced a set of 2071 high-quality SNPs and retained 543 individuals. No entire populations or OTUs were lost from the analysis through filtration.

The first three axes of the ‘genus-wide’ PCA (Fig. 7a, b) explained 20.6, 18.2, and 11.4% of the total variation respectively. OTUs segregated from Melichrus urceolatus sensu Paterson (represented by circles) separated into at least five distinct groups occupying the extremities of all axes (Fig. 7 clusters 1, 4, 5, 7 & 9). Cluster 1 included M. urceolatus, M. medius and M. sp. Inglewood, and M. sp. Silent Grove. In further analyses this is referred to as the ‘M. urceolatus cluster’. All OTUs attributed to M. procumbens sensu Paterson (1958) (represented by diamonds) clustered closely together (cluster 2; ‘M. procumbens cluster’) adjacent to the M. urceolatus cluster (cluster 1) and cluster 3 that contained all M. adpressus sensu Paterson (‘M. adpressus sensu Paterson cluster’) segregates (represented by squares). Note that the distinctions between clusters 1–3 are more easily observed in the network analysis and PCA of subset A (see below). Cluster 4 included all six M. sp. Mareeba populations, and cluster 5 included all eight M. sp. Isla Gorge and the single M. sp. Tara populations (‘M. sp. Isla Gorge cluster’). Melichrus hirsutus and M. gibberagee clustered closely together (cluster 6; ‘M. hirsutus cluster’). Axis 3 defined clusters corresponding to OTUs or groups thereof: cluster 7, consisting of M. sp. Galambary, M. sp. Gilgandra and M. adpressus sens. str. (‘M. sp. Galambary cluster’); cluster 8, comprising M. sp. Kaputar (except one individual that is displaced from the group); cluster 9, with M. sp. Gardens of Stone; and cluster 10, consisting only of M. erubescens samples (represented by triangles).

Fig. 7.

Melichrus genus-wide PCA and Neighbour-Net analysis of a DArTseq dataset of 2071 SNPs. (a) PCA axes 1 & 2. (b) PCA axes 1 & 3. (c) SplitsTree Neighbour-Net. In the ordinations (a, b) individuals are represented by single points, in the SplitsTree network (c) each terminus represents a single population. Point shape indicates membership to a sensu Paterson group. Numbered clusters on the SplitsTree network correspond with number clusters on the PCA ordinations. See Fig. 1 for a visual key and Table 1 for a textual description of symbols.

The groups observed in the SplitsTree Neighbour-Net analysis (Fig. 7c) mirrored the PCA, with some additional insights provided by the distance-based relatedness network. In a relationship network, the parallel lines reflect splits in the data, with the longer lines indicating more support for that particular split (Huson and Bryant 2006). The split separating the Melichrus erubescens cluster from the remainder of the genus is especially strong in this analysis, as are the splits defining the M. procumbens (cluster 2), M. sp. Isla Gorge (5) and the M. sp. Galambary (7) clusters. The network confirms that which is not necessarily distinct in the PCA, that the M. urceolatus (1), M. procumbens (2) and M. adpressus sensu Paterson (3) clusters are distinct groups. The higher number of parallel edges in the sections corresponding with the M. urceolatus (1) and M. adpressus sensu Paterson (3) clusters respectively indicate a larger degree of discordance in these sections of the dataset. The M. sp. Mareeba cluster (4) shares multiple parallel lines with both M. sp. Isla Gorge (5) and M. sp. Herberton + M. sp. Kroombit Tops (3) that may indicate a hybrid ancestry for M. sp. Mareeba. Parallel edges between M. sp. Silent Grove and M. procumbens may indicate some introgression, as do the parallel lines between M. sp. Kaputar and M. sp. Galambary.

PCA of subset A, that includes only the samples from the genus-wide PCA clusters 1–3 (see Table 3 for a summary of subsets used in analyses; Fig. 8a), exhibited three broad clusters distributed across axes 1 (28.4%) and 2 (14.1%). Each group matched with the Melichrus urceolatus (1), M. procumbens (2) and M. adpressus sensu Paterson (3) clusters as seen in the Neighbour-Net analysis.

Fig. 8.

The most informative principal component axes (PCs) plotted for select subsets of a DArTseq dataset of 2071 SNPs. (a) Subset A; PC1, PC2. (b) Subset E; PC1, PC2. (c) Subset D; PC1, PC2. (d) Subset J; PC1, PC2. (e) Subset F; PC1, PC2. (f) Subset F; PC1, PC3. (g) Subset I; PC1, PC2. (h) Subset I; PC1, PC3. (i) Subset H, PC1, PC2. (j) Subset G, PC1, PC2. Individuals are represented by single points. See Fig. 1 for a visual key and Table 1 for a textual description of symbols. The OTUs included in each subset are listed in Table 3. Codes used to indicate geographic distributions are abbreviations of botanical districts and follow (Centre for Australian National Biodiversity Research 2018), except ‘nNNT’ and ‘sNNT’ that are modified to indicate the northern and southern extents of the Northern Tablelands NSW (see Supplementary Table S1 for further geographic information).

In a PCA of the Melichrus procumbens cluster (subset E), the only axis that explained more than 10% of the variation in the dataset was PC1 (23.5%). Along this axis, the two populations of M. sp. Torrington were distinct from M. procumbens samples. Melichrus procumbens showed some geographic grouping corresponding with coastal and tableland populations (Fig. 8b).

PC1 (16%) from the PCA of subset D, defined a ‘northern cluster’ (Fig. 8c) containing Melichrus urceolatus specimens collected north of Henry River, New South Wales in the Northern Tablelands (nNNT) and Queensland Darling Downs (QDD), and the single populations of M. sp. Silent Grove and M. sp. Inglewood respectively. Along PC2 (10.4%) three additional clusters were apparent: the first containing all M. medius (Southern Tablelands, NST) populations; the second, populations collected on the Northern Tablelands south of Henry River (sNNT), the third contained M. urceolatus samples from the North Western Slopes (NNS), Central Tablelands (NCT), Central Coast (NCC) and South Coast (NSC) botanical districts of New South Wales. The single M. urceolatus plant found in Conimbla National Park (NCT) clustered most closely to the M. medius samples. In a PCA of only the ‘northern cluster’ (subset J), M. sp. Silent Grove was distinct (PC1, 21.7%) from M. sp. Inglewood and the northern M. urceolatus specimens (Fig. 8d).

Subset F returned several distinct clusters corresponding with hypothesised species limits (Fig. 8e, f: PC1 = 36.3%; PC2 = 17.7%; PC3 = 11.9%). Melichrus sp. Herberton samples clustered with M. sp. Kroombit Tops, despite a vast geographic disjunction and all M. sp. Gurulmundi samples formed a distinct cluster. The remaining samples in this subset formed four discrete clusters corresponding with the OTU designations, M. adpressus sensu Paterson, M. sp. Boonoo Boonoo, M. sp. Yuraygir and M. sp. Wahlmoorum.

Melichrus adpressus sens. str. clustered with the populations of M. sp. Galambary to which this is geographically closest (NCT) across the first three axes of a PCA of subset I (Fig. 8g, h: PC1 = 15.6%; PC2 = 13.3%; PC3 = 8.6%). Other populations of M. sp. Galambary were spread across the three axes in clusters apparently corresponding with the geographical proximity or environmental continuity. The single M. sp. Gilgandra population included in the analysis clustered near to but always separate from the Central Tablelands M. sp. Galambary and M. adpressus sens. str. clusters.

In subset H, Melichrus gibberagee populations were strongly separated from the M. hirsutus populations across PC1 (42.3%: Fig. 8i). Melichrus sp. Tara samples clustered with the two M. sp. Isla Gorge populations to which this is geographically closest across all major PCA axes (subset G; Fig. 8j). All other Melichrus sp. Isla Gorge samples clustered by population, according to geographic proximity to each other.

For subset A, an ancestry model with eight populations (K = 8) had the highest mean log-likelihood (Supplementary Fig. S5a, b) and the individual log-likelihood values for all 10 runs at K = 8 were similar. The major (5/10) and minor (3/10 and 2/10) modes showed similar patterns of admixture (Supplementary Fig. S6). For K = 9, two runs with outlier values of log-likelihood were removed (Supplementary Fig. S5b), producing a slightly higher mean log-likelihood than that of K = 8. This indicated that an admixture model with eight or nine ancestral populations likely best represents this subset of the data (Fig. 9). Log-likelihood values for runs at a higher value of K (K = 10–13) did not converge. Ancestry models with K = 8 and K = 9 were best supported but across all values of K, consistent patterns of ancestry and admixture emerged (Supplementary Fig. S6).

Fig. 9.

STRUCTURE bar plots for Subset A, a dataset of 2636 DArTseq SNPs. The major clustering modes for the best supported number of ancestral populations, K = 8 and K = 9 are shown. White lines separate individuals. Black lines separate populations. Each OTU is labelled and the botanical district of each population is given (Centre for Australian National Biodiversity Research 2018). Within each OTU, populations are arranged by latitude. Ancestral populations are coloured by the OTU (see Table 1 and Fig. 1) that these primarily but not necessarily exclusively represent. The asterisk (*) indicates M. sp. Wahlmoorum.

Corroborating results observed in the PCA, all OTUs from the Melichrus procumbens cluster shared a common ancestry distinct from other OTUs in the STRUCTURE analysis of subset A (Fig. 9: green). Samples of M. procumbens from the southern extent of the species range drew exclusively from this ancestral population, with populations of M. procumbens further north showing admixture with an additional ancestral population. Melichrus sp. Torrington and M. sp. Silent Grove that co-occur in Torrington National Park (NNT) also drew substantially from this population (Fig. 9: purple).

All samples of Melichrus medius shared a single ancestral population in the STRUCTURE analysis of subset A (Fig. 9: yellow). Melichrus urceolatus from the southern extent of the Northern Tablelands (sNNT) partially shared ancestry with M. medius and some with M. urceolatus (Fig. 9: dark blue). Mirroring results observed in the PCA (Fig. 8c), M. urceolatus samples showed a major divide in genetic clustering corresponding to a north-south geographic division. Melichrus urceolatus from most of New South Wales, ranging from the South Coast through to the North Western Slopes shared most of the ancestry (Fig. 9: dark blue). Melichrus urceolatus from north of this divide – occurring in the northernmost extent of the Northern Tablelands and the North Coast in New South Wales and the Darling Downs in Queensland – derived most ancestry from a separate population that the species shares with M. sp. Inglewood (Fig. 9: light blue). There is a small but notable amount of admixture between these two ancestral populations in samples from the North Western Slopes. The possibility that this pattern of differentiation is a product of isolation by distance was explored through conStruct analysis (see below).

The STRUCTURE analysis of subset A provided strong support for Melichrus sp. Herberton and M. sp. Kroombit Tops as conspecific, distinct from the other M. adpressus sensu Paterson segregates (Fig. 9: bluish-green). Melichrus sp. Gurulmundi was also supported as a distinct entity (Fig. 9: grey), though this shared some ancestry with samples of M. adpressus sensu Paterson from the Moreton district (Fig. 9: pink). At K = 8, M. sp. Boonoo Boonoo and M. sp. Yuraygir also shared this cluster but at K = 9, an additional ancestral population, primarily corresponding with M. sp. Yuraygir but also with M. sp. Boonoo Boonoo is introduced (Fig. 9: brown). As in other samples from the Northern Tablelands, M. sp. Boonoo Boonoo also shared the ancestral population that defines M. sp. Silent Grove (Fig. 9: purple). As expected, the program did not return a distinct grouping for the single, isolated specimen of M. sp. Wahlmoorum.

For a STRUCTURE analysis of subset B, an ancestry model with four ancestral populations (K = 4) had the highest mean log-likelihood (Supplementary Fig. S5c, d) and the individual log-likelihood values for all 10 runs at K = 4 were similar. The major (6/10) and minor (4/10) modes showed similar patterns of admixture in the plots, with one important difference (Fig. 10). The ancestry of the only known population of M. sp. Kaputar was derived from the population’s own exclusive cluster in the major mode, whereas most of the ancestry in the minor mode was shared with M. sp. Galambary. Two runs for K = 5 were removed as these had very different log-likelihood values from the remaining eight runs (Supplementary Fig. S5d). Once the outliers were removed, K = 5 had a slightly higher mean log-likelihood than K = 4, indicating that an admixture model with four or five admixing populations likely best represents this subset of the SNP data.

Fig. 10.

STRUCTURE bar plots for Subset B, a dataset of 2710 DArTseq SNPs. The major and minor clustering modes for the best supported number of ancestral populations, K = 4 (major and minor clustering modes) and K = 5 (major mode: 6/10) are shown. White lines separate individuals. Black lines separate populations. Each OTU is labelled and the botanical district of each population is given (Centre for Australian National Biodiversity Research 2018). Within each OTU, populations are arranged by latitude. Ancestral populations are coloured by the OTU (see Table 1 and Fig. 1) that these primarily but not necessarily exclusively represent.

Melichrus sp. Gardens of Stone was represented primarily by a single exclusive ancestral population in the STRUCTURE analysis of subset B (Fig. 10: red). All samples of M. erubescens shared either one (Fig. 10: brown: K = 4 major mode, K = 5 all modes) or two ancestral populations (Fig. 10: brown and dark blue: K = 4 minor mode). At K = 4, M. sp. Galambary, M. sp. Gilgandra and M. adpressus sens. str. shared a single population (Fig. 10: dark green). At K = 5 an additional population was introduced, from which northern populations generally derived more ancestry (Fig. 10: yellow). The possibility that the pattern of differentiation observed in this subset was a product of isolation by distance was explored through conStruct analysis (see below). STRUCTURE plots for all values of K are included in Supplementary Fig. S7.

For a STRUCTURE analysis of subset C, an ancestry model with four populations (K = 4) had the highest mean log-likelihood (Supplementary Fig. S5e, f) and the individual log-likelihood values for all 10 runs at K = 4 were similar. A model with four admixing ancestral populations likely best represents this subset of the SNP data. For the major mode (9/10), all populations of M. sp. Mareeba showed approximately equal admixture from ancestral populations corresponding to M. sp. Isla Gorge and M. sp. Herberton + M. sp. Kroombit Tops respectively (Fig. 11: olive-green and bluish-green). This pattern, characteristic of hybridisation, was evident across most values of K (Supplementary Fig. S8). At all levels of parametrisation in this analysis, populations of M. hirsutus and M. gibberagee shared a single ancestral population (Fig. 11: blue). A STRUCTURE run consisting of only these two OTUs (subset H) returned a result where K = 2 was decisively the best supported number of populations (Supplementary Fig. S5g) and the two ancestral populations corresponded with each species, with negligible admixture detected (Fig. 11: blue and dark green; all plots included in Supplementary Fig. S9).

Fig. 11.

STRUCTURE bar plots for subset C and subset H, datasets of 2915 and 1568 DArTseq SNPs respectively. Left: subset C, best supported number of ancestral populations, K = 4 major clustering mode. Right: subset H, best supported number of ancestral populations, K = 2, major clustering mode. White lines separate individuals. Black lines separate populations. Each OTU is labelled and the botanical district of each population is given (Centre for Australian National Biodiversity Research 2018). Within each OTU, populations are arranged by latitude. Ancestral populations are coloured by the OTU (see Table 1 and Fig. 1) that these primarily but not necessarily exclusively represent. The asterisk (*) indicates M. sp. Wahlmoorum.

To test how isolation by distance may have contributed to the formation of the northern Melichrus urceolatus genetic group observable in the PCA and STRUCTURE results, subset D (see Table 3) was analysed with conStruct. The spatial model had a higher predictive accuracy overall than the non-spatial model (Supplementary Fig. S10) and the spatial model with six layers (K = 6) had the highest mean predictive accuracy of any model. For the non-spatial model, K = 8 had the highest mean predictive accuracy. The layer contribution analysis showed that for the spatial model only the first five layers contributed significantly to the total explained variation (Supplementary Fig. S11). As models with four or five layers had similar mean predictive accuracy, these were approximately equally suitable models for the dataset (K = 5, spatial (a) and non-spatial (b) are shown in Fig. 12).

Fig. 12.

conStruct analysis results for subsets D and B of the DArTseq SNP dataset (see Table 3). (a) Subset D; K = 5 spatial model. (b) Subset D; K = 5 non-spatial model. (c) Subset B; K = 4 spatial model. (d) Subset B; K = 5 spatial model. Each pie chart represents a single population. For figures (a) and (b) all unlabelled populations are M. urceolatus. For figures (c) and (d) the populations that are not labelled are either M. sp. Galambary and M. adpressus sens. str. (primarily drawn from the dark green layer) or M. erubescens (primarily drawn from the brown layer). The scale represents geographic distance. Layers are coloured by the OTU (see Table 1 and Fig. 1) that these primarily but not necessarily exclusively represent.

Across several values of K, including the best supported (K = 5, Fig. 12a), the spatial model delivered a result with a clear distinction between the northern and southern M. urceolatus groups. This pattern was also evident in the non-spatial model but was less clear in the noise of several other layers describing differentiation across the geographic range. This indicates that the ancestry pattern observed in the STRUCTURE and PC analyses is likely substantially independent of the influence of isolation by distance. Both the spatial and non-spatial models found that all M. medius populations drew nearly exclusively from a single layer, across several levels of parameterisation (including K = 5). The spatial model showed a stronger relationship between the M. medius and M. urceolatus populations from the southern extent of the Northern Tablelands than was found by the non-spatial model.

Subset B was subjected to a conStruct analysis to test if the distinction between Melichrus sp. Kaputar, M. sp. Galambary and M. sp. Gardens of Stone observed in the PC and STRUCTURE analyses could be attributed to isolation by distance. The spatial model had a higher predictive accuracy overall than the non-spatial model (Supplementary Fig. S12) and the spatial model with six layers (K = 6) had the highest mean predictive accuracy of any model. For the non-spatial model, K = 9 had the highest mean predictive accuracy. The layer contribution analysis showed that for the spatial model only the first four or five layers contributed significantly to the total explained variation (Supplementary Fig. S13). The model with five layers had a higher mean predictive accuracy than K = 4 but contributed less in the layer contribution analysis, therefore on balance spatial models with K = 4 or 5 are likely the most suitable models for the dataset (Fig. 12c, d).

Across several values of K, including the best supported, the layer composition of Melichrus sp. Kaputar differed only slightly from that of M. sp. Galambary populations across both the spatial and non-spatial analyses (Fig. 12c, d). Melichrus sp. Gardens of Stone, however, while showing some amount of admixture with M. sp. Galambary, drew substantially from an independent layer (Fig. 12c, d, shown in red).

Filtration of the dataset for population statistics produced a set of 4936 high-quality SNPs and retained 524 individuals across 95 populations. Observed heterozygosity (H_O) was generally slightly lower than expected heterozygosity (H_E) and F_IS was moderately positive or negative for most populations (Fig. 13). The notable exception was Melichrus sp. Mareeba, for which all populations had H_O values approximately twice the magnitude of the H_E values and strongly negative F_IS values (Supplementary Table S4).

Fig. 13.

Observed heterozygosity (H_O), expected heterozygosity (H_E) and inbreeding coefficient (F_IS) for all populations of Melichrus with three or more individuals for a DArTseq dataset of 4936 SNPs. Each tick represents one population. Populations are labelled by OTU (Table 1).

Pairwise F_ST values estimated for all pairs of populations in the study ranged from 0.01–0.93 (Supplementary Table S5). Intraspecific F_ST values for Melichrus sp. Mareeba populations (0.02–0.06) were an order of magnitude lower than interspecific F_ST values between M. sp. Mareeba and all other populations (0.43–0.75). Comparison of F_ST values between populations was consistent with the pattern of diversity observed in the PCA and STRUCTURE results.