Polyploidy breaks speciation barriers in Australian burrowing frogs Neobatrachus

Abstract

Polyploidy has played an important role in evolution across the tree of life but it is still unclear how polyploid lineages may persist after their initial formation. While both common and well-studied in plants, polyploidy is rare in animals and generally less understood. The Australian burrowing frog genus Neobatrachus is comprised of six diploid and three polyploid species and offers a powerful animal polyploid model system. We generated exome-capture sequence data from 87 individuals representing all nine species of Neobatrachus to investigate species-level relationships, the origin and inheritance mode of polyploid species, and the population genomic effects of polyploidy on genus-wide demography. We describe rapid speciation of diploid Neobatrachus species and show that the three independently originated polyploid species have tetrasomic or mixed inheritance. We document higher genetic diversity in tetraploids, resulting from widespread gene flow between the tetraploids, asymmetric inter-ploidy gene flow directed from sympatric diploids to tetraploids, and isolation of diploid species from each other. We also constructed models of ecologically suitable areas for each species to investigate the impact of climate on differing ploidy levels. These models suggest substantial change in suitable areas compared to past climate, which correspond to population genomic estimates of demographic histories. We propose that Neobatrachus diploids may be suffering the early genomic impacts of climate-induced habitat loss, while tetraploids appear to be avoiding this fate, possibly due to widespread gene flow. Finally, we demonstrate that Neobatrachus is an attractive model to study the effects of ploidy on the evolution of adaptation in animals.

Fig 1. Independent origins of Neobatrachus tetraploids and high levels of reticulation.

(A) Gene trees, colored by clade, for 361 nuclear loci based on 2 individuals per species show considerable incongruence and differ from the species trees (bold black topology). Conflict between gene tree clusters (S3 Fig) and the nuclear species tree suggest non-bifurcating relationships between the species. (B) Pie charts represent summarised admixture proportions for each species (summing assignments for each individual, S1 Fig, Fig 2) at optimal clustering with K = 7. Tetraploids (N. sudellae, N. aquilonius and N. kunapalari) show highly admixed ancestries. (C) Dated diploid-only species tree. Colors represent consistency levels between gene genealogies with red being most conflicted and blue most consistent. Grey bars represent 95% confidence intervals on the ages of nodes, noted in millions of years before present.

Fig 2. ADMIXTURE results (K = 7) shown separately for each species.

According to the geographical locations of the sampled individuals, pie charts show the probability of the assignment of the individual to one of the 7 individually colored clusters. Overlapping pie charts on the map have been moved just enough to appear separate. Diploid Neobatrachus species (top 6: N. pelobatoides, N. albipes, N. wilsmorei, N. sutor, N. pictus, N. fulvus) are each assigned to separate clusters, while all three tetraploid species (bottom 3: N. kunapalari, N. sudellae, N. aquilonius) show inter-species admixture.

Fig 3. Widespread introgression between Neobatrachus species.

(A) Bifurcating maximum likelihood tree produced by TreeMix. (B) Example of a graph produced by TreeMix with 5 allowed migration events. (C) Scaled residual fit between observed data and predicted model in (A). Plot shows half of the residual covariance between each pair of populations divided by the average standard error across all pairs. Positive residuals represent populations where the model underestimates the observed covariance, meaning that populations are more closely related to each other in the data than in the modeled tree. Such population pairs are candidates for admixture events. Similarly, negative residuals indicate pairs of populations where the model overestimates the observed covariance. Overall, the residual plot of the model suggested that model fit could be improved by additional edges (migration events). (D) Scaled residual fit between observed data and predicted model in (B). Compared to Fig 3C this suggests that, although the complexity of the species relatedness is not fully represented by the model, major gene flow events and their direction were probably captured. (E) Box plots of 30 runs of TreeMix (each started with a different seed for random number generation) likelihood at different numbers of allowed migration events; saturation starts after 3 additional migration edges. (F) Bar plot showing the number of times a particular directional migration event was inferred in 30 TreeMix runs with 5 migration events allowed. We show only the events which were inferred more than twice.

Fig 4. Diversity and differentiation of Neobatrachus species and geographical suitability estimates.

(A) Example of the estimation of the suitable distribution area for N. sutor, based on occurrence data and current climate. (B) Example of the projection of the suitable distribution area for N. sutor based on the past climate at around 20Kya at LGM (last glacial maximum). Note that the scales in A and B are the same; Australian continent here is larger due to lower sea levels at LGM. (C) Scatter plot showing relative change of the predicted suitable area at the LGM and current conditions for each species as a function of Tajima’s D estimator. Diploid species show high correlation between Tajima’s D and distribution area change (blue line, Pearson’s correlation -0.88 (R² = 0.72, p-value = 0.02); (D) Hierarchical clustering of Neobatrachus species based on mean nucleotide diversity within and between the species.