Open Access This article is
- freely available
Genes 2019, 10(4), 317; https://doi.org/10.3390/genes10040317
Mitochondrial Introgression, Color Pattern Variation, and Severe Demographic Bottlenecks in Three Species of Malagasy Poison Frogs, Genus Mantella
CIBIO Research Centre in Biodiversity and Genetic Resources, InBIO, Universidade do Porto, Rua Padre Armando Quintas, N° 7, 4485-661 Vairão, Portugal
School of Biosciences, Cardiff University, Sir Martin Evans Building, Museum Avenue, Cardiff CF10 3AX, UK
Mention Zoologie et Biodiversité Animale, Faculté des Sciences, Université d’Antananarivo, BP 906, Antananarivo 101, Madagascar
Zoological Institute, Technische Universität Braunschweig, Mendelssohnstr. 4, 38106 Braunschweig, Germany
Author to whom correspondence should be addressed.
Received: 6 February 2019 / Accepted: 15 April 2019 / Published: 23 April 2019
Madagascar is a biodiversity hotspot particularly rich in amphibian diversity and only a few charismatic Malagasy amphibians have been investigated for their population-level differentiation. The Mantella madagascariensis group is composed of two rainforest and three swamp forest species of poison frogs. We first confirm the monophyly of this clade using DNA sequences of three nuclear and four mitochondrial genes, and subsequently investigate the population genetic differentiation and demography of the swamp forest species using one mitochondrial, two nuclear and a set of nine microsatellite markers. Our results confirm the occurrence of two main mitochondrial lineages, one dominated by Mantella aurantiaca (a grouping supported also by our microsatellite-based tree) and the other by Mantella crocea + Mantella milotympanum. These two main lineages probably reflect an older divergence in swamp Mantella. Widespread mitochondrial introgression suggests a fairly common occurrence of inter-lineage gene flow. However, nuclear admixture seems to play only a limited role in this group, and the analyses of the RAG-1 marker points to a predominant incomplete lineage sorting scenario between all five species of the group, which probably diverged relatively recently. Our demographic analyses show a common, severe and recent demographic contraction, inferred to be in temporal coincidence with the massive deforestation events that took place in the past 1000 years. Current data do not allow to conclusively delimit independent evolutionary units in these frogs, and we therefore refrain to suggest any taxonomic changes.
Keywords:hybridization; Mantella aurantiaca; M. crocea; M. milotympanum; M. madagascariensis; M. pulchra; cytochrome b; RAG-1; microsatellites
Madagascar is home to an almost unparalleled diversity of flora and fauna, with the majority of its native species being endemic to this island . Amphibians are represented on Madagascar by five frog clades of independent origin [2,3]. At least some of these clades are very rich in species, such as the endemic family Mantellidae and the endemic microhylid subfamily Cophylinae [4,5], but a large portion of this species diversity remains to be taxonomically described [6,7]. While the taxonomic knowledge on Malagasy amphibians has substantially progressed in the last thirty years [2,8,9] and a large research effort in more recent times has been directed to study the processes of species diversification acting on these radiations [4,10,11,12,13], the population-level differentiation of Malagasy amphibians is still poorly explored. Madagascar is home to many endemic and endangered frog species  and this approach can provide useful information to improve the management of the many threatened species. So far, this important aspect of conservation genetics has been investigated mainly in a few charismatic species such as the large tomato frogs , the rainbow toad , and some Malagasy poison frogs of the genus Mantella [16,17].
Phylogeographic studies in amphibians often identify limited gene flow among populations and a genetic structure that correlate with their geography and distribution [18,19,20], which usually is interpreted as a signal of low vagility. In Madagascar poison frogs, this pattern has been observed in some species (e.g., Mantella bernhardi, M. expectata and, to a shallower extent also in Mantella viridis [15,16,21]). On the contrary, other studies have identified a surprisingly low degree of phylogeographic structure and common mitochondrial haplotype sharing in Mantella species of the eastern rainforest [22,23]; rapid range expansions in other Mantella species and in the tomato frogs of the genus Dyscophus [14,17], and widespread co-occurrence of deeply divergent mitochondrial haplotype lineages within the same population that seems to characterize an arid adapted Mantella species [24,25].
Due to their usually bright aposematic coloration (a case of convergent evolution with South American poison dart frogs [26,27,28], diurnal activities, and high prevalence in the pet-trade [29,30] Mantella are probably among the best-studied group of Malagasy amphibians [17,21,22,23,24,31,32,33]. Yet, information on their distribution and habitat preferences is still lacking for many species [34,35,36,37], and several studies suggest the genus needs a taxonomic revision [16,19,21,24,25].
Within the genus Mantella it is striking that often, sympatric species share very similar color pattern, which might indicate phenomena of Müllerian mimicry . This is most obvious in the species pair M. baroni and M. madagascariensis which belong to different species groups; but also the syntopic M. laevigata and M. manery share a dorsal pattern similar to each other, as do the uniformly orange M. aurantiaca and M. milotympanum which live in areas only few km apart. On the other hand, species that apparently are very closely related to each other can strongly differ in coloration [17,21,22] pointing to a probably fast and repeated evolution of chromatic characteristics . Most of the studies to date were entirely based on mitochondrial genes or single nuclear markers which might be affected by introgression phenomena that could blur true relationships; in fact, hybridization has already been confirmed between M. cowani and M. baroni , and recent gene flow was detected between the morphologically distinct M. nigricans, M. baroni and M. aff. baroni .
In this study we focus on the Mantella madagascariensis group that is characterized by a particularly high color variation, both among species and among supposedly conspecific populations. We deem this characteristic favoured the description of different color morphs as distinct species (e.g., in the case of the poorly differentiated M. crocea Pintak and Böhme, 1990 and M. milotympanum Staniszewski, 1996) and we investigate the population genetics of this group to explore this hypothesis. According to current taxonomy, this group consists of two rainforest species, M. madagascariensis and M. pulchra, and three swamp forest species M. aurantiaca, M. crocea and M. milotympanum . Mantella aurantiaca is characterized by a uniform and translucent orange coloration (Table S1). Mantella milotympanum from the type locality near Fierenana is opaque red orange with a black tympanic region (the characteristic that is also giving the name to this species) and a ventral surface uniformly orange. Mantella crocea shows a lateral black coloration, a frenal stripe, a ventral black color with a variable number of markings that can be grey, bluish-white or yellowish, and a dorsal color being yellow-brown in the more southern populations and bright green in the populations of the western and northern portion of its range. Populations with intermediate patterns between Mantella crocea and M. milotympanum (Ambatombolana, Andriambe, Savakoanina and Andaingo ) were assigned to M. milotympanum (Ambatombolana, Andriambe, Savakoanina) based on the presence of the black tympanic region and a uniformly orange ventral surface and on the lack of both lateral black coloration and frenal stripe. The population of Andaingo was assigned to M. crocea based on the presence of a lateral black coloration, a frenal stripe and a ventral black coloration with sparse grey, yellowish-white markings (see Table S1 for all details). Despite some overlap in habitat specialization of rainforest vs. swamp forest species, the latter are most common in a stretch of the northern central east of Madagascar characterized by extended swamps surrounded by forests, often with many Pandanus screw palms (Figure 1). These species are characterized by distinct alkaloid profiles , probably due to low prey electivity and specific arthropod sources in their habitat . As far as known, species in this group usually do not occur syntopically, but their ranges have various contact zones. These frogs are known to hybridize in captivity and haplotype sharing has been reported among Mantella crocea and M. aurantiaca, and among M. crocea and M. milotympanum [22,23]. Mantella milotympanum and M. aurantiaca are currently categorized as Critically Endangered following the red-listing criteria of the International Union for Conservation of Nature (IUCN) . This is mostly due to their restricted and severely fragmented distribution (with an area of occupancy of less than 10 km2), the ongoing decline in the extent and quality of their habitat and their possible over-exploitation in the pet-trade. On the contrary, M. crocea is listed as of Least Concern due to its relatively wide extent of occurrence and presumed large population . Due to the widespread mitochondrial haplotype sharing between M. crocea and M. milotympanum (and to a lesser extent between M. crocea and M. aurantiaca) and their radically different IUCN Red List categorization it is important to investigate the degree of differentiation of these species and assess their taxonomic status, given that most studies to date focused on only mitochondrial genes and did not cover the entire known distributional range of these species.
Although there is still no consensus on the time of human colonization of Madagascar ([42,43] and literature cited therein for more details), a recent study finds no support for human activity before 1500 y B.P., and otherwise finds diverse evidence of human activity from about 1350 y B.P. . It is expected that deforestation increased sharply after human colonization and that this markedly affected the demographic history of Madagascan native species. To test this hypothesis, we investigated the demographic history of the populations of the three swamp forest species of the M. madagascariensis group for which we had a complete extended molecular dataset.
To summarize, in this study we investigate the genetic differentiation and demography of species in the M. madagascariensis group, in particular targeting the three swamp forest species. We first ascertain monophyly of this this focal group of species by a genus-wide analysis of multiple nuclear and mitochondrial DNA sequences, and subsequently analyse population genetic differentiation in several species of the group based on an extended data set with samples from new localities, and with the inclusion of new molecular markers, RAG-1, RAG-2 and a set of nine microsatellites .
2. Materials and Methods
2.1. Tissue Sampling, DNA Extraction, and Amplification
Two-hundred and fifteen tissue samples of individuals of the M. madagascariensis group (sensu ) were analyzed in this study, including: 68 samples from four localities of M. aurantiaca, 44 samples from four localities of M. milotympanum, 72 samples from eight localities of M. crocea, 22 samples from three localities of M. madagascariensis, and 9 samples from one locality of M. pulchra (Figure 1; Tables S1 and S2 for more details). These covered a large proportion of the known distribution of M. aurantiaca, M. milotympanum and M. crocea [2,35], with the northernmost population of M. crocea in the Makira Reserve (locality Hevirina) extending the known distribution of this species substantially. Twenty-three samples of Mantella spp. representing all the currently described species and the recently identified cryptic lineages were used for the phylogenetic analyses (details below).
Tissue samples were either toe clips from subsequently released individuals, saliva swabs (M. crocea from Zahamena) or femur muscle from collected voucher specimens, all collected between 2002 and 2009 during several field expeditions along the eastern coast of Madagascar and stored in 96% EtOH. Sampling and tissue sample export permits were approved by the Ministère de l’Environnement et des Forêts of Madagascar.
Total genomic DNA was extracted using a proteinase K (10 mg/µL) digestion following a standard saline method . For the phylogeographic analyses (dataset available in Tables S1 and S2), polymerase chain reactions (PCR) were used to amplify a fragment of the mitochondrial cytochrome b (COB) gene of ca. 600 nucleotides, a fragment of the nuclear recombination activating gene 1 (RAG-1) of ca. 750 nucleotides, and a fragment of the nuclear recombination activating gene 2 (RAG-2) of ca. 500 nucleotides (Table S3). In addition to COB, RAG-1 and RAG-2, a fragment of the brain-derived neurotrophic factor (BDNF), and of the mitochondrial cytochrome oxidase subunit 1 (COX1), 12S rRNA and two fragments of 16S rRNA genes were sequenced for all known Mantella species to ascertain the monophyly of the M. madagascariensis group within the genus (for primer sequences, see Table S4).
PCRs were performed in a final volume of 11 µL and using 0.3 µL each of 10 pmol primer, 0.25 µL of total dNTP 10 mM (Promega, Madison, WI, USA), 0.08 µL of 5 U/mL GoTaq®, and 2.5 µL 5x Green GoTaq® Reaction Buffer (Promega). Primers and PCR amplification conditions in Table S4. The successfully amplified double-stranded PCR products, treated with Exonuclease I (New England Biolabs, Ipswich MA, USA) and Shrimp Alkaline Phosphatase (Promega) to inactivate remaining primers and dNTPs, were directly used for the cycle-sequencing reaction using dye-labeled dideoxy terminators (Applied Biosystems, Foster City, CA, USA) with the amplification primers on an ABI 3130xl automated DNA sequencer. All nuclear genes were sequenced in both directions.
The same tissue samples that were used for the phylogeographic analyses were genotyped for nine microsatellite loci (Table S1). Details on the primer combinations used, the amplification, and the genotyping protocols are described in . For a complete list of species, collecting localities, coordinates, gene fragments and microsatellite used see Tables S1 and S2.
2.2. Sequence Alignment, Statistics, and Phasing
Chromatographs were checked and sequences were edited and aligned using Codon-Code Aligner (v. 2.0.6, Codon Code Corporation, Centerville, MA, USA). A homologous COB sequence of M. baroni (AY862201) was added to the COB alignment. A Neighbor-Joining (NJ) tree based on uncorrected p-distances was constructed with MEGA, version 7.0.9  to gain a first overview of differentiation among sequences (Figure S1). We also added to this alignment 17 unpublished shorter COB sequences of M. aurantiaca from the captive-bred colony housed in the breeding facility managed by Association Mitsinjo in Andasibe) (Figure S1). All newly obtained sequences ware deposited in GenBank (MK422202-MK422436; MK440330-MK440540).
We discarded the hypothesis of pseudogenes occurring in the COB sequence dataset by confirming the absence of stop codons in the analyzed fragment. Because the RAG-1 and RAG-2 alignments contained heterozygous positions we phased genotypes to identify haplotypes for further analyses using the PHASE algorithm (version 2.1.1) with default settings  as implemented in the software DnaSP (version 6.12.01; ). PHASE parameters were 1000 iterations, one thinning interval and 100 burn-in iterations and a posterior probability threshold of 0.9 to determine the most probable inferred haplotypes for each nuclear sequence. Analyses were repeated three times using different seed values. We determined the number of base substitutions using DnaSP for the COB, RAG-1 and RAG-2 fragments and computed Tajima’s average number of pairwise differences (∏; ) and haplotype diversity (Hd; ) for each locality and for all fragments used for the phylogeographic analyses (Table S2). We also computed Tajima’D , computed as the difference between the mean number of pairwise differences and the number of segregating sites, for the populations that were used to investigate the demographic history of the three swamp forest species (details below).
2.3. Phylogenetic Analysis
DNA sequences of four mitochondrial and three nuclear genes were sequenced from one individual (voucher number available in Figure 2) of each Mantella species, to verify monophyly of the M. madagascariensis group, which contains the focal taxa of this study. Sequences were aligned as described above (nuclear sequences not phased) and concatenated for analysis. Because mitochondrial introgression has previously been detected in Mantella  we analyzed nuclear and mitochondrial genes separately. The alignment of the three nuclear genes was analyzed as a single partition to avoid overparametrization, under a TIM1+I substitution model selected by the AIC criterion in jModeltest . For the concatenated mitochondrial genes, a partition scheme with seven partitions was selected with PartitionFinder , as follows: (1) 1st positions of COB and 3rd position of COX1, HKY+I; (2) 2nd position of COB, GTR+G; (3) 12S and 3rd position of COB, GTR+I+G; (4) 1st position of COX1, GTR+G; (5) 2nd position of COX1, SYM+I; (6) 16S, GTR+I; (7) intervening tRNAs (tRNA Leut, tRNA Val) between 16S and 12S, GTR+G. Bayesian inference (BI) phylogenetic analysis was subsequently run in MrBayes , with 50 million generations and a burn-in of 25% (Figure 2).
2.4. Haplotype Network Reconstruction
Individual COB sequences and phased allelic sequences for the RAG-1 and RAG-2 gene fragments were collapsed into haplotypes using the online web tool DnaCollapser 1.0 available on the FaBox platform (http://users-birc.au.dk/biopv/php/fabox/). Prior to calculating haplotype relationships using networks, we analyzed the nuclear RAG-1 and RAG-2 datasets with the Phi test of Bruen et al.  as implemented in SplitsTree4  to assess the presence of recombination in the data. Neither of the genes presented significant evidence of recombination (p-values > 0.1) and thus they were used for network reconstruction. We used TCS version 1.21  to analyse the relationships of the COB haplotypes and of the RAG-1 and RAG-2 inferred haplotypes (Figure 3, Figure 4 and Figure 5). With these networks we show instances of haplotype sharing and highlight haplotype differences between and within analyzed populations, providing an intuitive graphic representation of the observed differentiation between populations, but being aware that phased haplotypes might not in all cases be correct due to the probabilistic nature of the Phase algorithm, and due to the large haplotype variation, the reconstructed haplotype relationships might not be fully accurate. The method underlying TCS  calculates the number of mutational steps by which pairwise haplotypes differ, computing the probability of parsimony for pairwise differences until the probability exceeds 0.95. The manual adjustment of the threshold was necessary to infer the number of mutational steps that were needed to join haplotypes into one single network for the analyzed COB fragment. Loops that occurred in the network were resolved by (1) treating connections with singletons and/or rare haplotype to be less likely than connections with central and frequent haplotypes  and (2) using the geographical criterion, favouring connection between haplotypes of the same or close population over connection with haplotypes from distant populations .
2.5. Microsatellite Analyses
The presence of null alleles was tested with MICRO-CHECKER . Linkage disequilibrium between loci was tested using ARLEQUIN 220.127.116.11  with 104 permutations to assess significance. The average number of alleles per locus, allelic richness (AR), expected (HE) and observed (HO) heterozygosities and inbreeding coefficients (FIS) were calculated using MSA v. 4.05  (Table 1). Allelic richness was calculated based on a minimum of 11 individuals per population. GENEPOP  was used to assess the Hardy–Weinberg equilibrium of each locus and population using 1000 dememorizations, 100 batches and 1000 iterations per batch.
2.6. Population Structure
The mitochondrial DNA based analyses reflect only the maternal genealogy relating these Mantella populations; therefore, we carried out an analysis of population structure using the nine nuclear microsatellites to attempt to identify the best way of grouping the individuals with a dataset with higher statistical power and representative of both parental lineages. This analysis was carried out with the software Bayesian Analysis of Population Structure (BAPS)  (Figure 6). Only individuals with ~10% or fewer missing data were kept for this analysis to reduce potential noise caused by the missing data (e.g., spurious clusters). This reduced the microsatellite dataset to 129 samples (Table 1). A range of possible values (1 to 10) for the optimal number of clusters in the data (K) was tested using the group-based approach in BAPS, and a-priori grouping individuals according to their sampling localities. The analysis was repeated three independent times to assess repeatability of the results. A NJ tree describing the similarity between sampling localities was built using Phylip  with 100 bootstrap replicas to assess branching pattern statistical support on the basis of the genetic distance (Figure 7). Proportion of Shared Alleles  was calculated in MSA. Pairwise FST between sampling localities was also computed with MSA with 104 permutations to assess significance (Table 2). We also used BayesAss v. 3.0.4  to estimate migration patterns between populations. BayesAss uses an MCMC algorithm to estimate the proportion of individuals form a population that derives from a different source population. This is achieved by estimating the likelihood of ancestry for each individual from each population analyzed using assignment methods that use the information on the individual’s genotype and the allele frequency distributions in the different populations. BayesAss’ MCMC running parameters were burn-in for 25 million steps, followed by 25 million steps of data collection saving each 1000th step. BayesAss was run in triplicate to assess repeatability of the results inferred with the MCMC algorithm. Convergence of the Markov chains was determined using Tracer from the BEAST suite .
2.7. Demographic History
Due to the conservation concerns on these species we carried out analyses of the demographic history of these populations to determine whether they present robust effective populations sizes that have been stable across time and thus maintained their evolutionary potential. For this purpose, the demographic history of the nine localities with the largest sample sizes and no more than 10% missing data was estimated with Msvar 1.3 [70,71]. Msvar uses the counts of alleles for each locus in a population to estimate whether the current effective population size (N0) is explained by a demographic change (t) generations ago from an ancestral effective population size (Nt), e.g., an expansion or a bottleneck. Msvar uses coalescent simulations conditioned on prior distributions for each of the demographic parameters, as well as the mutation rate of the microsatellites. For these analyses we carried out runs under different prior combinations representing bottlenecks, expansions and stable demography (no change in effective population size) for each population to determine that the posterior estimates of the parameters were not dependent on the underlying demographic model. The prior and hyperprior distributions used are shown in Table S5. A generation length of one year was assumed for all species and an average vertebrate microsatellite mutation rate of 10−4 [72,73,74] allowing it to vary between 10−2 and 10−6. Each Msvar was run for each demographic scenario for 400 million steps of the MCMC algorithm. For each sampling locality three independent runs (one under each demographic prior) were analyzed for convergence using Gelman and Rubin’s statistic  implemented in the library CODA  in R , after discarding the initial 50% of the MCMC steps as burn-in phase.
3.1. Monophyly of the M. madagascariensis Group
To understand if the target species of this study, defined as the M. madagascariensis group (M. aurantiaca, M. crocea, M. madagascariensis, M. milotympanum, M. pulchra) together form a monophyletic group within the genus Mantella, we performed a phylogenetic analysis based on DNA sequences of multiple genes. To exclude effects of mitochondrial introgression, we analyzed mitochondrial and nuclear genes separately. The phylogenetic trees obtained from BI analysis of the concatenated nuclear and of the mitochondrial data sets (Figure 2a,b) both resolved the deep phylogenetic relationships in the genus concordantly. Mantella bernhardi was placed as sister group to the remaining species of Mantella; of these, the M. cowani group (M. baroni, M. cowani, M. haraldmeieri, M. nigricans) was sister to a clade containing the M. betsileo group (M. betsileo, M. ebenaui, M. expectata, M. viridis), M. laevigata group (M. laevigata, M. manery) and M. madagascariensis group. The nuclear DNA analysis failed to unambiguously resolve relationships among the M. laevigata and M. betsileo groups, probably due to missing data for M. manery for which only a RAG-1 sequence was available. Similarly, relationships within the species groups did not receive significant support in the nucDNA analysis. However, both analyses provided high support for monophyly of the M. cowani group and of the M. madagascariensis group.
3.2. DNA Sequence Diversity and Differentiation within the M. madagascariensis Group
The total alignment length (after trimming the ends of the alignment to remove missing data) for the three genes was: COB, 215 individual sequences and a total length of 552 nucleotides; RAG-1, 211 individual sequences (422 after phasing) and a total length of 732 nucleotides; RAG-2, 21 individual sequences (42 after phasing) and a total length of 531 nucleotides. In all the three alignments we did not detect any gaps, and the translation into amino acids resulted in neither nonsense nor stop codons. In the COB alignment we found 81 (14.7%) variable sites (out of 552), among which 69 were parsimony-informative and 12 were singletons, and a total of 49 haplotypes were identified (C01-C49; see Tables S1–S3 for more details). In the phased RAG-1 alignment 59 sites (out of 732) were variable (8.1%), among which 51 were parsimony-informative and 8 were singletons, and a total of 78 haplotypes were identified (R01–R78; Tables S1–S3). In the phased RAG-2 alignment 27 sites (out of 531) were variable (5.1%), among which 16 were parsimony-informative and 11 were singletons, and a total of 17 haplotypes were identified (01–17; Tables S1 and S3). In the analyzed COB fragment the Hd values of the populations of M. crocea were high (ranging from 0.6 and 1 in M. crocea), while the Hd values of the analyzed populations of M. aurantiaca and M. milotympanum were very variable, ranging from 0 to 0.684 in M. milotympanum and from 0 to 0.752 in M. aurantiaca (Table S2). The nucleotide diversity (π) values for the localities provide an estimate on the genetic variability in a given data set. In the COB alignments the M. aurantiaca populations from Andranomandry and Besariaka show the highest values of π (see Table S2). In the RAG-1 and RAG-2 fragment estimated Hd values were on average high, ranging from 0.533 and 1 in all analyzed populations (Table S2). RAG-1 nucleotide diversity (π) values were on average lower than the estimated π at the COB fragment for M. aurantiaca and M. crocea and higher for M. milotympanum, M. madagascariensis and M. pulchra (see Table S2). Although the number of analyzed sequences per locality differed across the different localities (ranging from 2 to 52 in the phased RAG-1 alignment), it is worth noting that in this alignment the highest π was observed in the population from Ampangadimbolana of M. crocea (see Table S2). Nucleotide and haplotype diversity values for each locality per each analyzed locus are provided in Table S2. The values of Tajima’s D computed for the COB fragment are mostly negative (see Table S2), suggesting population size expansion although this result is likely to reflect the occurrence of widespread mitochondrial introgression between the analyzed populations. On the other side, the values of Tajima’s D computed for the RAG-1 and RAG-2 fragments are mostly positive (see Table S2), and are in overall agreement with the results obtained in the demographic analyses performed with the microsatellite data (see details in Section 3.6), suggesting a decrease in population size.
The NJ tree (Figure S1) based on the COB sequences suggests several instances of introgression of the M. crocea mitochondrial genome into the populations of M. aurantiaca of Andranomandry and Besariaka, and of introgression of M. pulchra mitochondrial genome into the population of M. aurantiaca of Torotorofotsy (see also Figure S1). A probable introgression of an M. aurantiaca mitochondrial genome was identified for the population of M. crocea from Ambohitantely, but the respective sequence differed by four mutations from the nearest haplotype identified in the extant populations of M. aurantiaca sampled by us.
The COB sequences of M. milotympanum and M. crocea were poorly differentiated from each other, and in this group the most differentiated clade was composed by sequences originated from the recently identified population of M. crocea in Hevirina (in the Makira area), the northernmost population known for this species. This analysis also shows that the northern populations of M. madagascariensis share the mitochondrial genome with M. pulchra (Figure S1). The 17 shorter sequences of M. aurantiaca coming from the captive bred colony showed almost no genetic differentiation to the population of M. aurantiaca from Totorotofotsy (Figure S1). As already observed in some individuals of M. aurantiaca from Andranomandry and Besariaka, also one individual of this pool carried the mitochondrial genome of M. crocea (Figure S1), suggesting that introgression has also taken place in these northern populations of M. aurantiaca, and not only in the southern ones.
3.3. Haplotype Networks
The haplotype network of the COB alignment led to the identification of three unconnected subnetworks plus one additional disconnected haplotype (C41) (Figure 3). One network was mostly composed of sequences of M. madagascariensis and M. pulchra, a second network contains all the sequences of M. milotympanum, most of the sequences of M. crocea (both from the brown and green phenotype populations) and some sequences of M. aurantiaca from Andranomandry and Besariaka. The third network is mostly composed of sequences on M. aurantiaca from its entire distributional range, but it also includes one haplotype of M. crocea from Ambohitantely (green phenotype; C17; see Table S3). Haplotype C41 requires at least 11 mutational steps to be linked to the M. aurantiaca network (haplotype C10), and at least 23 mutational steps to be linked to the M. crocea/M. milotympanum network (haplotype C13; Figure 3). Haplotype C41 is shared among all the analyzed individuals of M. madagascariensis from Ranomafana suggesting that this population is quite differentiated from the other analyzed populations, although the limited number of included populations and specimens of this population limits our conclusions. Twenty-three mutational steps are needed to link the network of the M. pulchra/M. madagascriensis to the network of the M. crocea/M. milotympanum. In the M. crocea/M. milotympanum network, the individuals of the population from Hevirina (that share the green phenotype with specimens from Ambohitantely, Zahamena and Andaingo) group together, and 10 mutational steps link this group of haplotypes to all the others. See Tables S1–S3 and Figure 1 and Figure 3 for more details on haplotype frequencies and distribution.
Only one network was retrieved for the phased RAG-1 alignment, which included all the five species of the M. madagascariensis group (Figure 4). The three individuals of M. aurantiaca from Andranomandry that show sign of mitochondrial introgression with M. crocea (see Figure 3, Tables S1 and S3) share two of the most widespread RAG-1 haplotypes (R1, shared by all species of the M. madagascariensis group, and R3, shared by all species of the M. madagascariensis group with the exception of the southern population of M. madagascariensis, see Table S1 and Figure 4), but also haplotype R6, which is exclusive of M. aurantiaca populations of Andranomandry and Besariaka. In the population of M. aurantiaca from Besarika, nine of the 21 analyzed individuals show mitochondrial introgression with M. crocea (see Figure 3, Tables S1 and S3). Most of these introgressed individuals, for RAG-1, bear haplotype R2, which is shared by individuals of all the analyzed localities of M. aurantiaca, but also by M. milotympanum, M. crocea (green phenotype populations) and M. madagasariensis. Haplotype R10 and R15 are uniquely found in these individuals. Other RAG-1 haplotypes found in these introgressed individuals are: R1, widely shared by all analyzed species; R3, shared by all species of the M. madagascariensis group with the exception of the southern population of M. madagascariensis; R4, mostly found in the brown populations of M. crocea but also found in M. milotympanum, M. aurantiaca and in the southern population of M. madagascariensis; R5, mostly found in M. aurantiaca populations from Andranomandry and Besariaka, but also found in the M. crocea population from Andasibe and in the population of M. madagascariensis from Ranomafana; R6, exclusive of M. aurantiaca populations of Andranomandry and Besariaka; R12, mostly found in M. milotympanum but also found in the M. aurantiaca population of Besariaka, in the population of M. crocea of Ambohitantely (green phenotype) and in the southern population of M. madagascariensis (Figure 4). The individual of M. aurantiaca from Torotorofotsy which shows signs of mitochondrial introgression from M. pulchra, shares RAG-1 haplotype R30 with another M. aurantiaca individual from Torotorofotsy, with M. crocea from Ambodivoasary and Andasibe (brown phenotype) and with M. pulchra (Figure 4). In the RAG-1 network, the individual of M. crocea from Ambohitantely with introgressed M. aurantiaca mitochondrial genome shares haplotype R1 with all the other species of the M. madagascariensis group, and R28 mostly with individuals of M. crocea from Ambohitantely (green phenotype) and Ampangadimbolana (brown phenotype), but also with another individual of M. aurantiaca from Torotorofotsy (Figure 4).
The specimens of M. madagascariensis from the North (Besariaka) mostly share haplotype R18 (found also in one population of M. aurantiaca and M. milotympanum) and R19 (found also in one population of M. aurantiaca and in the population of M. madagascariensis of unkown origin). Other RAG-1 haplotypes found in this population are: R21, also found in one population of M. aurantiaca and M. milotympanum and in the population of M. madagascariensis of unkown origin); R67, also found in the population of M. madagascariensis of unkown origin; and R68, uniquely found here (Figure 4). In the population of M. madagascariensis of unknown origin (but showing a pattern of genetic variability similar to the population of M. madagascariensis from Besariaka) we found the following RAG-1 haplotypes: R1, R2, R3, R19, R21 and R67 as discussed above; R65, found also in the population of Andriambe of M. milotympanum and M. pulchra; and R73, uniquely found in this population (Figure 4).
All the RAG-1 haplotypes group together in a single network, but there are a few groups of haplotypes that are both uniquely found in one population and are separated by at least 3 mutational steps by all the rest: (i) haplotypes R36, R40, and R41 are uniquely found in the population of M. crocea from Ambohitantely; and (ii) R57 and R56 in the M. crocea population of Hevirina (both of them showing the green phenotype); (iii) haplotype R47 was found in the population of M. crocea from Ampangadimbolana (with the brown phenotype); (iv) haplotypes R10, R13 and R17 plus haplotypes R49 and R50, and haplotype R59, are found in the population of M. aurantiaca from Besariaka, M. crocea from Ampangadimbolana (brown phenotype) and M. crocea from Zahamena (green phenotype), respectively. Similar to what was observed in the COB network, the population of M. crocea showing a more consistent genetic differentiation from the others is the northernmost known population of Hevirina (cf. Figure 3 and Figure 4; RAG-1 haplotypes R56 and R57). See Tables S1 and S3 and Figure 4 for more details on RAG-1 haplotype frequencies and distribution.
The RAG-2 fragment was amplified only for a small subset of individuals from six populations (some samples of M. aurantiaca from Andranomandry, of M. crocea from Ampangadimbolana, of M. milotympanum from Fierenana-Sahamarolambo, Andriambe and Savakoanina, and of M. madagascariensis from Besariaka), mostly to verify if other nuclear genes display a pattern of differentiation comaparable to that in RAG-1. The phased RAG-2 sequences grouped into a single network, although the genetic differentiation of the RAG-2 haplotypes might be slightly higher compared to RAG-1 (Figure 5). See Tables S1 and S3 and Figure 5 for more details on RAG-2 haplotype frequencies and distribution.
3.4. Population Genetic Variation in Microsatellites
We tested for Hardy–Weinberg Equilibrium (HWE) in the nine microsatellite markers studied by sampling locality. At all localities included in the microsatellite analysis, except for the M. crocea population from Ambodivoasary, the data revealed absence of HWE (Fisher Exact test p-value < 0.001). In all cases this lack of HWE was due to a deficit of heterozygotes, in particular in loci G10X, B8 and M14. Such pattern could be indicative of the presence of null alleles on these loci, however, we found no significant difference in estimates of genetic variation and inbreeding coefficients when excluding them (Welch corrected t-tests p-value > 0.1), thus we retained these loci for further analyses. No significant linkage disequilibrium within any of the sampling localities was detected after Bonferroni correction for multiple testing, with the exception of the pairwise comparison between loci A9 and B7 in the M. milotympanum population of Andriambe (p-value < 0.0001). Excluding the M. crocea sampling site of Ambodivoasary due to its small sample size, the number of alleles per locus ranged between a minimum of 6.11 and 12.78 across localities, the observed and expected heterozygosity were relatively high, ranging between 0.46 and 0.74, and 0.77, and 0.88, respectively (Table 1). The inbreeding coefficient ranged between a minimum of 0.14 and a maximum of 0.45 (Table 1) in the different localities.
3.5. Population Structure
BAPS identified seven clusters as the most probable clustering solution in the data using prior information on the sampling locality of the individuals (PP, posterior probability of each of the analysis replicates > 0.99). The same results were obtained for three independent BAPS runs indicating strong support for the pattern of population structure identified. These included three clusters of M. aurantiaca, i.e., (1) Andranomandry and Besariaka, (2) Sahasarotra, and (3) Torotorofotsy. The M. crocea samples formed three additional clusters: (4) Ambodivoasary, (5) Ambohitantely, and (6) Ampangadimbolana. Lastly, cluster (7) was formed by the M. milotympanun samples from Fierenana-Sahamarolambo, Andriambe and Savakoanina (Figure 6). With the exception of one individual of M. aurantiaca from Besariaka, no other sample presented significant admixture (Figure 6). The genetic variation of this admixed individual resembles to 50% cluster 1, 32% cluster 6, 11% cluster 7, 4% cluster 2, 2% cluster 4, and 1% cluster 3 (Figure 6), however, this observation is likely to be an outlier genotype relative to the allelic frequencies expected under Hardy–Weinberg Equilibrium for these Mantella populations, rather than the outcome of multiple past hybridisations. These observations are supported by the partition results for progressively larger values of K (i.e., from 2 until 7; Figure S2), which show that from an earlier level on (i.e., K = 2) the localities form well defined clusters, and single individuals appear admixed for the different values of K: the same individual from Besariaka found admixed in the K = 7 solution (Figure 6) was also found admixed in K = 5 and K = 6 solutions (Figure S2); however, different individuals were identified as possibly admixed in the K = 4 and K = 3 solutions. An individual-based clustering approach (without prior locality information) in BAPS rendered incongruent results as no stable clustering solution was identified, further reflecting the effect of the little divergence observed between these localities (e.g., Table 2). Consistent with this hypothesis, a correspondence factor analysis using Genetix  of either the individuals, or groups of individuals by localities, resulted in a large overlap between samples and localities, respectively (Figures S3 and S4).
Pairwise divergence between sampling localities as measured by the FST showed that most of the comparisons were significant after Bonferroni correction (Table 2). The FST values were overall low both within and between species, with the average within species FST (0.05) slightly smaller than the average between species FST (0.07). The observed low differentiation between sampling localities was reflected by the average low frequency of private alleles, namely 0.0819. These values reflect either recent gene flow between populations, or rather little divergence from the ancestral population, also reflected by the low statistical support in the dendrogram depicting the relationships between sampling localities (Figure 7). Migration patterns estimated by BayesAss revealed that the majority of the individuals from one locality can be assigned to its sampling locality (between ca. 65–90%), while the remaining individuals correspond to putative migrants from the other localities, each contributing at most up to 5% of the individuals.
3.6. Demographic History
For each of the populations included in the microsatellite analysis we used Msvar to estimate its past demographic history using as priors a bottleneck, a population expansion and stable demography (Table S6). The three independent runs of Msvar were assessed for convergence using Gelman and Rubin’s statistic that for each parameter of the analysis presented values smaller than 1.1 indicating convergence of the demographic model for each locality independently of the prior used. Overall all populations showed evidence of a drastic demographic bottleneck that reduced effective population sizes in some cases from around 100,000 individuals to sometimes a couple of hundreds (e.g., the population of M. aurantiaca in Andranomandry; Figure 8, Table S6). This bottleneck was estimated in very recent times, with the closest estimate to the present at about 87 years ago and the oldest approximately 600 years ago, and an average estimate of approximately 250 years ago (Figure 8).
4.1. M. aurantiaca, M. crocea and M. milotympanum—Separate Species or Color Morphotypes?
Among the swamp frog species of the M. madagascariensis group, M. aurantiaca is the taxon with the most stable coloration, being characterized by a bright orange color on a translucent base (see Figure 1, Figure 6, Figure 7 and Figure 8). Based on their overall coloration, M. milotympanum and M. aurantiaca were previously thought to be closely related to each other, until genetic data confirmed M. milotympanum to be more closely related to M. crocea from a mitochondrial perspective [22,23]. Based on these results it has been hypothesized that the orange color might have evolved convergently in M. aurantiaca and M. milotympanum . Here we confirm that there are two main mitochondrial lineages of which one is prevalent in M. aurantiaca and the other in M. crocea + M. milotympanum, and we detected widespread mitochondrial haplotype sharing between the southern populations of M. aurantiaca and the M. crocea/M. milotympanum group, in COB (with two haplotypes, C9 and C5, being shared by M. aurantiaca and both brown and green populations of M. crocea), RAG-1 (with widespread haplotype sharing between these taxa) and RAG-2 (with 1 widespread shared haplotype between M. aurantiaca and M. milotympanum) (Figure 3, Figure 4 and Figure 5).
In addition to this widespread haplotype sharing in the analyzed mitochondrial marker, it is worth noting the existence of at least two other cases of mitochondrial introgression, with some COB haplotypes uniquely identified in the southern populations of M. aurantiaca that cluster in the M. crocea/M. milotypmpanum phylogroup (haplotypes C3, C4 and C8), and one haplotype uniquely identified in the Torotorofotsy population of M. aurantiaca that clusters with the M. pulchra/M. madagascariensis (haplotype C11) network (Figure 3). For the first time, we also identify an instance of mitochondrial introgression of M. aurantiaca in M. crocea (green population from Ambohitantely), with haplotype C17 being nested within the M. aurantiaca network, and therefore confirming mitochondrial introgression happening in both directions between these taxa (Figure 3). Mantella aurantiaca was until recently known from 21 localities  and its distribution has been currently updated to include 139 breeding sites . The four analyzed populations of M. aurantiaca come from both its southern (Andranomandry, Sahasarotra and Besariaka) and northern (Torotorofotsy) edge of the known distribution. Two of the southern populations of M. aurantiaca (Andranomandry and Besariaka) formed a single cluster in our microsatellite analysis (Figure 6) suggesting the occurrence of gene flow between these two geographically close populations. Despite this, in our microsatellite-based tree the cluster including the four analyzed populations of M. aurantiaca turned out to be the only grouping that received relevant bootstrap support, therefore indicating that the populations of orange translucent color (i.e., M. aurantiaca) may form a genetically coherent unit also from a nuclear DNA perspective (Figure 6 and Figure 7).
Different from M. aurantiaca, populations of both Mantella crocea and M. milotympanum are characterized by a great chromatic variability. Individuals of M. crocea show a brown-yellowish coloration in the southern and a greenish coloration in the northern part of their range (including the populations of Hevirina, Andaingo, Ambohitantely, and Zahamena analyzed in this study), whereas M. milotympanum populations vary greatly in the extent of the black color around the tympanic area and in the overall coloration (predominantly greenish in the southern and orange in the northern part of their range), despite their restricted distributions. Here we confirm fairly common COB haplotype sharing between M. milotympanum and M. crocea. Although the majority of the analyzed individuals of M. milotympanum have haplotypes that are poorly differentiated, they are different from the haplotypes found in the populations of M. crocea (Figure 3). The analysis of the RAG-1 marker reveal a deeply interconnected M. milotympanum, with only four exclusive haplotypes identified in the analysis of the entire RAG-1 M. milotympanum dataset (Figure 4), although probably this marker suggests an overall widespread incomplete lineage sorting between these taxa.
At the interspecific level, our microsatellite-based population structure analyses recover a single M. milotympanum population (Figure 6), confirming that the variability of their color morphotypes is not indicative of taxonomic differentiation. In contrast, the three populations of M. crocea were recovered as distinct clusters in the BAPS analysis, therefore suggesting that gene flow between color morphs, and more in general between populations, is somewhat limited in this other taxon (Figure 6).
In this study we analyzed the population genetics of three swamp forest poison frogs of Madagascar using a varied panel of molecular markers to infer both recent and more ancient genetic processes of differentiation. We interpret the two main mitochondrial lineages in the three swamp forest species (one prevalent in M. aurantiaca the second in M. crocea/milotympanum; obvious also from our multi-gene phylogenetic analysis; Figure 2), as reflecting an original, probably allopatric divergence. However, we think that inter-lineage gene flow has been a common phenomenon in the subsequent history of these groups, with probably more ancient mitochondrial introgressions (of aurantiaca mtDNA into the Ambohitantely population of M. crocea) and recent gene flow (multiple introgressions between neighbouring populations, of M. crocea mtDNA into M. aurantiaca) (Figure 3). At the same time, we observe that in some cases genetic divergence increases with geographic distance, e.g., in COB the haplotypes of the geographically distant, green-colored populations of M. crocea, are most divergent (Figure 1, Figure 3 and Figure 4). Despite these evidences, it seems that populations are in part maintaining their genetic identity, with admixture playing only a limited role (Figure 6), and recent migration among populations being rare according to the microsatellite data.
4.2. Taxonomic Conclusions
A taxonomic resolution of closely related species complexes such as the M. madagascariensis group is inevitably a highly controversial matter. In the studies available to date ([22,23,31,79]; and this study) different kinds of data sets suggested different relationships within this group. Compared to most other species and candidate species of Madagascar frogs [6,7] the genetic distances between valid species of Mantella spp. are exceedingly low; and this seems to be the case also in the M. madagascariensis group; e.g., in the 16S rRNA, the sequenced specimens of M. crocea and M. milotympanum have identical haplotypes, M. aurantiaca and M. crocea/milotympanum differ by 0.8%, and their distance to M. pulchra and M. madagascariensis is 0.2–1.4%. These values are much lower than the 3% that was set as threshold for initial identification of candidate species for taxonomic revision in neobatrachian frogs [6,80]. Along with the evidence for widespread hybridization between these taxa, this indeed casts doubts on the species status of all these lineages. While we feel that at least some of these deserve being considered as distinct species (e.g., M. aurantiaca), the current data do not allow to conclusively delimit independent evolutionary units in the entire complex, and therefore we do not see a simple answer in immediate reach to clarify the uncertain taxonomic status of M. crocea and M. milotympanum. Given the conflicting genetic results, but also the importance of these species in conservation and their listing in the Convention on the International Trade in Endangered Species (CITES), a certain stability of the (necessarily inconsistent) current taxonomy is important and we therefore here choose not to propose—and strongly discourage—any premature nomenclature and/or taxonomic changes in this group .
Independent from their taxonomic status, our data herein suggest that various populations of the M. madagascariensis group should be managed as separate, genetically differentiated units for conservation. This refers to the taxonomically named units, but also within M. aurantiaca and M. crocea, there are genetically divergent clusters of populations that warrant attention in efforts to maintain the genetic variation of these swamp forest frogs.
While the present paper was under final review, Klonoski et al.  published a population genomic analysis of M. aurantiaca, M. crocea and M. milotympanum using restriction site-associated DNA (RAD) sequencing. Similar to our results, their study found clear genetic clusters that did not align with current species designations, and supported a substantial genetic differentiation between yellow- and green-colored populations of M. crocea. However, the RADseq data —obtained from different individuals and partly different populations than our microsatellite data—found strong indication of admixture at several sites, and an intriguing pattern of genetic clustering of yellow-colored M. crocea specimens with M. aurantiaca, not recovered by our study. These results reinforce our conclusion that the taxonomy of this complex of frogs is convoluted but without an obvious and simple resolution in sight, despite an increasing body of detailed population genetic and population genomic data.
4.3. Mitochondrial Introgression in the Northern Populations of M. madagascariensis?
Samples from the three analyzed populations of M. madagascariensis did not form a single mitochondrial network. The southern populations were extremely homogeneous (only one single haplotype was identified from Ranomafana), whereas, the samples of M. madagascariensis from the central east (morphologically and chromatically indistinguishable from the population from the South) show haplotype sharing and shallow mitochondrial differentiation with the chromatically different M. pulchra (Figure 1 and Figure 3), in apparent agreement with allozyme analyses  and karyological data . On the contrary, the haplotype network of the analyzed RAG-1 fragment show that haplotype sharing between M. pulchra and M. madagascariensis is extremely rare in this gene and only affects a single haplotype (haplotype R65; Figure 4). The observed pattern may reflect either widespread introgression of divergent haplotypes after speciation, masking the true phylogenetic relationships among taxa and populations , or the populations of M. madagascariensis from the North can instead be a different color morph of M. pulchra (considering that the Type locality of M. madagascariensis is “Forêt d’Ambalavatou, entre Mananzarine et Fianarantsoua” it is in fact likely that the population from Ranomafana will be the one maintaining the name M. madagascariensis). The Ranomafana population could not be included in the multi-gene phylogenetic analysis (Figure 2) because no more material (tissue, DNA) from these specimens was left in our collection. Additional sampling of M. madagascariensis, especially from the Ranomafana population, and of M. pulchra is needed to infer the degree of differentiation of these populations and possibly propose a new taxonomic status.
4.4. The Effect of Deforestation on the Demographic History of the Swamp Forest Species
Inferring the demographic history of one species is crucial to better understand how present-day genetic diversity is distributed . A recently published study reveals that both climate change and human colonization triggered habitat loss and fragmentation in Northern Madagascar, two factors that affected the demography of two lemur species . Although it is still widely debated, recent estimates find support for first human activities in Madagascar having started only about 1350 B.P. . Similarly, recent palynological evidence suggests that large scale deforestation in Madagascar occurred later than argued previously, with rises in charcoal abundance older that 2000 y. B.P. being available only from the arid southwest, where natural firing is expected; whereas in wetter regions (such as along the eastern rainforest belt) steep rises in charcoal abundance are mostly found from around 1150–950 y B.P. [86,87]. Similarly interesting, it seems that almost three-quarters of the remaining primary forest of Madagascar was cleared during Colonial Madagascar times (from 1895 to 1925) when population growth was slow and shifting cultivation banned , an imposed change in land use that forced many people to find refuge in the forests and oblige them to survive there as shifting cultivators . Finally, Harper et al.  reports that Madagacar’s ecosystems have been deeply affected by human-mediated habitat loss and fragmentation in the recent decades, with about the 50% remaining forest cover loss since the 1950s [89,90]. Our demographic analyses identified a severe pattern of demographic contraction affecting the three analyzed swamp forest species (Figure 8) and suggest that this contraction took place in the recent past. Although we cannot exclude that our estimates of effective population size may have been inflated by alleles introduced into the populations by migrants from other populations, it is unlikely that this pattern would have led to a concordant reconstruction of bottlenecks in all populations as found by our analyses. Our probability densities are mostly distributed between 1000 and 100 y ago in almost coincidence with early wide-scale deforestation, and maintained at least until 100 y BP, in accordance with other evidences for increased trends in forest loss and fragmentation in more recent times [88,89], suggesting that forest loss was probably a major trigger for the demographic contraction observed in these species. In addition to this, it is known that these frog species are affected by collecting for the pet trade; about 230,000 specimens of Mantella were exported from Madagascar from 1994 to 2003 [29,30], with M. aurantiaca being the third worldwide-most traded amphibian species . However, the onset of large-scale collection for the pet trade is more recent than the estimated demographic bottlenecks in the three swamp frog species, which instead are in temporal coincidence with the massive events of forest loss in the past centuries [86,87,88,89]. While also today, habitat loss is certainly the most important and imminent threat to Mantella populations (and biodiversity in general), pet trade over-collecting can exacerbate this threat and therefore requires careful regulation. Our data exemplify how the massive past forest reduction has drastically reduced the effective population size of amphibian species, and highlight the importance of fine-scale conservation genetic assessments of the remaining fragmented populations of these species.
The following are available online at https://www.mdpi.com/2073-4425/10/4/317/s1, Figure S1: Neighbor-Joining tree reconstructed using the p-distance method and using complete COB dataset; Figure S2: Bayesian analysis of the selected populations of the M. madagascariensis group used for the analysis of population structure. Figure S3: Results of the correspondence factor analysis that shows overlap between localities analyzed by a group-based approach. Figure S4: Results of the correspondence factor analysis that shows overlap between the localities analyzed by an individual-based approach; Table S1: Sample ID, field number, taxon, sampling locality and ID, note on the phenotype and respective COB, RAG-1 and RAG-2 haplotypes and the raw scoring for each of the 9 analyzed microsatellites. Table S2: Taxon, locality names and ID, geographical location and number of individuals sequenced for the COB, RAG-1 and RAG-2 fragments per population. Table S3: List of frequencies of different haplotypes (numbered as in the haplotype networks) per population. Table S4: Primer name and sequences and PCR conditions used in the present study. Table S5: Prior distributions of the parameters estimated in Msvar. Table S6: Highest Posterior Density intervals and mode of posterior distributions of demographic parameters in Msvar.
Conceptualization of the work, A.C., and M.V.; A.C., F.C.E.R. and J.S.H. generate the data, A.C., P.O.-t. and M.V. analyzed and interpreted the data and wrote the first draft of the manuscript; F.C.E.R. and J.S.H. critically revised the manuscript. All authors contributed to drafting the manuscript and have approved the final version of the article.
Financial support was provided by the Volkswagen Foundation and NWO/WOTRO. The Portuguese National Funds through FCT (Foundation for Science and Technology) support the Investigador FCT grant to AC (IF/00209/2014). This work is funded by National Funds through FCT within the IF/00209/2014/CP1256/CT0011 Exploratory Research Project.
We are grateful to numerous colleagues, field and laboratory assistants, and students for contributing to sample collection and laboratory work, in particular Parfait Bora, Devin Edmonds, Euan Edwards, Frank Glaw, Meike Kondermann, Olga Jovanovic, Roger-Daniel Randrianiaina, Fanomezana Ratsoavina, Tokihery Razafindrabe, Fatina Siwczak, Paul Tröger, David R. Vieites, and Katharina C. Wollenberg Valero. This work was carried under collaboration accords between the University of Antananarivo (Department of Animal Biology) and the Technische Universität Braunschweig. We are grateful to the Malagasy authorities for issuing research, export and CITES export permits, and to the Bundesamt für Naturschutz in Germany for issuing CITES import permits.
Conflicts of Interest
The authors declare no conflict of interest.
- Goodman, S.M.; Benstead, J. The Natural History of Madagascar, 1st ed.; The University of Chicago Press: Chicago, IL, USA, 2004; ISBN 0226303071. [Google Scholar]
- Glaw, F.; Vences, M. A Fieldguide to the Amphibians and Reptiles of Madagascar, 3rd ed.; Vences & Glaw Verlag: Cologne, Germany, 2007; ISBN 392944903X. [Google Scholar]
- Crottini, A.; Madsen, O.; Poux, C.; Strauß, A.; Vieites, D.R.; Vences, M. Vertebrate time-tree elucidates the biogeographic pattern of a major biotic change around the K-T boundary in Madagascar. Proc. Natl. Acad. Sci. USA 2012, 109, 5358–5363. [Google Scholar] [CrossRef]
- Wollenberg, K.C.; Vieites, D.R.; Glaw, F.; Vences, M. Speciation in little: The role of range and body size in the diversification of Malagasy mantellid frogs. BMC Evol. Biol. 2011, 11, 217. [Google Scholar] [CrossRef]
- Scherz, M.D.; Vences, M.; Rakotoarison, A.; Andreone, F.; Köhler, J.; Glaw, F.; Crottini, A. Reconciling molecular phylogeny, morphological divergence and classification of Madagascan narrow-mouthed frogs (Amphibia: Microhylidae). Mol. Phylogenet. Evol. 2016, 100, 372–381. [Google Scholar] [CrossRef]
- Vieites, D.R.; Wollenberg, K.C.; Andreone, F.; Köhler, J.; Glaw, F.; Vences, M. Vast underestimation of Madagascar’s biodiversity evidenced by an integrative amphibian inventory. Proc. Natl. Acad. Sci. USA 2009, 106, 8267–8272. [Google Scholar] [CrossRef]
- Perl, R.G.B.; Nagy, Z.T.; Sonet, G.; Glaw, F.; Wollenberg, K.C.; Vences, M. DNA barcoding Madagascar’s amphibian fauna. Amphib.-Reptil. 2014, 35, 197–206. [Google Scholar] [CrossRef]
- Glaw, F.; Vences, M. A Fieldguide to the Amphibians and Reptiles of Madagascar, 2nd ed.; Vences & Glaw Verlag: Cologne, Germany, 1994; ISBN 3929449013. [Google Scholar]
- Köhler, J.; Vieites, D.R.; Bonett, R.M.; García, F.H.; Glaw, F.; Steinke, D.; Vences, M. New amphibians and global conservation: A boost in species discoveries in a highly endangered vertebrate group. BioScience 2005, 55, 693–696. [Google Scholar] [CrossRef]
- Vences, M.; Wollenberg, K.C.; Vieites, D.R.; Lees, D.C. Madagascar as a model region of species diversification. Trends Ecol. Evol. 2009, 24, 456–465. [Google Scholar] [CrossRef]
- Pabijan, M.; Wollenberg, K.C.; Vences, M. Small body size increases the regional differentiation of populations of tropical mantellid frogs (Anura: Mantellidae). J. Evol. Biol. 2012, 25, 2310–2324. [Google Scholar] [CrossRef]
- Wollenberg Valero, K.C.; Garcia-Porta, J.; Rodríguez, A.; Arias, M.; Shah, A.; Randrianiaina, R.D.; Brown, J.L.; Glaw, F.; Amat, F.; Künzel, S.; et al. Transcriptomic and macroevolutionary evidence for phenotypic uncoupling between frog life history phases. Nat. Commun. 2017, 8, 15213. [Google Scholar] [CrossRef]
- Marshall, J.C.; Bastiaans, E.; Caccone, A.; Camargo, A.; Morando, M.; Niemiller, M.L.; Pabijan, M.; Russello, M.A.; Sinervo, B.; Sites, J.W., Jr.; et al. Mechanisms of speciation in reptiles and amphibians: A synopsis. PeerJ Prepr. 2018, 6, e27279v1. [Google Scholar] [CrossRef]
- Orozco-Terwengel, P.; Andreone, F.; Louis, E., Jr.; Vences, M. Mitochondrial introgressive hybridization following a demographic expansion in the tomato frogs of Madagascar, genus Dyscophus. Mol. Ecol. 2013, 22, 6074–6090. [Google Scholar] [CrossRef]
- Crottini, A.; Chiari, Y.; Mercurio, V.; Meyer, A.; Vences, M.; Andreone, F. Into the canyons: The phylogeography of the Malagasy frogs Mantella expectata and Scaphiophryne gottlebei in the arid Isalo Massif, and its significance for conservation (Amphibia: Mantellidae and Microhylidae). Org. Divers. Evol. 2008, 8, 368–377. [Google Scholar] [CrossRef]
- Vieites, D.R.; Chiari, Y.; Vences, M.; Andreone, F.; Rabemananjara, F.C.E.; Bora, P.; Nieto-Román, S.; Meyer, A. Mitochondrial evidence for distinct phylogeographic units in the endangered Malagasy poison frog Mantella bernhardi. Mol. Ecol. 2006, 15, 1617–1625. [Google Scholar] [CrossRef]
- Rabemananjara, R.C.E.; Chiari, Y.; Ravoahangimalala Ramilijaona, O.; Vences, M. Evidence for recent gene flow between north-eastern and south-eastern Madagascan poison frogs from a phylogeography of the Mantella cowani group. Front. Zool. 2007, 4, 1. [Google Scholar] [CrossRef] [PubMed]
- Avise, J.C. Phylogeography: The History and Formation of Species; Harvard University Press: Cambridge, MA, USA, 2000; ISBN 9780674666382. [Google Scholar]
- Vences, M.; Wake, D.B. Speciation, species boundaries and phylogeography of amphibians. In Amphibian Biology; Heatwole, H.H., Tyler, M., Eds.; Surrey Beatty & Sons: Chipping Norton, Australia, 2007; Volume 7, ISBN 9780980311310. [Google Scholar]
- Zeisset, I.; Beebee, T.J.C. Amphibian phylogeography: A model for understanding historical aspects of species distributions. Heredity 2008, 101, 109–119. [Google Scholar] [CrossRef] [PubMed]
- Crottini, A.; Brown, J.L.; Mercurio, V.; Glaw, F.; Vences, M.; Andreone, F. Phylogeography of the poison frog Mantella viridis (Amphibia: Mantellidae) reveals chromatic and genetic differentiation across ecotones in northern Madagascar. J. Zool. Syst. Evol. Res. 2012, 50, 305–314. [Google Scholar] [CrossRef]
- Chiari, Y.; Vences, M.; Vieites, D.R.; Rabemananjara, F.C.E.; Bora, P.; Ramilijaona Ravoahangimalala, O.; Meyer, A. New evidence for parallel evolution of colour patterns in Malagasy poison frogs (Mantella). Mol. Ecol. 2004, 13, 3763–3774. [Google Scholar] [CrossRef] [PubMed]
- Vences, M.; Chiari, Y.; Raharivololoniaina, L.; Meyer, A. High mitochondrial diversity within and among populations of Malagasy poison frogs. Mol. Phylogenet. Evol. 2004, 30, 295–307. [Google Scholar] [CrossRef]
- Rabemananjara, F.C.E.; Crottini, A.; Chiari, Y.; Anderone, F.; Glaw, F.; Duguet, R.; Bora, P.; Ravoahangimalala Ramilijaona, O.; Vences, M. Molecular systematics of Malagasy poison frogs in the Mantella betsileo and M. laevigata species groups. Zootaxa 2007, 1501, 31–44. [Google Scholar] [CrossRef]
- Cocca, W.; Rosa, G.M.; Andreone, F.; Aprea, G.; Eusebio Bergò, P.; Mattioli, F.; Mercurio, V.; Randrianirina, J.E.; Rosado, D.; Vences, M.; et al. The herpetofauna (Amphibia, Crocodylia, Squamata, Testudines) of the Isalo Massif, Southwest Madagascar: Combining morphological, molecular and museum data. Salamandra 2018, 54, 178–200. [Google Scholar]
- Daly, J.W.; Highet, R.J.; Myers, C.W. Occurrence of skin alkaloids in non-dendrobatid frogs from Brazil (Bufonidae), Australia (Myobatrachidae) and Madagascar (Mantellinae). Toxicon 1984, 22, 905–919. [Google Scholar] [CrossRef]
- Vences, M.; Glaw, F.; Böhme, W. Evolutionary correlates of microphagy in alkaloid-containing frogs (Amphibia: Anura). Zoologischer Anzeiger 1998, 236, 217–230. [Google Scholar]
- Tarvin, R.D.; Santos, J.C.; O’Connell, L.A.; Zakon, H.H.; Cannatella, D.C. Convergent substitutions in a sodium channel suggest multiple origins of toxin resistance in poison frogs. Mol. Biol. Evol. 2016, 33, 1068–1081. [Google Scholar] [CrossRef]
- Rabemananjara, F.C.E.; Ravoahangimalala Ramilijaona, O.; Raminosoa, N.R.; Andreone, F.; Bora, P.; Carpenter, A.I.; Glaw, F.; Razafindrabe, T.; Vallan, D.; Vieites, D.R.; et al. Malagasy poison frogs in the pet trade: A survey of levels of exploitation in the species in the genus Mantella. In A Conservation Strategy for the Amphibians of Madagascar; Andreone, F., Randriamahazo, H., Eds.; Museo Regionale di Scienze Naturali: Torino, Italy, 2008; pp. 277–300. ISBN 9788886041720. [Google Scholar]
- Carpenter, A.I.; Andreone, F.; Moore, R.D.; Griffiths, R.A. A review of the international trade in amphibians:the types, levels and dynamics of trade in CITES-listed species. Oryx 2014, 48, 565–574. [Google Scholar] [CrossRef]
- Vences, M.; Hille, A.; Glaw, F. Allozyme differentiation in the genus Mantella (Amphibia: Anura: Mantellinae). Folia Zool. 1998, 47, 261–274. [Google Scholar]
- Rabemananjara, F.C.E.; Bora, P.; Razafindrabe, T.; Randriamitso, E.; Ravoahangimalala Ramilijaona, O.; Raminosoa, N.R.; Rakotondravony, D.; Vieites, D.R.; Vences, M. Rapid assessments of population sizes in ten species of Malagasy poison frogs, genus Mantella. In A Conservation Strategy for the Amphibians of Madagascar; Andreone, F., Randriamahazo, H., Eds.; Museo Regionale di Scienze Naturali: Torino, Italy, 2008; ISBN 9788886041720. [Google Scholar]
- Tessa, G.; Mattioli, F.; Mercurio, V.; Andreone, F. Egg numbers and fecundity traits in nine species of Mantella poison frogs from arid grasslands and rainforests of Madagascar (Anura: Mantellidae). Madag. Conserv. Dev. 2009, 4, 113–119. [Google Scholar] [CrossRef]
- Vieites, D.R.; Rabemananjara, F.E.C.; Bora, P.; Razafimahatratra, B.; Ramilijaona Ravoahangimalala, O.; Vences, M. Distribution and population density of the black-eared Malagasy poison frog, Mantella milotympanum Staniszewski, 1996 (Amphibia: Mantellidae). In African Biodiversity: Molecules, Organisms, Ecosystems, Proceedings of the 5th International Symposium on Tropical Biology, Museum Koenig, Bonn, Germaby, 5–8 May 2004; Huber, B.A., Sinclair, B.J., Lampe, K.H., Eds.; Springer: Berlin, Germany, 2005; pp. 197–204. ISBN 978-0-387-24320-7. [Google Scholar]
- Bora, P.; Dolch, R.; Jenkins, R.; Jovanovic, O.; Rabemananjara, F.C.E.; Randrianirina, J.; Rafanomezantsoa, J.; Raharivololoniaina, L.; Ramilijaona, O.; Raminosoa, N.; et al. Geographical distribution of three species of Malagasy poison frogs of high conservation priority: Mantella aurantiaca, M. crocea and M. milotympanum. Herpetol. Notes 2008, 1, 39–48. [Google Scholar]
- Randrianavelona, R.; Rakotonoely, H.; Ratsimbazafy, J.; Jenkins, R.K.B. Conservation assessment of the critically endangered frog Mantella aurantiaca in Madagascar. Afr. J. Herpetol. 2010, 59, 65–78. [Google Scholar] [CrossRef]
- Piludu, N.; Dubos, N.; Razafimanahaka, J.H.; Razafindraibe, P.; Christian, J.; Jenkins, R.K.B. Distribution, threats and conservation of a Critically Endangered amphibian (Mantella aurantiaca) in Eastern Madagascar. Herpetol. Notes 2015, 8, 119–123. [Google Scholar]
- Schaefer, H.C.; Vences, M.; Veith, M. Molecular phylogeny of Malagasy poison frogs, genus Mantella (Anura: Mantellidae): Homoplastic evolution of color pattern in aposematic amphibians. Org. Divers. Evol. 2002, 2, 97–105. [Google Scholar] [CrossRef]
- Andriamaharavo, N.R.; Garraffo, H.M.; Spande, T.F.; Giddings, L.A.; Vieites, D.R.; Vences, M.; Saporito, R.A. Individual and geographic variation of skin alkaloids in three swamp-forest species of Madagascan poison frogs (Mantella). J. Chem. Ecol. 2015, 41, 837–847. [Google Scholar] [CrossRef]
- Woodhead, C.; Vences, M.; Vieites, D.R.; Gamboni, I.; Fisher, B.L.; Griffiths, R.A. Specialist or generalist? Feeding ecology of the Malagasy poison frog Mantella aurantiaca. Herpetol. J. 2007, 17, 225–236. [Google Scholar]
- International Union for Conservation of Nature. The IUCN Red List of Threatened Species. Version 2018-2. IUCN, Gland, Switzerland. Available online: http://www.iucnredlist.org (accessed on 1 February 2019).
- Douglass, K.; Zinke, J. Forging ahead by land and by sea: Archaeology and Paleoclimate Reconstruction in Madagascar. Afr. Archaeol. Rev. 2015, 32, 267–299. [Google Scholar] [CrossRef]
- Anderson, A.; Clark, G.; Haberle, S.; Higham, T.; Nowak-Kemp, M.; Prendergast, A.; Radimilahy, C.; Rakotozafy, L.M.; Ramilisonina; Schwenninger, J.-L.; et al. New evidence of megafaunal bone damage indicates late colonization of Madagascar. PLoS ONE 2018, 13, e0204368. [Google Scholar] [CrossRef] [PubMed]
- Hauswaldt, S.J.; Crottini, A.; Vences, M. Nine microsatellite loci for the Malagasy Poison Frogs (Mantella). Conserv. Genet. Resour. 2009, 1, 269–271. [Google Scholar] [CrossRef]
- Bruford, M.W.; Hanotte, O.; Brookfield, J.F.Y.; Burke, T. Single locus and multilocus DNA fingerprint. In Molecular Genetic Analysis of Populations: A Practical Approach; Hoelzel, A.R., Ed.; IRL Press: Oxford, UK, 1992; pp. 225–270. [Google Scholar]
- Kumar, S.; Stecher, G.; Tamura, K. MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 2016, 33, 1870–1874. [Google Scholar] [CrossRef] [PubMed]
- Stephens, M.; Smith, N.J.; Donnelly, P. A new statistical method for haplotype reconstruction from population data. Am. J. Hum. Genet. 2001, 68, 978–989. [Google Scholar] [CrossRef] [PubMed]
- Rozas, J.; Ferrer-Mata, A.; Sánchez-DelBarrio, J.C.; Guirao-Rico, S.; Librado, P.; Ramos-Onsins, S.E.; Sánchez-Gracia, A. DnaSP 6: DNA sequence polymorphism analysis of large datasets. Mol. Biol. Evol. 2017, 34, 3299–3302. [Google Scholar] [CrossRef] [PubMed]
- Tajima, F. The effect of change in population size on DNA polymorphism. Genetics 1989, 123, 597–601. [Google Scholar]
- Nei, M. Molecular Evolutionary Genetics, 1st ed.; Columbia University Press: New York, NY, USA, 1987; ISBN 9780231063210. [Google Scholar]
- Tajima, F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 1989, 123, 585–595. [Google Scholar]
- Darriba, D.; Taboada, G.L.; Doallo, R.; Posada, D. jModelTest 2: More models, new heuristics and parallel computing. Nat. Methods 2012, 9, 772. [Google Scholar] [CrossRef]
- Lanfear, R.; Calcott, B.; Ho, S.Y.W.; Guindon, S. PartitionFinder: Combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 2012, 29, 1695–1701. [Google Scholar] [CrossRef]
- Ronquist, F.; Teslenko, M.; Van der Mark, P.; Ayres, D.L.; Darling, A.; Höhna, S.; Larget, B.; Liu, L.; Suchard, M.A.; Huelsenbeck, J.P. MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 2012, 61, 539–542. [Google Scholar] [CrossRef]
- Bruen, T.C.; Philippe, H.; Bryant, D. A simple and robust statistical test for detecting the presence of recombination. Genetics 2006, 172, 2665–2681. [Google Scholar] [CrossRef]
- Huson, D.H.; Bryant, D. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 2006, 23, 254–267. [Google Scholar] [CrossRef]
- Clement, M.; Posada, D.; Crandall, K.A. TCS: A computer program to estimate gene genealogies. Mol. Ecol. 2000, 9, 1657–1659. [Google Scholar] [CrossRef]
- Templeton, A.R.; Crandall, K.A.; Sing, C.F. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 1992, 132, 619–633. [Google Scholar]
- Posada, D.; Crandall, K.A. Intraspecific phylogenetics: Trees grafting into networks. Trends Ecol. Evol. 2001, 16, 37–45. [Google Scholar] [CrossRef]
- Pfenninger, M.; Posada, D. Phylogeographic history of the land snail Candidula unifasciata (Poiret 1801) (Helicellinae, Stylommatophora): Fragmentation, corridor migration and secondary contact. Evolution 2002, 56, 1776–1788. [Google Scholar] [CrossRef]
- Van Oosterhout, C.; Hutchinson, W.F.; Shipley, P.; Wills, D.P.M. MICRO-CHECKER: Software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 2004, 4, 535–538. [Google Scholar] [CrossRef]
- Excoffier, L.; Lischer, H.E.L. Arlequin Suite ver 3.5, a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 2010, 10, 564–567. [Google Scholar] [CrossRef] [PubMed]
- Dieringer, D.; Schloetter, C. Microsatellite analyser (MSA): A platform independent analysis tool for large microsatellite data sets. Mol. Ecol. Notes 2003, 3, 167–169. [Google Scholar] [CrossRef]
- Raymond, M.; Rousset, F. Genepop (Version 1.2): Population-genetics software for exact tests and ecumenicism. J. Hered. 1995, 86, 248–249. [Google Scholar] [CrossRef]
- Corander, J.; Marttinen, P. Bayesian identification of admixture events using multilocus molecular markers. Mol. Ecol. 2006, 15, 2833–2843. [Google Scholar] [CrossRef]
- Felsenstein, J. PHYLIP (Phylogeny Inference Package), version 3.61; Department of Genome Sciences, University of Washington; Distributed by the author; Seattle, WA, USA, 1993. [Google Scholar]
- Bowcock, A.M.; Ruiz-Linares, A.; Tomfohrde, J.; Minch, E.; Kidd, J.R.; Cavalli-Sforza, L.L. High resolution of human evolutionary trees with polymorphic microsatellites despite a constraint in allele length. Nature 1994, 368, 455–457. [Google Scholar] [CrossRef]
- Wilson, G.A.; Rannala, B. Bayesian inference of recent migration rates using multilocus genotypes. Genetics 2003, 163, 1177–1191. [Google Scholar]
- Suchard, M.A.; Lemey, P.; Baele, G.; Ayres, D.L.; Drummond, A.J.; Rambaut, A. Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 2018, 4, vey016. [Google Scholar] [CrossRef]
- Beaumont, M.A. Detecting population expansion and decline using microsatellites. Genetics 1999, 153, 2013–2029. [Google Scholar]
- Storz, J.F.; Beaumont, M.A. Testing for genetic evidence of population expansion and contraction: An empirical analysis of microsatellite DNA variation using a hierarchical Bayesian model. Evolution 2002, 56, 154–166. [Google Scholar] [CrossRef]
- Brinkmann, B.; Klintschar, M.; Neuhuber, F.; Huhne, J.; Rolf, B. Mutation rate in human microsatellites: Influence of the structure and length of the tandem repeat. Am. J. Hum. Genet. 1998, 62, 1408–1415. [Google Scholar] [CrossRef]
- Brohede, J.; Primmer, C.R.; Moller, A.; Ellegren, H. Heterogeneity in the rate and pattern of germline mutation at individual microsatellite loci. Nucleic Acids Res. 2002, 30, 1997–2003. [Google Scholar] [CrossRef] [PubMed]
- Bulut, Z.; McCormick, C.R.; Gopurenko, D.; Williams, R.N.; Bos, D.H.; DeWoody, J.A. Microsatellite mutation rates in the eastern tiger salamander (Ambystoma tigrinum tigrinum) differ 10-fold across loci. Genetica 2009, 136, 501–504. [Google Scholar] [CrossRef] [PubMed]
- Brooks, S.P.; Gelman, A. General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 1998, 7, 434–455. [Google Scholar]
- Plummer, M.; Best, N.; Cowles, K.; Vines, K. CODA: Convergence diagnosis and output analysis for MCMC. R News 2006, 6, 7–11. [Google Scholar]
- R Development Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2009. [Google Scholar]
- Belkhir, K.; Borsa, P.; Chikhi, L.; Raufaste, N.; Bonhomme, F. GENETIX4.05, Logiciel sous Windows TM pour la génétique des Populations. Laboratoire génome, Populations, Interactions, CNRS UMR. 5000; Université de Montpellier II: Montpellier, France, 2004. [Google Scholar]
- Odierna, G.; Vences, M.; Aprea, G.; Lötters, S.; Andreone, F. Chromosome data for Malagasy poison frogs (Amphibia: Ranidae: Mantella) and their bearing on taxonomy and phylogeny. Zool. Sci. 2001, 18, 505–514. [Google Scholar] [CrossRef]
- Fouquet, A.; Gilles, A.; Vences, M.; Marty, C.; Blanc, M.; Gemmell, N.J. Underestimation of species richness in Neotropical frogs revealed by mtDNA analyses. PLoS ONE 2007, 2, e1109. [Google Scholar] [CrossRef]
- Pinna, P.H.; Fernandes, D.S.; Passos, P. If you choose not to decide you still have made a choice. Bionomina 2018, 13, 65–68. [Google Scholar] [CrossRef]
- Klonoski, K.; Bi, K.; Rosenblum, E.B. Phenotypic and genetic diversity in aposematic Malagasy poison frogs (genus Mantella). Ecol. Evol. 2019, 9, 2725–2742. [Google Scholar] [CrossRef]
- Albertson, R.C.; Markert, J.A.; Danley, P.D.; Kocher, T.D. Phylogeny of a rapidly evolving clade: The cichlid fishes of Lake Malawi, East Africa. Proc. Natl. Acad. Sci. USA 1999, 96, 5107–5110. [Google Scholar] [CrossRef]
- Keinan, A.; Clark, A.G. Recent explosive human population growth has resulted in an excess of rare genetic variants. Science 2012, 336, 740–743. [Google Scholar] [CrossRef]
- Salmona, J.; Heller, R.; Quéméré, E.; Chikhi, L. Climate change and human colonization triggered habitat loss and fragmentation in Madagascar. Mol. Ecol. 2017, 26, 5203–5222. [Google Scholar] [CrossRef] [PubMed]
- Burney, D.A.; Burney, L.P.; Godfrey, L.R.; Jungers, W.L.; Goodman, S.M.; Wright, H.T.; Jull, A.J.T. A chronology for late prehistoric Madagascar. J. Hum. Evol. 2004, 47, 25–63. [Google Scholar] [CrossRef]
- Burns, S.J.; Godfrey, L.R.; Faina, P.; McGee, D.; Hardt, B.; Ranivoharimanana, L.; Randrianasy, J. Rapid human-induced landscape transformation in Madagascar at the end of the first millennium of the Common Era. Quat. Sci. Rev. 2016, 134, 92–99. [Google Scholar] [CrossRef]
- Jarosz, L. Defining and explaining tropical deforestation: Shifting cultivation and population growth in colonial Madagascar (1896–1940). Econ. Geogr. 1993, 69, 366–379. [Google Scholar] [CrossRef]
- Harper, G.J.; Steininger, M.K.; Tucker, C.J.; Juhn, D.; Hawkins, F. Fifty years of deforestation and forest fragmentation in Madagascar. Environ. Conserv. 2007, 34, 325–333. [Google Scholar] [CrossRef]
- Schwitzer, C.; Chikhi, L.; Donati, G.; Irwin, M.; Johnson, S.E.; Mittermeier, R.A.; Peacock, H.; Ratsimbazafy, J.; Razafindramanana, J.; Louis, E.E., Jr.; et al. Protecting lemurs–response. Science 2014, 344, 358–360. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Map with sampling localities and cytochrome b (COB) haplotypic frequencies. Original layers from the background map from the Madagascar Vegetation Mapping project: orange, reddish and green colors mark remaining primary vegetation.
Figure 2. Phylogenetic trees of the genus Mantella, showing the position of the target species in the Mantella madagascariensis group, obtained by Bayesian Inference analysis of (a) 2531 bp of four mitochondrial genes (COB, cytochrome oxidase subunit 1 (COX1), 12S, 16S), and (b) 1710 bp of three nuclear genes (recombination activating gene 1 (RAG-1), recombination activating gene 2 (RAG-2), brain-derived neurotrophic factor (BDNF)). Values at nodes are Bayesian posterior probabilities (not shown if <0.9).
Figure 3. Haplotype network of 49 unique COB haplotypes of Mantella crocea (green; darker green highlights the haplotypes found in the populations that show the green phenotype), M. milotympanum (orange), M. aurantica (red), M. madagascariensis (blue) and M. pulchra (light blue), based on the analysis of a 552-bp DNA fragment of this mitochondrial gene. Size of circles is proportional to the number of individuals sharing a given haplotype (haplotype frequencies as in the inset legend). Shared haplotypes are indicated by circles, with the area being proportional to the number of individuals sharing that haplotype. Inferred intermediate haplotypes are indicated by small black dot.
Figure 4. Haplotype network of 78 unique haplotypes of Mantella crocea, M. milotympanum, M. aurantica, M. madagascariensis and M. pulchra (color code as in Figure 1 and Figure 3), based on the analysis of a 732-bp fragment of the RAG-1 gene. Size of circles is proportional to the number of individuals sharing a given haplotype. Shared haplotypes are indicated by circles, with the area being proportional to the number of individuals sharing that haplotype (haplotype frequencies as in the inset legend). Inferred intermediate haplotypes are indicated by small black dots.
Figure 5. Haplotype network of 17 unique haplotypes of Mantella crocea, M. milotympanum, M. aurantica, M. madagascariensis, and M. pulchra (color code as in Figure 1, Figure 3 and Figure 4), based on the analysis of a 531-bp fragment of the RAG-2 gene. Size of circles is proportional to the number of individuals sharing a given haplotype. Shared haplotypes are indicated by circles, with the area being proportional to the number of individuals sharing that haplotype (haplotype frequencies as in the inset legend). Inferred intermediate haplotypes are indicated by small black dots.
Figure 6. Bayesian analysis of the selected populations of the M. madagascariensis group used for the analysis of population structure. Barplot of the clustering solution with the highest posterior probability (K = 7). Each individual in the analysis is represented by a vertical bar, and groups of individuals of the same sampling site are grouped in boxes with an associated label below. The only significantly admixed individual (from the Besariaka population of M. aurantiaca) presents its genotype partitioned in colors according to the proportion of it that resembles each of the clusters (e.g., 32% brown reflecting that 32% of the genotype matches that of the M. crocea cluster of Ampangadimbolana).
Figure 7. Dendrogram of the sampling localities used for the population structure analyses with nine microsatellite loci. Only the four analyzed populations of M. aurantiaca form a substantially supported cluster (bootstrap support 91%).
Figure 8. Demographic history of the analyzed Mantella spp. populations. For each sampling locality the estimated current effective population size (solid lines; Ne graph), the ancestral effective population size (dashed lines; Ne graph) and the time of the bottleneck (dashed lines; t graph) are shown. Each population presents three curves corresponding to the posterior distribution using three different priors reflecting a bottleneck, expansion and stable demographic scenarios. Each sampling locality is represented with a different color (name color key is shown in the t graphs). Ne, effective population size; t, time in years.
Table 1. Distribution of sample size per site and species for the population genetic analyses. n: sample size; ANAP: average number of alleles per locus; HO: observed heterozygosity; HE: expected heterozygosity; FIS: inbreeding coefficient.
Table 2. Pairwise FST between sampling localities. Fields in white correspond to within species comparisons, and fields in grey to between species comparisons. FST values can be found above the diagonal, Bonferroni corrected p-values are provided below the diagonal with non-significant p-values (i.e., p-values > 0.05) indicated as n.s. aur—M. aurantiaca; croc—M. crocea; mil—M. milotympanum.
© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).