Reticulate Evolution, Ancient Chloroplast Haplotypes, and Rapid Radiation of the Australian Plant Genus Adenanthos (Proteaceae)

Introduction

Hybridization is important in the evolution of many plant groups (Arnold, 1992; Soltis and Soltis, 2009; Givnish, 2010). Examples of gene flow between species are common in plants across many different evolutionary and phylogenetic scales, from deep reticulate introgression events (Folk et al., 2017; García et al., 2017) to the evolutionary process of speciation by hybridization (Mallet, 2007; Soltis and Soltis, 2009). Reticulation can be indicated by discordance between organellar (plastid and mitochondrial) and nuclear molecular datasets, due to the different modes of inheritance and evolution between the two genomes (Birky, 1995; Soltis and Kuzoff, 1995; Small et al., 2004). However, cytonuclear discordance may also result from incomplete lineage sorting (ILS) or poor resolution among sampled loci (Willyard et al., 2009; Gurushidze et al., 2010). Addressing the causes of cytonuclear incongruence is increasingly realistic using next-generation sequencing (NGS) approaches including targeted hybrid capture (Lemmon et al., 2012; Lemmon and Lemmon, 2013; Weitemier et al., 2014). These methods, which can generate sequences from multiple nuclear and organellar loci, allow rigorous exploration of causes of cytonuclear incongruence, including hybridization, using robustly supported phylogenies (Howarth and Baum, 2005; Vargas et al., 2017).

Along with the developments in sequencing technology, there has been significant progress in analytical approaches to untangling the influence of hybridization and ILS on cytonuclear discordance. While studies applying these approaches cannot rule out the presence of ILS, they can confidently separate the signals of hybridization from ILS (e.g., Joly et al., 2009). Several recent studies have made inferences of deep reticulation from multiple introgression events throughout the evolutionary history of their study groups (Folk et al., 2017; García et al., 2017). Introgression can play an important role in plant evolution and has been linked to rapid radiations in some of these groups (Seehausen, 2004). Most studies to date have focused on Northern Hemisphere plants (Francisco-Ortega et al., 1996; Barrier et al., 1999; Stankowski and Streisfeld, 2015). Different evolutionary drivers may have been involved in the Southern Hemisphere due to the older age of its biota (Hopper, 2009 and references therein). Many prominent lineages in the Australian contemporary flora are thought to have originated in the Cretaceous (Crisp et al., 2011; Lamont and He, 2012; Crisp and Cook, 2013) and show a radiation pulse in the mid-Cenozoic (25–10 Ma) (Crisp et al., 2004) in response to increased seasonality initiated at the end of the Eocene (c. 33 Ma) and subsequent aridification after the mid-Miocene (c. 14 Ma) (Macphail, 2007). However, no studies to date have explored the link between large scale climatic change, radiation, and hybridization in the early evolution of Australian plants. Similarly, while adaptive introgression has been shown to have spurred radiations of many groups in other regions of the world (Barrier et al., 1999; Seehausen, 2004; Givnish, 2010), a conclusive link between the two has not yet been demonstrated in Australia.

Natural hybridization has been documented in a number of Australian plants (Ashton and Sandiford, 1988; Griffin et al., 1988; Sedgley et al., 1992; Holman and Playford, 2000; Walker et al., 2009). In a few cases it is extensive (Leach and Whiffin, 1978; Potts and Reid, 1985; McIntosh et al., 2014). Hybridization has been suspected in the endemic Australian plant genus Adenanthos Labill. (Proteaceae), based on morphology alone (Nelson, 1977), and has subsequently been confirmed using molecular data (Walker et al., 2018). Adenanthos comprises 31 extant species, the majority of which (29 of 31 species) occur in southwest Western Australia (SWA) (Figure 1). Two species are disjunct from the rest of the genus across the Nullarbor Plain and are restricted to the southern peninsulas of South Australia. The genus consists of perennial shrubs or in some cases trees and are thought to be bird pollinated (Keighery, 1982; Collins and Rebelo, 1987). High outcrossing rates associated with bird pollination and also general self-incompatibility found in most members of Proteaceae suggests that population dynamics of Adenanthos might be fundamentally different to the majority of plants that are insect pollinated (Keighery, 1982; Goldingay and Carthew, 1998). However, detailed population genetic studies on Adenanthos are currently lacking. All species of Adenanthos and its close relatives (e.g., Isopogon Knight, Petrophile R.Br. ex Knight, Leucadendron L.) that have chromosome counts are n = 13 (Ramsay, 1963; Stace et al., 1998); ploidy variation within genera and clades is relatively uncommon in Proteaceae.

Figure 1. Representative floral and leaf diversity of Adenanthos: (A) Adenanthos cygnorum subsp. chaemaephyton E.C.Nelson. (B) A. stictus A.S.George. (C) A. glabrescens subsp. exasperatus E.C.Nelson. (D) A. glabrescens E.C.Nelson subsp. glabrescens. (E) A. linearis Meisn. (F) A. venosus Meisn. (G) A. × cunninghamii Meisn. (H) A. obovatus Labill. (I) A. forrestii F. Muell. (J) A. sericeus [Labill. cultivated.] (K) A. macropodianus [E.C.Nelson.] (L) A. terminalis [R.Br.] Photos: F. J. Nge.

Adenanthos diverged from its sister-group (Leucadendrinae P. H. Weston and N. P. Barker) at the Eocene–Oligocene boundary (stem age c. 33.9 Ma) (Sauquet et al., 2009), thus providing an excellent case study to investigate potential reticulate patterns of evolution across deep timescales in the context of the Australian flora. Here, we use a NGS hybrid capture approach to infer nuclear and chloroplast phylogenies of Adenanthos to: (1) reconstruct its evolutionary and biogeographic history, and (2) assess for signs of hybridization and deep reticulate evolution within the genus.

Materials and Methods

Sampling

We included 44 samples (30 of the 31 recognized species and 2 putative hybrids) covering all infrageneric sections and subsections within Adenanthos according to the most recent taxonomic revision (Nelson, 1977). We included a natural hybrid between A. cuneatus and A. sericeus (A. × cunninghamii) in our study. Half of the samples were collected in the field with fresh leaf tissue dried in silica gel (Supplementary Table S1A). The remaining samples were sourced from recently collected herbarium specimens (after 1960) lodged in PERTH and AD (Supplementary Table S1B).

DNA Extraction, Library Preparation, and Sequencing

Approximately 20 mg of silica dried leaf material per sample was used for DNA extractions, performed by Intertek Group plc using magnetic bead-based chemistry. We used a set of 30–100 single-copy nuclear and 13 plastid loci developed as phylogenetic markers for angiosperms (Waycott et al., in preparation) using the MYBaits target enrichment system (MYcroarray, Ann Arbor, Michigan) for sequence capture of the selected loci. In brief, genomic DNA (normalized to 1 ng/uL) was sheared using a Diagenode Bioruptor Pico sonicator for seven cycles to fragment lengths of c. 400–600 bp. DNA libraries were constructed using a JetSeq Flex DNA Library preparation kit (Bioline). To enable bioinformatics processing following hybrid capture, two 8 bp synthetic barcodes were annealed at each end of the DNA fragments. During the hybrid capture protocol, for each 96-well plate, the first barcode is replaced every 48 samples (i.e., two barcodes, one for each half of the plate), while the second barcode is unique to each sample of each half-plate (i.e., 48 different barcodes). This ensured that each sample has a unique combination of the two barcodes for downstream identification. Libraries were pooled in equimolar concentrations and sent for Illumina paired-end sequencing (2 × 150) on a lane of a HiSeqX Ten at the Garvan Institute for Medical Research in Sydney.

Sequence Assembly

High-throughput 150 bp paired-end reads were processed using CLC Genomics Workbench v7.5.1¹. Following demultiplexing and quality trimming (Phred-score threshold of 20), we used de novo assembly of pooled Adenanthos samples to generate a set of reference contigs for each sample. In order to recover the targeted nuclear loci, the de novo assembly was converted to a BLAST database and reference genomic sequences of Aquilegia coerulea (downloaded from Phytozome v 12²) used as query sequences using an E-value ≤ 1E-20. The de novo contigs matching the Aquilegia genes were used as a mapping reference for each individual to generate a per sample assembly at each locus. From these, we extracted the majority rule consensus sequence inserting “Ns” when coverage was lower than 5.

The resultant mapping files were exported in BAM format and allele phasing was performed using SAMTools Phase with default parameters applied (Li et al., 2009). SAMTools calls heterozygous SNPs (single nucleotide polymorphisms) at one site and segregates the reads (which contain one or the other heterozygous SNP) into two new “phased” BAM files. Reads lacking the given SNP site (but in part overlapping the segregated reads) are segregated randomly to either BAM file. The phased BAM files were then imported into Geneious v.1.11.5 (Kearse et al., 2012), and a majority-rule consensus extracted using a 65% cut-off, then aligned using the MUSCLE (Edgar, 2004) plugin with default parameters.

In the majority of samples, we also recovered the 18S–26S nuclear ribosomal internal transcribed spacer (ITS) region. Although not specifically targeted, nuclear ribosomal DNA has a high copy number and can be recovered as by-catch. We used an ITS reference sequence (Isopogon sphaerocephalus, GenBank accession number AF508820.1) as a query sequence for BLAST and generated a per sample assembly as outlined above.

Plastid (chloroplast) targets were recovered using the chloroplast genome sequence of Macadamia integrifolia (GenBank reference number 34480) as a mapping reference. Reads from each sample were mapped to the reference using default parameters with a length fraction of 0.7 and a similarity fraction of 0.9. Consensus sequences were extracted as above. Consensus sequences for each individual and locus were imported into Geneious v1.11.5 (Kearse et al., 2012) and aligned using the MUSCLE (Edgar, 2004) plugin with default parameters, then manually checked and adjusted. Samples with more than 70% missing data in both nuclear and plastid alignments were excluded from our final dataset.

Phylogenetic Analyses and Divergence Time Estimation

Maximum Likelihood (ML) analyses were implemented in RAxML v.8.2.10 (Stamatakis, 2014) for two datasets: (1) 35 nuclear contigs phased and unphased (44 taxa, 25,646 bp), and (2) concatenated chloroplast sequences (43 taxa, 34,218 bp), using the GTR + I + G substitution model and bootstrap support obtained with 1,000 standard boostrap replicates. Single-gene trees were also estimated for a subset of our unphased nuclear dataset (19 nuclear contigs) that excluded potential paralogs, with 100 standard bootstrap replicates each. We used BLAST searches against de novo assemblies to screen for potential paralogs, assuming that divergent and overlapping contigs recovered for a single target gene represent paralogy. To assess for phylogenetic congruence and signal among loci, well-supported clades (>75% bootstrap) in each nuclear gene tree were compared with all other gene tree topologies manually. Bayesian analyses were conducted in BEAST v.2.4.7 (Bouckaert et al., 2014) for our concatenated datasets to obtain age estimates for Adenanthos for each dataset using a range of fossil calibration regimes (see Supplementary Tables S2–S4 for details). Available nuclear (ITS) and plastid (matK, rbcL) sequences for outgroups were sourced from GenBank (Supplementary Table S5). For these three gene regions, we used the fossil calibrations applied in the Proteaceae family-wide study of Sauquet et al. (2009) to obtain divergence estimates for Adenanthos. One fossil calibration point (Cranwellipollis palisadus; for stem of Franklandia) was available within subfamily Proteoideae, which includes Adenanthos. We also included five additional calibration points in other subfamilies within Proteaceae to increase the accuracy of these estimates, applying uniform calibration priors following recommended practice (Sauquet et al., 2012). Because NGS sequences were not available for the outgroups, we applied similar calibration regimes for our full NGS datasets and compared the divergence age estimates with those obtained from secondary calibrations derived from these estimates. Secondary calibrations included (i) the stem age of Isopogon, the sister genus of Adenanthos and Leucadendrinae (set as log-normal distribution, offset = 41.5 Ma, SD = 0.23), and (ii) Adenanthos and Leucadendrinae crown (set as log-normal distribution, offset = 16 Ma, SD = 0.23) obtained from the plastid (matK) BEAST run.

BEAUti v2.4.7 was used to create input files for BEAST. We used a GTR + I + G substitution model and a relaxed lognormal clock model. Three parallel BEAST runs were performed for each analysis with the number of Markov Chain Monte Carlo (MCMC) generations and sampling frequency dependent on the size of the dataset (Supplementary Table S6). The first 20% of runs were discarded as burn-in. Tracer v1.6.0 (Rambaut et al., 2015) was used to assess convergence of the posterior, which was determined when effective sample size (ESS) reached ≥200. Tree output files were combined using LogCombiner v2.4.7, summarized in TreeAnnotator v2.4.7, and visualized using FigTree v1.4.3 (Rambaut, 2012). Lineage-through-time plots were constructed from pruned nuclear and chloroplast BEAST trees where each species was represented with only one terminal, using the ‘phytools’ package (Revell, 2012) in R (R Core Team, 2016).

All RAxML and BEAST analyses were run on the CIPRES Science Gateway portal (Miller et al., 2010). Conflicts between the nuclear and chloroplast ML topologies were visualized using the tanglegram tool in Dendroscope v. 3.5.10 (Scornavacca et al., 2011; Huson and Scornavacca, 2012).

Hybridization Assessment

In order to distinguish nuclear and chloroplast topological discordance as a result of either hybridization or ILS, we simulated plastid gene trees using scaled nuclear trees to obtain estimates of ILS for the plastid dataset, following the approach of García et al. (2017) and Folk et al. (2017). The simulated trees were then compared with the actual plastid topology and hybridization events inferred by areas that are incongruent. We utilized ASTRAL v. 5.6.3 (Zhang et al., 2018) to obtain a species tree with coalescent branch lengths from the individual nuclear gene trees obtained through RAxML. ASTRAL utilizes unrooted gene trees to generate phylogenetic quartets, which is relevant for our dataset as our individual gene trees did not contain outgroups. This is beneficial as random rooting can mimic the coalescent process (Rosenfeld et al., 2012; Tian and Kubatko, 2014). Support was assessed through the final normalized quartet scores of the overall species trees and local posterior probabilities of each branch terminal as a measure of gene tree conflict. Branch lengths of the species tree were scaled by a factor of four to account for organellar inheritance, as maternal inheritance of the plastid genome is typical for flowering plants (Mogensen, 1996). We simulated 1,000 gene trees from the scaled ASTRAL species tree by applying a coalescent model using a python-based script with DendroPy (Mirarab et al., 2014). The simulated gene trees were visualized in DensiTree 2 (Bouckaert and Heled, 2014).

We used JML v.1.3.1 (Joly, 2012) following the approach of Joly et al. (2009) to assess the relative contributions of hybridization and ILS to discordance. This method uses the posterior distribution of species trees, population size, and branch lengths estimated in ^∗BEAST (Bouckaert et al., 2014) to simulate sequence data under coalescent scenarios with no migration. To achieve this the minimum pairwise distance between sequences of two extant species from the simulated dataset were compared with empirical data. Hybridization or introgression can be inferred when observed pairwise distances from empirical data are significantly smaller than the simulated dataset derived from JML analyses, rejecting ILS as the only cause for topological conflict between the datasets (Joly et al., 2009; Joly, 2012). A coalescent tree was inferred in ^∗BEAST using a combined nuclear and chloroplast dataset. Genetic distances were calculated using JML with only the chloroplast dataset, as Joly et al.’s (2009) approach assumes that the markers used to estimate the genetic distances are non-recombinant. JML analyses were conducted on the complete dataset, and also on individual clades which were shown to be incongruent between our nuclear and chloroplast topologies separately (Supplementary Table S7). We tested the performance of JML to assess deeper introgression events, evident by particular subclades showing incongruence among nuclear and chloroplast datasets. Analysis of subsets of the total data were conducted where only a fraction of the sampled taxa were included, including only one representative from each pair of conflicting clades across the two topologies. For each analysis, 1,000 simulations were computed for the chloroplast pairwise distance comparisons.

Haplotype and Splitstree Networks

As bifurcating trees may not accurately represent reticulate events among closely related taxa, we used network analyses to better represent relationships and assess for conflict between the nuclear and chloroplast datasets. Haplotype networks were constructed for the chloroplast dataset using TCS 1.13 (Clement et al., 2000) in PopART v.1.7 (Leigh and Bryant, 2015), classified into different subregions in SWA according to the Interim Biogeographic Regionalization for Australia (IBRA7) bioregional classification scheme³. Distance-based Neighbour-Nets were created in SplitsTree v.4.14.4 (Huson and Bryant, 2006) for (i) nuclear, (ii) plastid, and (iii) combined datasets using uncorrected p-distances.

Results

Our plastid alignment through reference mapping and BLAST for downstream analyses contained 13 curated contigs that were 34,218 bp in length and included 43 taxa. Total curated nuclear alignments had a length of 25,646 bp comprising 35 independent loci after potential paralogs were excluded, covering 44 sampled taxa. Of the 44 samples, 39 were present in both nuclear and plastid datasets after removal of taxa with poor quality sequences or missing data.

Phylogenetic Relationships and Conflict in Nuclear and Plastid Data

Incongruence is significant between the plastid and nuclear ML topologies, with the two datasets recovering a different number of clades and statistically well-supported conflicting relationships across species and clades (Figure 2) (for further details see Supplementary Results).

Figure 2. Tanglegram comparing the phylogenetic placement of Adenanthos taxa connected by gray lines from maximum-likelihood concatenated RAxML nuclear (left) and plastid (right) trees. Bootstrap support (BS) is indicated above branches: very strongly supported clades (BS = 100%) are represented by “*”, weakly supported branches with <50% BS are represented by “–”. Letters (A–F) and (A’–E’) refer to clades discussed in text. Terminal tips represent taxa from A. sect. Eurylaema (red circles), A. sect. Adenanthos subsect. Anaclastos (green squares), and A. subsect. Adenanthos (unlabeled).

Both the concatenated ML and coalescent ASTRAL analyses gave largely congruent results for the nuclear dataset (Figures 2, 3). Several clades are resolved with high support from the ML topology (Clade C: bootstrap BS = 92; Clade D: BS = 82; Clade E: BS = 90; Clade F: BS = 74), however, the backbone of the tree was unresolved (Figure 2). Similarly, the backbone of the plastid topology was largely unresolved.

Figure 3. ASTRAL species trees of Adenanthos inferred using the multispecies coalescent model and gene trees of nuclear loci. Support values for the ASTRAL tree represent the local quartet score of each branch. Letters (A–D) refer to clades discussed in text. Terminal tips represent taxa from A. sect. Eurylaema (red circles), A. sect. Adenanthos subsect. Anaclastos (green squares), and A. subsect. Adenanthos (unlabeled).

Divergence Age Estimates and Radiation of Adenanthos

Age estimates for Adenanthos obtained from our NGS nuclear trees were older than those from the plastid dataset, despite employing similar fossil calibration constraints (Supplementary Table S8). The divergence time estimates from the combined ITS, matK, and rbcL, as well as ITS-only topologies, are closer to those obtained from our NGS plastid topology (Supplementary Table S8). These differences in divergence time estimates were consistent across all fossil calibration schemes, including those obtained from secondary calibration of NGS data only, when outgroups were excluded due to missing data (Supplementary Table S8). We focus on the divergence age estimates of the NGS plastid and ITS combined topologies here, as the older age estimates from the nuclear NGS data likely reflect the lack of available NGS data for outgroup taxa used for the calibration regimes as well as missing data from our NGS nuclear dataset.

The stem age of Adenanthos was estimated at 36.1 Ma (95% CI: 15.3–33.2 Ma) for our NGS plastid topology and 38.4 Ma (95% CI: 15.6–36.6) for the ITS, matK and rbcL combined topology, employing the Proteaceae-wide calibration scheme (Figure 4A and Supplementary Table S8). The crown age was estimated at 24 Ma (95% CI: 15.3–33.2 Ma) for the plastid topology and 25.1 Ma (95% CI: 24.5–41.0 Ma) for the combined ITS and plastid topologies, respectively, in the late-Oligocene–Miocene (Figure 4A and Supplementary Table S1). Radiations of clades in the mid-Miocene (15–20 Ma) was consistent in both the NGS plastid and ITS-plastid combined topologies (Figure 4A and Supplementary Figure S1).

Figure 4. (A) Calibrated BEAST chronogram of Adenanthos inferred from 13 chloroplast markers, with branch lengths scaled according to time (Ma). Non-Proteoideae outgroups used in the family-wide calibration scheme were pruned. Divergence times (Ma) are given at each node, with 95% confidence intervals indicated as blue bars. Clade labels denote the clades from the plastid ML RAxML tree. Adenanthos stictus illustration denotes the crown of the genus. Star denotes the fossil calibration used in Proteoideae (Franklandia). Subspecies are abbreviated by the first two letters of the epithet. (B) Lineage-through-time plot of Adenanthos based on the BEAST plastid topology. (C) Lineage-through-time plot of Adenanthos based on the BEAST ITS nuclear topology.

Interestingly, both our chloroplast and nuclear trees (ITS and NGS) showed that the southeastern Adenanthos clade is strongly nested within other SWA clades, even despite the strong topological conflicts between the two datasets. The divergence of the southeastern Australian clade from one of the SWA subclades was estimated at 16.5 Ma (95% CI: 9.5–24.8 Ma) based on the plastid NGS topology (Figure 4A). In contrast, the crown age of the southeastern species was inferred to be significantly older based on the ITS topology, with divergence of A. macropodianus from A. terminalis estimated in the Miocene c. 8.4 Ma (95% CI: 2.8–14.7 Ma) compared with a Pleistocene divergence c. 1.3 Ma (95% CI: 0.29–2.9 Ma) inferred from the NGS plastid topology (Figure 4A and Supplementary Figure S2).

Hybridization Assessment

Reticulate evolution and introgression within Adenanthos were detected amongst our molecular datasets, indicated by the widespread discordance between the nuclear and plastid topologies. Evidence for hybridization was supported by both the gene tree simulations and JML approaches. The simulated plastid gene tree distribution derived from the nuclear ASTRAL species tree indicates that the discordance between the nuclear and plastid datasets is at least partly due to hybridization, as the simulated plastid topology did not match the actual plastid topology (Figure 5). The discordance between the simulated and actual plastid topology indicates that several of the observed plastid-nuclear discordances are almost never expected to occur under coalescence alone, indicating that they are unlikely to be caused by ILS alone.

Figure 5. Severe incongruence between observed empirical plastid topology (red, marked with a star) and simulated plastid consensus tree (blue) with incomplete lineage sorting but no introgression – simulated from ASTRAL nuclear species tree, with a branch scaling factor of 4. The distribution of simulated plastid trees are indicated by thin green lines drawn in DensiTree.

Instances of ancient and putative recent hybridization were supported by our JML analyses (Supplementary Table S9). In the dataset with 25 taxa, only 9 out of 253 pairwise comparisons had non-significant values (p value > 0.1); that is, 96% (244/253) of the pairwise chloroplast distance comparisons are significantly smaller than expected in a scenario with only ILS (p value < 0.1) (Supplementary Table S9A). This indicates that the model cannot accurately predict the observed minimum distances and that a strict bifurcating species tree model is inadequate, due to the presence of hybridization. Subsequent analyses on the three subsets all detected signals of introgression between different clades within Adenanthos (Supplementary Table S9B).

Haplotype and Reticulate Networks

The majority of species have chloroplast haplotypes that were exclusive to each taxon, with unique haplotypes for sampled individuals within species (Figure 6). Only four instances where haplotypes were shared across multiple species were detected in our genus-wide haplotype network. Twenty haplotypes inferred by TCS were not present in analyzed individuals, and represent missing individuals or extinct lineages in the network. Interspecific geographic structure was evident in the network, with instances of shared haplotypes in species occurring in the same region (Figure 6 and Supplementary Figure S3). In contrast, less infraspecific geographic patterning was noted with different populations of species exhibiting unique haplotypes across their geographic range, and in some instances scattered across the network. The two southeastern Australian species share a single haplotype, as do sympatric populations of A. ellipticus and A. cuneatus in the South Coast region of SWA. Southeastern Australia has less chloroplast haplotype diversity than SWA, containing only one haplotype shared between its two species. This haplotype is also less divergent than other haplotypes in SWA, being separated from its nearest extant haplotype by two extinct haplotype lineages. Potentially long-persisting haplotypes were detected across multiple species; in some cases the age of the chloroplast haplotype pre-dates the radiation of the lineage i.e., these haplotypes persisted in extant lineages from their most recent common ancestor (Table 1, Figure 4, and Supplementary Figure S2).

Figure 6. (A) TCS network of Adenanthos chloroplast haplotypes colored according to subregions of southwestern and South Australia. Circle size represents number of taxa (1–3) sharing chloroplast haplotypes. Black dots represent missing, unsampled, or extinct intermediate haplotype lineages. (B) Map of Australia with box highlighting southwestern Australia and colored regions in South Australia indicate distribution of Adenanthos terminalis and macropodianus. (C) Map of southwestern Australia showing each of the subregions.

Table 1. Divergence age estimates for long-persisting chloroplast haplotypes in Adenanthos, with nuclear divergence estimate obtained from the ITS topology.

The splitstree networks suggest that reticulate relationships in both nuclear and plastid data sets are largely confined to the backbone (Supplementary Figures S4–S6). Both nuclear and plastid networks resolved distinct clades which were in conflict, resulting in the combined dataset having high levels of reticulate relationships throughout the network (Supplementary Figure S6).

Discussion

We present the first well-resolved, densely sampled phylogeny of Adenanthos. Our results indicate that extensive hybridization is present throughout the evolution of this genus, including deep reticulation events coinciding with the radiation of the genus in the Miocene. A revised infrageneric classification is also warranted to better reflect evolutionary relationships within the genus.

Phylogenetic Incongruence and Reticulate Evolution

Extensive cytonuclear discordance in Adenanthos is best explained by multiple introgression events throughout its evolution that are particularly evident across deeper timescales. Several instances of recent introgression between closely related extant species were also detected. While we cannot rule out the presence of ILS as a factor contributing to the observed incongruence in our plastid versus nuclear data, our analyses indicate that these conflicts cannot be solely due to ILS, as they are never expected to occur under the coalescent model alone.

At least four independent ancient introgression events in Adenanthos were detected in our plastid simulation and JML analyses. In theory, the JML method was developed for detecting hybrids between extant species pairs (Joly et al., 2009; Joly, 2012), and hence might not be optimal for detecting hybridization events that are ancestral in a clade (i.e., from non-extant lineages) (García et al., 2017; Vargas et al., 2017). However, in our study, this method has been able to detect introgression events between clades, by including a representative subset from each clade in the analyses instead of including all taxa. Indeed, this approach has been recommended to increase statistical power for detecting introgression (Joly, 2012). The detection of ancestral hybridization through JML may also be dependent on the strength of the signal. In Adenanthos, hybridization is widespread and deep reticulation events were detectable using pairwise comparisons of extant taxa, whereas in other groups such as rain-lilies (Hippeastreae; Amaryllidaceae), the signal was insufficient or no longer present across sampled taxa (García et al., 2017). Nevertheless, the presence of introgression across clades in Adenanthos and rain-lilies was corroborated with other gene simulation analyses. We encourage the use of multiple methods for detecting hybridization events throughout the evolution of a study group. We also acknowledge potential limitations in detecting introgression through the use of nuclear data in simulating a plastid topology with only ILS, as nuclear genes may also show signals of introgression. This is not the case in our study, where the majority of species showed no signs of introgression in the nuclear topology except for A. × cunninghamii. Further studies with additional population sampling and methological advancements in distinguishing ILS from introgression should provide us with a greater understanding on this topic.

Different plant groups are prone to different degrees of hybridization (Whitney et al., 2010). While hybrids have been identified in some genera of Proteaceae (Lamont et al., 2003; Pharmawati et al., 2004; Milner et al., 2012; McIntosh et al., 2014; Mitchell and Holsinger, 2018), hybridization appears to be uncommon within the family as a whole. Apart from Adenanthos, only the eastern Australian genus Lomatia shows extensive signals of interspecific hybridization (McIntosh et al., 2014). To our knowledge, our study is the first to demonstrate deep reticulation events within a genus in Proteaceae. The biological mechanisms that maintain species boundaries in the face of such extensive past and present introgression in Adenanthos are currently unknown.

Radiation of Adenanthos and Long-Persisting Haplotypes

The unresolved backbone of Adenanthos found in both our NGS nuclear and chloroplast topologies is suggestive of a rapid radiation in the Oligocene-Miocene. Organismal groups that have undergone a rapid radiation appear to be particularly prone to reticulate evolution (Anderson and Stebbins, 1954; Seehausen, 2004, 2013; Mallet et al., 2007; Genner and Turner, 2011; Vargas et al., 2017). Species boundaries may be more porous during a radiation event (Dilley et al., 2000; Smith et al., 2008), potentially driving reticulate evolutionary patterns and incongruence between nuclear and plastid genealogies, as observed in Adenanthos.

We show that, in some cases, sister species or sampled individuals within species have highly divergent chloroplast haplotypes, each of which is most similar to haplotypes found in phylogenetically distant extant species according to the nuclear topology. For example, the haplotypes of A. obovatus and A. barbigera, which are sister species in the nuclear phylogeny, are separated by at least six extinct haplotypes across the backbone of the network. We interpret this as resulting from introgression between these species and extinct lineages. Introgression may have occurred at any time between the divergence of the sister pair at c. 11.6 Ma and the divergence of populations within A. obovatus and A. barbigera (c. 2 Ma) (Supplementary Figure S2: ITS BEAST topology). The older estimate of c. 11.6 Ma for these potentially long-persisting haplotypes is significantly older than ancient chloroplast haplotypes noted in other plants, for example, c. 4 Ma in Jakob and Blattner (2006).

Complex geographic patterning of Adenanthos chloroplast haplotypes is seen in many areas of southwestern Western Australia, where multiple distinct and highly divergent haplotypes are present in localized areas. The South Coast subregion contains the highest haplotype diversity across SWA, followed by the Esperance, Perth, and Northern Sandplains subregions, indicating that highly divergent, old chloroplast haplotypes have persisted in these areas (Byrne, 2008). Several of these regions (Perth and Northern Sandplains) are also centres of floristic species richness (Hopper and Gioia, 2004) and phylogenetic diversity (Rosauer et al., 2009) in SWA, and hence are of high conservation value. Sniderman et al. (2013) and Nge et al. (2020) have hypothesized that relatively low extinction rates in SWA compared to other regions of Australia, due to climatic buffering over the course of multiple large-scale Eocene–Pleistocene climatic events, is one of the main drivers for these patterns.

The lower chloroplast haplotype diversity in southeastern Australia can be attributed to either higher local extinction rates compared with SWA, or a founder effect where a lineage was dispersed from SWA to southeastern Australia. Our divergence estimates based on the ITS topology for the disjunction of Adenanthos species across southern Australia postdates or coincides with the uplift of the Nullarbor Plain c. 14–13 Ma, which is a strong climatic and edaphic barrier for plant migration between the two southern temperate mesic regions. However, we caution against attributing this divergence solely to this vicariance event, as the 95% CI of this divergence event in Adenanthos (4.8–15.3 Ma) is too wide to discriminate between divergence as a direct result of the uplift of the Nullarbor Plain or post-uplift dispersal (see also Crisp and Cook, 2007). Further studies on the population genetics of the southeastern Australian species and additional sampling of outgroups with NGS nuclear data should provide us with more precise divergence age estimates of the clade in relation to its SWA sister groups. Not only does the southeastern Australian clade contain lower haplotype diversity, the clade is also relatively depauperate, containing only two species compared to more species-rich clades found in SWA. Future research on drivers of this disparity in species richness linking it with genetic mechanisms and results of this study would be especially promising.

The radiation of Adenanthos coincides with that inferred in many other Australian plant radiations during the Miocene (Crisp et al., 2004; Cardillo and Pratt, 2013; Puente-Lelièvre et al., 2013; Jabaily et al., 2014; Mast et al., 2015; Thornhill et al., 2019). Intensification of aridity and seasonality of rainfall across the continent at this time resulted in the retreat of mesic vegetation types and expansion of sclerophyllous and xeromorphic vegetation (Byrne et al., 2011; Crisp and Cook, 2013). These changes opened new niches, potentially allowing Adenanthos and other sclerophyllous groups to diversify. Hybridization between distinct lineages spurring adaptive radiations has been demonstrated for Hawaiian silverwords (Barrier et al., 1999) and African lake cichlids (Meier et al., 2017), with introgression providing genetic novelty leading to diversification into new niches. A recent review (Berner and Salzburger, 2015) has suggested that many adaptive radiations exhibit signals of hybridization, highlighting an important link between novel genetic variation derived from hybridization or introgression of adaptive genes and evolutionary radiations (Seehausen, 2004). Ours is the first study to assess these links in the context of the Australian flora; testing whether deep reticulate evolution is common or detectable in other Australian plant groups with similar radiations is an important next step. Explicitly testing for the adaptive function of introgressed genes in groups that experienced a rapid radiation and show signals of hybridization would further advance our understanding of the role that hybridization plays in the evolution of such groups (Suarez-Gonzalez et al., 2018).

Introgression between extant Adenanthos species in this study and others (Walker et al., 2018) show that reticulate evolution is ongoing. Porous species boundaries may still be evolutionarily advantageous in Adenanthos, resulting in intermediate phenotypes that can occupy different niches to the parental species (Givnish, 2010). This has been observed in Adenanthos, with hybrids observed to occupy intermediate or disturbed habitats (Nelson, 1977) and some taxa are known to be disturbance specialists (Groom and Lamont, 2015). Further studies are required to investigate whether ongoing gene flow across extant species plays an important role in the speciation of the genus (e.g., through homoploid hybrid speciation).

Recent Introgression and Hybridization Events

Several instances of recent introgression between extant Adenanthos species, while mainly isolated to closely related clade-specific lineages, were also detected in addition to strong incongruence predominantly shown across deeper scales across the backbone of the genus. These include potential chloroplast capture events where sympatric populations of one species share their plastid with another species with an overlapping geographic range (Rieseberg and Soltis, 1991). Examples include A. macropodianus–A. terminalis, A. ellipticus–A. cuneatus, and A. dobagii–A. oreophilus (Figures 2, 6). The chloroplast capture event for the southeastern Australian species occurred relatively recently in the Pleistocene (c. 1.3 Ma, 95% CI: 0.3–2.9 Ma) compared with the species divergence of the pair (A. macropodianus–A. terminalis) estimated at 15.3 Ma (95% CI: 8.2–23.0 Ma) based on the nuclear topology. Signals of ancient chloroplast capture events were also evident for A. stictus–A. cygnorum, and A. pungens/A.cuneatus–A. linearis, which have close geographic proximity and exhibit closely related chloroplast haplotypes but are phylogenetically distant in the nuclear phylogeny (Figures 2, 6). In these cases, the chloroplast reflects more the geographic patterning of these taxa instead of species relationships. Other hybridization events between extant species include Adenanthos × cunninghamii which shows conflicting placements in the nuclear and plastid topologies (sister to A. cuneatus and A. sericeus, respectively), corroborating the study by Walker et al. (2018) that it is a hybrid between A. cuneatus and A. sericeus. The cultivated A. sericeus is sister to A. dobagii, A. oreophilus and A. sericeus in both our nuclear and plastid topologies. We hypothesize that it is likely of a hybrid origin, with a parent species (likely maternal based on the plastid topology) from the A. sericeus clade and another undetermined parental species. Indeed the A. sericeus cultivar might be of multiple hybrid origins resulting from repeated back-crossing events, as two unsampled (or extinct) plastid haplotypes link it with the wild A. sericeus samples, and one unsampled haplotype with A. dobagii and A. oreophilus (Rieseberg and Brunsfeld, 1992). The putative A. cuneatus × ellipticus hybrid is sister to A. ellipticus in both our nuclear and plastid topologies but shares the same chloroplast haplotype with sympatric A. ellipticus and A. cuneatus individuals. Further studies applying extensive population sampling to assess for introgression across populations for these putative hybrids are required to confirm their status as well as the identities of their parent species.

Species Tree and Infrageneric Classification of Adenanthos

Extensive hybridization, as detected by our analyses, best explains the observed incongruence between nuclear and chloroplast phylogenies of Adenanthos. The nDNA topology is largely congruent with taxonomic concepts for the genus derived from morphology, whereas the plastid topology contains strong signals of multiple introgression events. This finding coupled with little to no detectable signal of introgression from our nuclear data warrants further discussion. It is possible that selection for adaptive organellar introgression or prevention of nuclear introgression could explain our results (Bonnet et al., 2017). Based on the simulaton study of Bonnet et al. (2017), these scenarios are the main drivers for this pattern where there is little evidence for nuclear introgression despite strong discordance between nuclear and organellar genomes. Local selection for different chloroplast genomes have been linked to different environmental performance of these genomes (Sambatti et al., 2008; Sloan et al., 2017). In sunflowers, for example, local adaptation to drier or wetter parts of a species’ range has contributed to multiple organellar introgression events across the genus (Sambatti et al., 2008; Lee-Yaw et al., 2019). It would be interesting to test in future studies whether the diverse chloroplast haplotypes and introgression events found within Adenanthos in SWA are the result of strong selection pressure for adaptation to local environmental conditions.

Many other studies have demonstrated strong cytonuclear discordance and separate evolutionary histories of plastid/mitochondrial vs. nuclear DNA (Soltis and Kuzoff, 1995; Yoo et al., 2002; Barrett et al., 2015) and some (e.g., Acosta and Premoli, 2010; Othman et al., 2010; Barrett et al., 2015; Folk et al., 2017; Vargas et al., 2017) have documented cases where plastid topologies are not reflective of species trees in comparison with nuclear data. Because the plastome is non-recombining and uniparentally inherited, a chloroplast lineage in one species can be replaced by an alien one following a single hybridization event (chloroplast capture), with the newly acquired chloroplast inherited by descendant lineages and persisting over long evolutionary timescales. This is expected to lead to plastome gene trees that are highly discordant with the species tree. Tree inference under the coalescent model using multiple nuclear loci is expected to provide a more accurate estimate of the species tree as compared to organellar genes that exhibit uniparental inheritance (Birky, 1995). For this reason, a cautious approach is needed when interpreting evolutionary signals between organellar and nuclear data. In particular, combining these datasets when they are in strong conflict will most likely compromise interpretations of evolutionary history.

While our nDNA topology is largely consistent with morphology in Adenanthos, nevertheless the infrageneric classification proposed by Nelson (1977) is partially inadequate. Nelson recognized two sections in Adenanthos, sect. Eurylaema and sect. Adenanthos, based on anther and style morphology, and two subsections within sect. Adenanthos based solely on perianth length. In our study, A. sect. Eurylaema was resolved as monophyletic in the nuclear topology but not the plastid topology, likely due to an ancient introgression event between A. obovatus and/or A. barbigera (sect. Eurylaema) with an extinct lineage from sect. Adenanthos. Both subsections of A. sect. Adenanthos were recovered as polyphyletic in both our plastid and nuclear datasets, from concatenated and coalescent analyses. We do not support the recognition of subsections within A. sect. Adenanthos and recommend they be merged.

Conclusion

Our study used complementary simulation approaches to detect introgression events across multiple scales within Adenanthos, and linked deep reticulate evolution to a rapid radiation in the Miocene coinciding with widespread aridification of the Australian continent. Dense sampling within Adenanthos allowed us to infer the extent and timing of introgression events within the genus. Reticulate signals were detected in a complex pattern of long-persisting haplotypes scattered across phylogenetically distant extant species. Some of these ancient chloroplast haplotypes are estimated to have diverged up to 12 Ma and may have persisted in southwestern Western Australia due to the relative stability of the landscape and buffering from major extinctions. Important open questions are the degree to which other Australian plant radiations show similar signals of reticulate evolution, and the effects of hybridization and introgression on their diversification.

Data Availability Statement

Datasets presented in this study can be found on the Dryad Digital Repository: doi: 10.5061/dryad.zw3r22876. Nge et al, (2020), data for the manuscript: reticulate evolution, ancient chloroplast haplotypes, and rapid radiation of the Australian plant genus - Adenanthos (Proteaceae), Dryad, Dataset: doi: 10.5061/dryad.zw3r22876.

Author Contributions

FN, EB, KT, and MW conceived the ideas for this study. FN designed the study. FN, EB, and KT collected the data. FN and EB complied the data. FN analyzed the data. FN wrote the manuscript with contributions from EB, KT, and MW. All authors contributed to the article and approved the submitted version.

Funding

Funding for this project was supported by the South Australian Department of Environment, Water and Natural Resources (D0004335204). FN was supported by an Australian Government Research Training Program (RTP) scholarship.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Hans Lambers, Marion Cambridge, Roberta Dayrell, Graham Zemunik, Rosalie and Robert Lawrence, Anne Rick and members of the Wildflower Society of Western Australia for their assistance in the field, Greg Keighery for sharing his knowledge on the current phrase-named taxa as well as hybrids within the genus, Dan Duval for providing material from the South Australian Seed Conservation Centre, Korjent van Dijk for guidance and help in laboratory work, Ryan Folk and Nicolás García for insightful discussions and help with the DendroPy simulation script, and Oscar Vargas for helpful discussions and clarification on the technicality of the JML software. We acknowledge the Department of Biodiversity, Conservation and Attractions (Western Australia; permit no. SW 019140), and Department of Environment, Water and Natural Resources (South Australia; permit no. G25787-4) for their permission to collect leaf samples on land under their administration. We also thank Amanda Shade and Fernanda Veraldo for permission and help in obtaining fresh leaf material from Adenanthos specimens grown at the Kings Park Botanic Gardens in Perth. Lastly, we would also like to thank the curation teams from both the Western Australian (PERTH) and South Australian Herbarium (AD) for their assistance and for providing access to herbarium specimens.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.616741/full#supplementary-material

Footnotes

1. ^ https://www.qiagenbioinformatics.com
2. ^ https://phytozome.jgi.doe.gov
3. ^ http://www.environment.gov.au/topics/land/national-reserve-system/science-maps-and-data/australias-bioregions-ibra

References