Phylogenetic analysis and species delimitation.
The evolutionary model and partition scheme (in the case of coding
genes) that best fit each dataset were estimated in ModelFinder
(Kalyaanamoorthy et al. 2017)
according to the Bayesian Informative Criterion (Table 1). Maximum
Likelihood (ML) and Bayesian Inference (BI) reconstruction methods were
used to estimate phylogenetic relationships for each gene and the
concatenated datasets. ML analyses were run in IQ-Tree
(Nguyen et al. 2015) using 10000
Ultrafast Bootstrap replicates (UFBoot;
Minh et al. 2013) and the GENESITE
resampling strategy to estimate branch support. The BI analyses were run
in MrBayes v3.2.6 (Ronquist and
Huelsenbeck 2003) using eight chains in two independent runs (10 million
generations each, sampling once every 1000 generations). We used Tracer
v1.7.1 (Rambaut et al. 2018) to
verify the convergence and seasonality of each run. The posterior
probability (pP) was obtained for individual nodes by constructing a
majority-rule consensus after discarding the trees prior to the
stationarity phase as burn-in. These analyses were implemented within
the CIPRES Science Gateway portal
(Miller et al. 2012).
We performed three molecular species delimitation methods: multi-rate
Poisson Tree Processes (mPTP;
Kapli et al. 2017), Bayesian
General Mixed Yule-Coalescent Model (bGMYC;
Reid and Carstens 2012), and
Species Tree and Classification Estimation, Yarely (STACEY,
Jones 2017). These methods were
selected because they do not require a priori assignment of individuals
to groups. This allowed us to define cryptic lineages within R.
mexicanus sensu lato without assumptions regarding their phylogenetic
relationships.
The cytb gene was used for the two single-locus methods mPTP and bGMYC.
The non-coalescent mPTP method was implemented in the Exelisis Lab
platform (http://www.exelixis-lab.org )
using the BI tree and default parameters. The bGMYC coalescent method
was implemented in the bGMYC package
(Reid and Carstens 2012) of the R
library (R Core Team 2018) using
the time-calibrated tree obtained in BEAST2 and the following input
parameters: mcmc = 100000, burnin = 90000, thinning = 100, t1 = 11, t2 =
16 (based on the upper range of suggested species with mPTP, considering
the outgroup), py1 = 0.5, py2 = 1.5, pc1 = 0.1, pc2 = 0.5, start = c
(1.0, 0.1, 11), scale = c(20, 10, 5.00). The multi-locus STACEY method
was run using each of the combined datasets (cytb + Fgb-I7 and cytb +
IRBP). This coalescent method is implemented as a package within BEAST2
and uses a modified birth-death-collapsed model for the species tree.
The species.tree file generated by STACEY was used as input in the
SpeciesDelimitationAnalyzer program (speciesDA.jar,www.indriid.com ; burnin = 1000
collapse height = 0.0001, and similarity cutoff = 1.0) to summarize the
posterior tree distribution and calculate the frequency with which each
pair of taxa were assigned to the same clade.
The criterion used to define species limits was congruence among the
greatest number of delimitation methods, which supports the correct
recognition of putative species
(Carsten et al. 2013). Genetic
distances between lineages proposed as species were estimated using the
cytb gene in MEGAX (Kumar et al.
2018) under the Kimura two-parameter evolutionary model (K2P). The 5%
cytb distance value associated with mammalian sister species recognition
(Bradley and Baker 2001) was set as the genetic distance threshold.