Reticulate evolutionary history of the three subspecies
The pattern of genomic divergence provides clues to infer how the three subspecies evolved. A previous phylogenetic analysis indicated that the three subspecies split in a bifurcated manner, with australasicadiverging first at about 2.7 Mya and marina diverging fromeucalyptifolia at about 1.8 Mya (X. Li et al., 2016). However, given the high variability of genetic divergence across the genome we described above, such phylogenetic analyses using a handful of markers are not reliable in resolving taxa below species. Our Approximate Bayesian Computation modeling supports a simultaneous split of the three subspecies. Such trifurcate split is probably a nascent event of speciation radiation, although the number of diverging lineages is not that large.
Mantel tests indicate that geographical distance might have played a role in isolating populations of A. marina . Considering that the three subspecies are distributed along a continuum, we expect to observe a significant positive correlation between genetic divergence and geographical distance in between-subspecies population pairs. However, the divergence cannot be completely explained by isolation by distance. Obvious geographical barriers are inferred between Southeast Asia and Australasia, roughly consistent with the boundary of the marinaand eucalyptifolia ranges. Based on the clines described above, we hypothesize a scenario of marinaeucalyptifolia split: New Guinea was connected to Australia during glacial ages when sea level was low (Duke et al., 1998) and split A. marina populations east and west of the Torres Strait, followed by westward expansion ofeucalyptifolia during periods of high sea level and opening of the Torres Strait (Gordon, 2005; Hall, 2009). On the other side,eucalyptifolia may have differentiated from australasicabetween Rockhampton and Brisbane on the east coast of Australia via the bifurcation of the North Caledonian Jet into the North Queensland and the East Australian Currents (Ganachaud et al., 2007; Schiller et al., 2008) and the latitudinal change in environmental conditions such as temperature. Further studies may clarify these hypotheses.
The discordance between morphology and genetics in some populations (e.g. BB) hints at gene flow among subspecies and is supported by the TreeMix inferences. The BB population on the west coast of Australia is morphologically diagnosed as subspecies marina (Duke, 1991) but shows high genetic divergence from other populations and is genetically closer to eucalyptifolia than marina . Remarkably, all subspecies-informative loci we surveyed show either marina oreucalyptifolia haplotypes in BB. Hence, BB is highly likely an admixed population. Our simulation confirmed this. In addition, the CA and DW populations show very different patterns of SNP allele frequencies but have not diverged much. This contrast is probably due to both populations harboring composite alleles introduced from eithermarina or australasica . CA is likely admixed withaustralasica and DW with marina, since they neighbor each other. The sharing of haplotypes with respective subspecies of the two populations appears to support this speculation (Figure 4).
The three subspecies appear to have evolved reticulately, with gene flow among them. This is consistent with their subspecies designation, which indicates they are in the continuum of speciation. Gene flow during speciation has been reported in many taxa, especially at the early stages (Brandvain, Kenney, Flagel, Coop, & Sweigart, 2014; Clarkson et al., 2014; Harr, 2006; Poelstra et al., 2014; Wang, He, Shi, & Wu, 2020). How gene flow recedes as speciation nears completion remains to be addressed. Gene flow is recently proposed to exist even at later stages of speciation (Wang, He, Shi, & Wu, 2020). The role of gene flow in the evolution of populations at the subspecies stage would be interesting to pursue in further studies.