Phylogenomic species delimitation in microendemic frogs of the Brazilian Atlantic Forest

The advent of next-generation sequencing allows researchers to use large-scale datasets for species delimitation analyses, yet one can envision an inflection point where the added accuracy of including more loci does not offset the increased computational burden. One alternative to including all loci could be to prioritize the analysis of loci for which there is an expectation of high informativeness, such as those with higher numbers of parsimony-informative sites. Here, we explore the issue of species delimitation and locus selection with species from two anuran genera: Melanophryniscus (Bufonidae) and Brachycephalus (Brachycephalidae). Montane species in these genera have been isolated in sky islands across the southern Brazilian Atlantic Forest, which led to the formation of a number of microendemic species. To delimit species, we obtained genetic data using target enrichment of ultraconserved elements from 32 populations (13 for Melanophryniscus and 19 for Brachycephalus), and we were able to create datasets that included over 800 loci with no missing data. We ranked loci according to their corresponding number of parsimony-informative sites, and we performed species delimitation analyses using BPP in each genus based on the top 10, 20, 40, 80, 160, 320, and 640 loci. We also conducted several additional analyses using 10 randomly sampled datasets containing the same numbers of loci to discriminate the relative contribution of increasing the number of loci from prioritizing those with higher informativeness. We identified three types of node: nodes with either consistently high or low support regardless of the number of loci or their informativeness, and nodes that were initially poorly supported, but their support became stronger with more data. Adding more loci had a stronger impact on model support than prioritizing loci for their informativeness, but this effect was less apparent in datasets with more than 160 loci. When viewed across all sensitivity analyses, our results suggest that the current species richness in both genera might have been underestimated. In addition, our results provide useful guidelines to the use of different sampling strategies to carry out species delimitation with phylogenomic datasets.

alone (Bornschein et al. 2015(Bornschein et al. , 2016bPie and Ribeiro 2015;Ribeiro et al. 2015), which can lead to underestimates of species diversity (Bickford et al. 2007). A recent delimitation study using a handful of loci suggested that there could be a number of Brachycephalus and Melanophryniscus species that remain undescribed, although variability between loci in delimitation schemes (Firkowski et al. 2016) suggested that more extensive data were needed to firmly establish species boundaries.
The main goal of this study is to delimit species of Brachycephalus and Melanophryniscus under the MSC using ultraconserved elements (UCEs, Faircloth et al. 2012). Specifically, we used a Bayesian Markov chain Monte Carlo program for species delimitation (BPP 3.3;Yang 2015) that relies on the MSC model to compare different models of species delimitation while accounting for incomplete lineage sorting due to ancestral polymorphism and gene tree-species tree conflicts Rannala 2010, 2014;Rannala and Yang 2013). The lineages investigated in this study are particularly suitable for the use of BPP because the candidate Pie et al. page: 7 species are almost always allopatric, such that gene flow among populations is likely low to nonexistent (simulation studies suggest that the migration threshold for genetic isolation is about 1 individual per 10 generations [Zhang et al. 2011]). In addition to investigating species limits in Brachycephalus and Melanophryniscus, we also explored the importance of locus variation (informativeness) and locus number on the ideal species delimitation schemes.

METHODS
We obtained tissue samples from field-collected specimens of 13 populations of Melanophryniscus and 19 populations of Brachycephalus (Table 1) tridactylus. Brachycephalus sulfuratus does not belong to the pernix species group, but was included in the analyses to improve the rooting of the guide tree. Also, these samples include all known records of montane Melanophryniscus, except those ascribed to M. biancae and M. vilavelhensis, including several novel records in relation to Bornschein et al. (2015). We extracted genomic DNA using PureLink Genomic DNA kit (Invitrogen, USA) and we fragmented the obtained DNA using a BioRuptor NGS (Diagenode). We prepared Illumina libraries using KAPA library preparation kits (Kapa Biosystems) and custom sequence tags unique to each sample (Faircloth and Glenn 2012). To enrich targeted UCE loci, we followed an established workflow (Gnirke et al. 2009;Blumenstiel et al. 2010)  We filtered reads for adapter contamination, low-quality ends, and ambiguous bases using an automated pipeline (https://github.com/faircloth-lab/illumiprocessor) that incorporates Trimmomatic (Bolger et al. 2014). We assembled reads for each taxon using Trinity (Grabehrr et al. 2011). We used the PHYLUCE software package (Faircloth 2015) to align assembled contigs back to their associated UCE loci, remove duplicate matches, and create a taxon-specific database of contig-to-UCE matches. We then generated two alignments: all Brachycephalus samples using Melanophryniscus alipioi as their outgroup, and all Melanophryniscus species using Brachycephalus sulfuratus as their outgroup. We selected loci to create 100% complete data sets for both genera, leading to 820 loci in the Brachycephalus data set and 1227 loci in the Melanophryniscus data set. We aligned data for each individual in each data set using MAFFT (Katoh 2013), and we trimmed resulting alignments using GBlocks (Castresana 2000) with default parameters.
Ideally one would determine the informativeness of a given locus based on how much information it contributes to a given analysis (e.g. Townsend 2007;Gilbert et al. 2015).
However, that would involve first running the entire analysis to then determine their informativeness a posteriori. A simple alternative is to calculate the absolute number of Pie et al. page: 9 parsimony-informative sites (PIS -a single nucleotide polymorphism that is present in more than one individual) for each locus and to use this measure as a proxy for informativeness during species delimitation analyses, as recently proposed by Hime et al. (2016). Other studies used PISs as a measure of informativeness in the context of phylogenetic inference as well (e.g. Hosner et al. 2016;Meiklejohn et al. 2016). We calculated the number of PIS for each locus using IPS 0.0-7 (Heibl 2014) and compared them between datasets using a Spearman's test with the COR.TEST function in R 3.3.2 (R Core Team 2017). We then ranked loci in each dataset according to their corresponding number of PIS, and we created subsets of the top 10, 20, 40, 80, 160, 320, and 640 most informative loci for subsequent analyses.
We used BPP 3.3 (Yang 2015) in all species delimitation analyses. Important assumptions of the MSC model implemented in BPP are no recombination within a locus, free recombination between loci, no migration (gene flow) between species and neutral evolution, and all of which seem reasonable given our model system of allopatric populations assessed using UCEs. It is also important to note that, although the core of UCEs tend to be highly conserved due to strong selection, flanking regions tend to behave similarly to neutral DNA from introns (see Jarvis et al. 2014). Every population was considered initially as a potential species as a conservative starting point (see Olave et al. 2014), allowing for the possibility of later being lumped together by BPP depending on the analyzed data. Three different sets of priors for population size (θ) and divergence time at the root of the species tree (τ 0 ) being assigned gamma priors were used: (1) small ancestral population sizes and shallow interspecific divergences: θ ~ Γ(2, 1000), τ 0~ Γ(2, 1000); (2) large ancestral population sizes and shallow interspecific divergences: θ ~ Γ(2, 1000), τ 0~ Γ(1, 10); and (3) large ancestral population sizes and deep interspecific divergences: θ ~ Γ(1, 10), τ 0~ Γ(1, 10). These sets of priors were used to test their effect on the chosen species Pie et al. page: 10 delimitation scheme, but given what is known about the microendemic distribution and recent divergence times of the studied species, the first set of priors seems more realistic in this case and will be given priority when conflicting results are detected. Other divergence time parameters are assigned the Dirichlet prior (Yang and Rannala 2010: equation 2). We used the A10 model (speciesdelimitation=1, speciestree=0) for species delimitation using a user-specified guide tree (Yang and Rannala 2010;Rannala and Yang 2013). For computational reasons, we assumed that all loci have the same mutation rate ("locusrate=0"). We used the uniform rooted tree prior (speciesmodelprior=1) and the same θs were assumed across loci (heredity=0). We used the option of automatic fine-tuning the MCMC searches, which provided adequate acceptance proportions (i.e. between 20 and 80%). Analyses were repeated using rjMCMC algorithm 0 (ε=2) and algorithm 1 (α=2, m=1), as well as including/excluding sites with gaps (cleandata=0,1), given that the latter has been recently shown as potentially affecting species delimitation using phylogenomic data (Domingos et al. 2017). Each analysis was run at least twice to confirm consistency between runs, given that it seems to be the most effective method to guarantee the reliability of the results (Yang 2015). Given the number of tips on the guide tree, in each analysis there were 768 potential species delimitation models for Brachycephalus and 145 for Melanophryniscus. In total, we ran 336 analyses (2 genera X 7 datasets with varying numbers of loci X 2 algorithms X 3 sets of priors X 2 replicates X 2 treatments of gaps). Each analysis was run for 10,000 generations, sampling every 10 th generation and with the first 10% used as burnin. Although we are aware that BPP can carry out simultaneous species tree and species delimitation analyses (i.e. the A11 model), preliminary analyses using a variety of phylogenetic inference methods showed consistent topologies, such that the computational demands of coestimating species trees (e.g. Caviedes-Solis et al. 2015) would not compensate for the Pie et al. page: 11 substantial increase in run time. Moreover, BPP has been shown to be robust to errors in guide trees (Zhang et al. 2014;Caviedes-Solis et al. 2015).
Finally, in order to discriminate the relative contribution of increasing the number of loci from prioritizing those with higher informativeness, 10 randomly-sampled datasets with the same numbers of loci were analyzed under model A10, algorithm 0, cleandata=0, and priors θ ~ Γ(2, 1000), τ 0~ Γ(2, 1000), and the obtained results were contrasted with the corresponding previous analysis using ranked loci.

RESULTS
We found considerable variation among UCE loci in their informativeness, as indicated by the distribution of the number of PISs for each locus ( Figure S1). In particular, most loci have few PIS and thus are relatively less able to provide information about intra-and interspecific genetic variation. Interestingly, when the same UCE locus was compared between the two datasets, we found a significant correlation in the number of PISs (ρ=0.31, p=1.88e-15), meaning that an informative locus from one genus tended to be informative in the other genus as well.
This relationship was skewed, such that a poorly informative locus from one genus was usually also poorly informative in the other genus; however, a highly informative locus from one genus might not be equally informative in the other (Figure 1). Nevertheless, these results suggest that one could, in principle, select UCE loci that are most informative based on similar datasets involving lineages that are not very closely related (e.g. within Anura). usually those closest to the base of the guide trees (2) nodes with low (or at least highly inconsistent) support across analyses, which were more often found in some regions near the tips of the guide trees; finally, (3) nodes ranging from low to very high support, which are indicated in Figures 2 and 3 with black circles. In the latter case, a clear pattern emerges of a monotonic increase in support with more loci (Figure 4). There seems to be a general monotonic increase in support with more data, which is expected in statistically consistent methods, although there was a slight decrease in the posterior probabilities between 20 and 40 for some analyses (e.g. Figure 4 E, F, L). Also, excluding sites with gaps (clear=1) does not seem to provide important insight, given that they simply led to more unstable posterior probability estimates (possibly because problematic sites had already been removed using Gblocks prior to the species delimitation).
A comparison of the results from ranked and randomly selected loci is shown in Figure 5. To facilitate their comparison, we present only the mean PPs between replicates for the ranked loci, and the mean PPs among the 10 randomly selected loci datasets. In general, the effect of adding more loci had a stronger impact on the variation in PPs than ranking them, given that the resulting PPs were consistently closer to the results of the 640-locus dataset. This effect seems to be mostly apparent up to the 160-loci datasets ( Figure 5), but this threshold has to be interpreted with caution, given that randomly selected datasets tend to share increasingly more similar with the ranked dataset as more loci are added simply as a result of the finite total dataset.

DISCUSSION
We carried out an extensive exploration of phylogenomic UCE data for species delimitation. randomly selected loci from the full set of informative loci. In the first dataset, only one set of priors led to convergence (runs with small ancestral population size, i.e. θ = 2, 2000), whereas on the second dataset did not reach convergence, even after 300 k iterations. Oswald et al. (2016) used Bayes factor species delimitation using SNPs drawn from UCEs with the BFD* method (Grummer et al. 2014;Leaché et al. 2014). Given that only one SNP was used per locus to minimize genetic data linkage, the difference in their chosen method means that their results are not directly comparable to those in this study. Nevertheless, as UCEs and other types of genomic tool become increasingly more available (e.g. Starrett et al. in press; Branstetter et al. in press, Faircloth in press), we anticipate that phylogenomic data might become a valuable tool for species discovery, delimitation, and diagnosis. In particular, our results are consistent with Hime et al. (2016) in that, although species often can be delimited with relatively few loci, sampling larger numbers of loci might be necessary to ensure that recent divergences are correctly identified. In that study, Hime et al. explored variation from 2 to 10 loci using SpedeSTEM (Ence and Carstens 2011) and from 10 to 89 loci in analyses using BPP, in 10-locus increments.
In those analyses, for each subsampling increment, they generated 10 replicate jackknifed datasets. They also ranked loci based on their informativeness, with increasing numbers of loci (as in the present study). Both the SpedeSTEM and BPP analyses showed more fine-scale delimitation with increasing number of loci, as expected based on simulation studies (e.g. Camargo et al. 2012;Hird et al. 2010). Indeed, our results strongly suggest that increasing the number of loci had a greater effect on delimitation than ranking them based on informativeness Pie et al. page: 14 and focusing on the most informative loci, which suggests that uninformative loci might be more problematic for phylogenetic inference than for species delimitation analyses (at least in the case of BPP).
The use of genetic data for species delimitation has been recently criticized by Sukumaran and Knowles (2017), who argued that species delimitation based on the MSC, particularly as implemented in BPP, diagnoses genetic structure, not species, and that it cannot statistically distinguish structure associated with population isolation vs. species boundaries. The main concern raised by Sukumaran and Knowles (2017) is the possibility of taxonomic inflation if species are described based on genetic data alone (see also Olave et al. 2014). We believe that their concern is highly relevant in theory, but a relatively moot point in practice for two main reasons. First, although many studies in recent years have used genetic species delimitation methods, very few actual species descriptions have been based on genetic data alone. Second, although conservation and management efforts ordinarily use species numbers as their main currency, such that taxonomic inflation could potentially lead to flawed conclusions, it is important to note that the opposite problem (the failure to recognize cryptic species as distinct Pie et al. page: 15 entities) is at least equally as important (Bickford et al. 2007). Molecular species delimitation is often used precisely when diagnosis using phenotypic traits alone is difficult. It is therefore ironic that the recommendation is to go back to the phenotype, which is often (tacitly) assumed to be more reliable. Yet one might easily forget that phenotypic traits are often plagued with the same kinds of uncertainty as those indicated for molecular data, such as the difficulty in finding fixed diagnostic traits among potentially distinct species (Wiens and Servedio 2000), and the challenge of discriminating between differentiated populations from "good" species. It is possible that some of the criticism of molecular species delimitation methods has been due to the naive optimism of some early efforts that it would finally provide an objective means of defining species (as in the early days of molecular taxonomy), but the practice of taxonomic research has already demonstrated that species delimitation will never be an automated task, with or without phenotypic traits. In the end, it is important to keep in mind that cryptic species do exist, and that their biological reality should not be negated simply because fixed phenotypic diagnostic traits have not been identified.
Our results strongly corroborate the evolutionary distinctiveness of anuran species that had been previously described using phenotypic data alone (e.g. Bornschein et al. 2015Bornschein et al. , 2016bPie and Ribeiro 2015;Ribeiro et al. 2015). This is significant, given that many of those speciesparticularly in the case of Brachycephalus -have been diagnosed to a large extent based only on coloration, a trait that is often considered unreliable for anuran taxonomy. In addition, our analyses found evidence for 3 undescribed species of Brachycephalus and that three currently verrucosus X B. olivaceus -see Figure 3) did not reach 1.0 posterior probability in all analyses, despite having clear morphological diagnoses (including color variation), suggesting that BPP might in fact be conservative in its species delimitation compared to phenotypic evidence. The confirmation of the genomic distinctiveness of the montane anuran species in our study that had been previously described using phenotypic data alone, and the possibility of additional new species that require further study, has important implications for their conservation. Cloud forests, which are the habitat of most of the species in the pernix group, are among the most threatened ecosystems globally (Doumenge et al. 1995;Aldrich et al. 1997;Toledo-Aceves et al. 2011). This is of particular concern, not only due to the key role played by these forests in hydrological cycle maintenance, but also because they are reservoirs of endemic biodiversity (Toledo-Aceves et al. 2011). Many of the species in this study are categorized as either threatened or as data deficient, yet any level of formal governmental protection necessarily involves the availability of a species name. Given their microendemic distribution and highly threatened habitats, one could argue that commission errors (considering them as different species) would be preferable in relation to omission errors (lumping them as single species). In practice, urgent management efforts should be enforced to ensure their long-term survival.