Population genomics and sexual signals support reproductive character displacement in Uperoleia (Anura: Myobatrachidae) in a contact zone

Abstract

When closely related species come into contact via range expansion, both may experience reduced fitness as a result of the interaction. Selection is expected to favour traits that minimize costly interspecies reproductive interactions (such as mismating) via a phenomenon called reproductive character displacement (RCD). Research on RCD frequently assumes secondary contact between species, but the geographical history of species interactions is often unknown. Population genomic data permit tests of geographical hypotheses about species origins and secondary contact through range expansion. We used population genomic data from single nucleotide polymorphisms (SNPs), mitochondrial sequence data, advertisement call data and morphological data to investigate a species complex of toadlets (Uperoleia borealis, U. crassa, U. inundata) from northern Australia. Although the three species of frogs were morphologically indistinguishable in our analysis, we determined that U. crassa and U. inundata form a single species (synonymized here) based on an absence of genomic divergence. SNP data identified the phylogeographical origin of U. crassa as the Top End, with subsequent westward invasion into the range of U. borealis in the Kimberley. We identified six F₁ hybrids, all of which had the U. borealis mitochondrial haplotype, suggesting unidirectional hybridization. Consistent with the RCD hypothesis, U. borealis and U. crassa sexual signals differ more in sympatry than in allopatry. Hybrid males have intermediate calls, which probably reduces attractiveness to females. Integrating population genomic data, mitochondrial sequencing, morphology and behavioural approaches provides an unusually detailed collection of evidence for reproductive character displacement following range expansion and secondary contact.

Keywords: mitochondrial genome; range expansion; reproductive interference; speciation; unidirectional hybridization.

FIGURE 1

(a) Maximum‐likelihood phylogeny of SNP data for all nonhybrid individuals of Uperoleia borealis and Uperoleia crassa, plus representatives of all other closely related species. Values on branches indicate support (SH‐like approximate likelihood ratio test [SH‐aLRT]/ultrafast bootstrap [UFboot]). (b) Result of the coalescent Bayesian snapp analysis, showing high support for species monophyly, but is unable to resolve the interspecific relationships in the clade. Values on branches represent the posterior probability

FIGURE 2

Bayesian coalescent phylogeny of the 16S and ND2 mitochondrial phylogeny (a). Values on nodes represent the support for each node from both the beast coalescent analysis and the maximum‐likelihood analysis (posterior probability [PP], SH‐like approximate likelihood ratio test [SH‐aLRT] and ultrafast bootstrap [UFboot]). Individuals identified as F₁ hybrids in the SNP analysis are marked with an asterisk. Distributions of mitochondrial clades are shown for Uperoleia borealis (b) and Uperoleia crassa (c). The asterisks on the maps correspond to the shaded clades in (a). These data identify only two well‐supported U. borealis clades in the mitochondrial DNA (vs. three in the SNP data set), and these clades overlap geographically (b). Two mitochondrial clades are identified in U. crassa, but these clades overlap geographically (c) and are differently distributed from populations in the SNP data set

FIGURE 3

Distribution maps, structure results and mitochondrial lineage assignments of Uperoleia borealis and Uperoleia crassa samples. The “Top End,” “Kimberley” and “Victoria River District” are identified as referenced in the main text. In the map in (a), U. borealis is presented in blue numbers, U. crassa is presented in yellow numbers and F₁ hybrids are presented in red numbers. Map numbers identify individuals in the structure result below. The barchart in (a) is the result of the Bayesian clustering of individuals from the program structure, which is based on allele frequencies of the nuclear DNA. Below this is the colour‐coded mitochondrial DNA clade of each individual at the species level (see Figure 2; black represents no mitochondrial sequence data). In (a), all hybrid individuals have a U. borealis mitochondrial DNA haplotype. Separate structure analyses are shown for U. borealis (b) and U. crassa (c). For each of these species, the within‐species mitochondrial lineage distributions differ substantially from the structure identified by the nuclear DNA

FIGURE 4

Geographical patterns of heterozygosity and range expansion in Uperoleia crassa. The observed heterozygosity values for individuals from the Kimberley and top end populations are shown as violin plots (left) and were significantly different (p < .001). The results of the rangeexpansion analysis (right), which identify the probable origin of expansion, show the distribution of individuals and their relative heterozygosity (dark circle is low, light circle is high). The probable expansion origin, as indicated by yellow colour and the X, is identified as the eastern top end

FIGURE 5

Variation in advertisement signals within the species complex. Waveforms are shown in (a)–(e); the scale bar in (a) applies to all waveforms. (b–d) Variation in Uperoleia crassa signals, which differ between allopatric (b) and sympatric (c and d) populations in the top end and Kimberley, respectively. Within sympatric populations, some individuals produce calls with a characteristic “final pulse” while other calls lack such a final pulse (compare panels c and d). (e) the advertisement call of a male confirmed to be a U. borealis/U. crassa hybrid using nuclear DNA. Violin plots show variation among groups (blue squares: U. borealis; yellow triangles: Sympatric U. crassa, green circles: Allopatric U. crassa) in pulse number (f), pulse rate (g) and call duration (h), with horizontal black bars showing group means. Star notation in (f)–(h) denotes when the starred group is statistically different from both other groups based on Tukey's range tests in Table S2. (i) the first two principal components from the PCA cumulatively explained 94.63% of the variation in temperature‐corrected advertisement call traits. Advertisement call data generally supported the prediction that U. crassa and U. borealis calls should differ more from each other in sympatry than in allopatry, and these differences were mainly driven by pulse characteristics

FIGURE 6

Geographical variation in advertisement call characteristics supports the reproductive character displacement hypothesis. (a) Map of advertisement call recording locations by species group. (b) PC1 (loaded largely by pulse characteristics) in relation to longitude. Uperoleia crassa and U. borealis differ more in sympatry than in allopatry. (c) PC2 (loaded predominately by call duration) in relation to longitude. (d) K‐means cluster analysis naïvely assigned U. crassa individuals to broad “western” (yellow) and “eastern” (green) clusters based on their advertisement call characteristics. The geographical “transition zone” at which frogs were more likely to be assigned to the eastern cluster than expected by chance was 127.65 degrees of longitude. Key at top of figure applies to (a)–(c)