Introduction

The family Stenopelmatidae Burmeister, 1838 is a group of cricket-like orthopterans whose species occur in North and Central America, South and Central Africa, Southeast Asia, and the Indian subcontinent (^{Cigliano et al., 2022}). This family was recently proposed by ^{Gorochov (2021)} to contain 3 subfamilies: Stenopelmatinae Burmeister, 1838, Schizodactylinae Blanchard, 1845, and Gryllacridinae Blanchard, 1845. However, only the former subfamily is currently recognized within Stenopelmatidae in the Orthoptera Species File (^{Cigliano et al., 2022}). Following the latter classification, the Stenopelmatinae comprises 5 tribes: Maxentiini Gorochov, 2021, Oryctopini Kevan, 1986, Oryctopterini Gorochov, 1988, Siini Gorochov, 1988, and Stenopelmatini Burmeister, 1838 (^{Cigliano et al., 2022}; ^{Gorochov, 2021}). Of these tribes, Stenopelmatini is the only one present in the American continent, with species distributed from British Columbia in Canada along the west coast of the USA to the Mesoamerican highlands from northern Mexico to Panama, and with one species from northern South America in Ecuador (^{Cigliano et al., 2022}; ^{Weissman et al., 2021}).

Members of Stenopelmatini, commonly known as Jerusalem crickets in the U.S. and as ‘cara de niño’ or ‘niño de la tierra’ in Mexico, are amongst the most recognizable insects in North and Mesoamerica, (^{Weissman, 2005}). They are mainly nocturnal, hiding under rocks or logs during the day (^{Weissman, 2001}, ²⁰⁰⁵). Species of this tribe mainly feed on roots and are considered important organic matter decomposers (^{Weissman, 2001}), although they also feed on small invertebrates (^{Weissman, 2001}).

In a molecular phylogenetic study for the Stenopel- matoidea based on 2 nuclear and 1 mitochondrial (mt) gene markers (^{Vandergast et al., 2017}), the Stenopelmatini was recovered as monophyletic, although, only members of another tribe, Siini, were included. In a more recent integrative phylogenetic study that focused on the Stenopelmatini, ^{Weissman et al. (2021)} recovered 3 strongly supported main clades for the subfamily based on both anchor hybrid enrichment and mt data. One of these clades had the micropterous S. piceiventris Walker, 1869 from Oaxaca in southeast Mexico, whereas the remaining 2 clades included species from western North America, Canada to northern Mexico, and Mesoamerica, respectively. Based on these relationships and examination of morphological, behavioral, and chromosomal features, the tribe was restricted to contain only 2 genera: Ammopelmatus Tinkham, 1965 for the North American clade and Stenopelmatus Burmeister, 1838 for the 2 clades with Mesoamerican species.

The genus Stenopelmatus currently contains 19 described macropterous, micropterous, and apterous species that mainly inhabit highlands with pine-oak and cloud forests (^{Cigliano et al., 2022}; ^{Weissman et al., 2021}). A recent phylogenomic study showed that the species richness in this genus is considerably higher, with several undescribed species with a marked geographic phylogenetic structure (^{Gutiérrez-Rodríguez et al., 2022}). Based on these results, 4 species-groups were proposed: the faulkneri, talpa, Central American and piceiventris species-groups, whose geographic distribution and molecular clock estimates were concordant with some geological processes that took place in the Mexican Transition zone (^{Gutiérrez-Rodríguez et al., 2022}).

The S. talpa species-group was proposed by ^{Gutiérrez-Rodríguez et al. (2022)} to contain 2 described (S. talpa Burmeister, 1838 and S. typhlops Rehn, 1903) and 5 undescribed species that occur across the Trans-Mexican Volcanic Belt (TMVB) and adjacent areas of the Sierra Madre Oriental (SMO) morphotectonic provinces in central Mexico. However, despite the effort of the above cited authors to assess the species boundaries in this group, additional information from different sources is still required to confirm its species limits, phylogenetic relationships, and biogeographic history.

Mt sequence data have been widely employed to elucidate evolutionary relationships, establish species boundaries, and carry out population genetics studies (^{Avise et al., 1987}; ^{Ballard & Pichaud, 2014}). The metazoan mt genome typically consists of 15-18 kilobases, comprising 13 protein-coding genes, 22 transfer RNAs (tRNAs), and 2 ribosomal RNAs (rRNAs) (^{Boore, 1999}), with this composition being generally conserved across bilaterian metazoans (^{Cameron, 2014}; Shen et al., 2016). The use of this locus has given access to a set of orthologous genes with little or no recombination (^{Ballard & Rand, 2005}) and has proved to provide valuable phylogenetic resolution at both shallow and deep evolutionary scales (^{Ballard & Rand, 2005}; ^{Cameron, 2014}; ^{DeSalle et al., 2017}; ^{Lumley & Sperling, 2010}; ^{Song et al., 2015}).

With the advent of next-generation sequencing (NGS) techniques, phylogenetic studies of a number of biological groups, including insects, have enormously benefited due to the facilitated generation and use of genomic information, including the assembly of complete mt genomes (^{Maddock et al., 2016}). In this respect, it has been shown that phylogenetic analyses of complete mt genome data increase resolution for resolving species diversification events (^{Bae et al., 2004}; ^{Song et al., 2015}). Here we generated and characterized 16 complete and partial mt genomes of Stenopelmatus species, 14 of which belong to species of the S. talpa species-group that were delimited in ^{Gutiérrez-Rodríguez et al.’s (2022)} study, and whose validity still needs to be further corroborated. We employed this mt genome data to assess the species limits among the examined taxa, estimate their phylogenetic relationships, and investigate their biogeographical history in central Mexico. Our results support the existence of at least 9 species within the S. talpa species-group, 2 of which have already been described, whose geographic structure is congruent with the geological history of the TMVB and adjacent morphotectonic provinces.

Materials and methods

A total of 16 Stenopelmatus specimens were included in this study, 14 of which belong to the S. talpa species- group. These specimens were collected in localities situated in 7 Mexican states across the TMVB and SMO provinces (Hidalgo, Mexico City, Michoacán, Puebla, Querétaro, State of Mexico, and Veracruz; Fig. 1). Our sampling included specimens assigned to the only 2 described species of the S. talpa species-group, S. talpa and S. typhlops, which were collected from their type locality. We also included 2 undescribed species from 2 other species-groups: one belonging to the faulkneri species- group, S. sp. aff. mineraldelmonte (El Chico, Hidalgo, CNIN3922), and a second to the Central America species- group, rooting the trees with the latter species. A list with the examined species, their locality, DNA voucher, and GenBank accession numbers is provided in Table 1. All samples were collected between the years 2012 to 2018 and preserved in 96% ethanol at -20°C until they were processed for DNA sequencing.

Figure 1 Map with the localities of the specimens examined in this study. The map shows the Trans-Mexican Volcanic Belt morphotectonic province. 

Genomic DNA was extracted from smooth muscle of hind femur with the EZ-10 Spin Kit minipreps DNA Genomic Column Kit (BIOBasic, Toronto, Ontario, Canada) following the manufacturer’s protocol. DNA fragmentation was performed with the Bioruptor Pico equipment (Diagenode Inc.) employing 3 cycles of sonication (15-90 sec. on-off pulse). We used 100 ng of input DNA for library preparation with the Kapa Hyper Prep kit (Kapa Biosystems Inc., Wilmington, MA, U.S.A), and the TruSeq-style dual-indexing adapters (^{Glenn et al., 2016}). Genomic DNA quantitations, pre- and post-PCR libraries, were performed using a Qubit 2.0 fluorometer (Invitrogen, Life Technologies, CA, USA). Whole-genome shotgun sequencing was conducted with an Illumina HiSeq X instrument at the Department of Environmental Health Science, University of Georgia, Athens, GA, USA.

Raw reads were trimmed and filtered for all samples using the program Geneious version 10.2.6 (^{Kearse et al., 2012}). We carried out mt genomes de novo assemblies in the program GetOrganelle (^{Jin et al., 2018}) using the mtDNA database. For 3 samples (CNIN 3669, 3741, 3922) we could not recover complete mt genomes in the de novo assembly. In these cases, we performed by-reference assemblies with Geneious to obtain longer contigs, using as reference either mt sequence genomes of geographically close samples or their cytochrome oxidase I (cox1) sequence previously obtained by Sanger sequencing. Mt genome annotations were performed with the MITOS 2 web server (^{Bernt et al., 2013}) and verified using the protein-coding genes signal from the “protein plots” generated by MITOS. We used Geneious to confirm the accuracy of our assemblies and annotations. We also registered the order of the mt genome genes to identify possible gene arrangements, and then compared them to the Pancrustacea ground pattern, which is the proposed Crustacea/Hexapoda common ancestor (^{Boore et al., 1995}, ¹⁹⁹⁸).

We extracted the sequences of the 13 protein-coding genes from all samples and independently aligned them using the program MAFFT version 7 (^{Katoh & Standley, 2013}). We then verified each protein-coding gene alignment with respect to the reading frame (invertebrate mt genetic code) in Geneious. The 2 ribosomal genes (rRNA) and the 22 transfer genes (tRNA) from each sample were also extracted and manually aligned. The following 3 matrices were obtained: 1) a cox1 matrix comprising the DNA barcoding locus (658 bp), 2) a protein-coding gene matrix (10,953 bp), and 3) a complete mt genome matrix (14,599 bp). The 3 matrices are available in figshare (https://doi.org/10.6084/m9.figshare.21506583).

Uncorrected genetic distances were calculated both for the complete mitogenome and cox1 matrices with the program PAUP version 4.0a (^{Swofford, 2002}). Mt species delimitation among the samples was assessed for the cox1 data set following a 2% pairwise distance criterion. The use of this marker has been proven to be an efficient tool for species identification of most animal taxa (^{Hajibabaei et al., 2007}; ^{Hebert et al., 2003}). We included in this analysis 4 previously published cox1 sequences belonging to specimens of the talpa species-group: 3 from the state of Morelos and one from Mexico City (^{Gutiérrez-Rodríguez et al., 2022}). The uncorrected cox1 distances were visualized building a neighbor-joining (NJ; ^{Saitou & Nei, 1987}) tree with PAUP.

The best-fit partitioning scheme for the protein-coding gene and complete mt genome matrices was obtained with ModelFinder (^{Kalyaanamoorthy et al., 2017}) in the program IQTree version 2 (^{Minh et al., 2020}). For both matrices we conducted maximum likelihood (ML) phylogenetic analyses using the same version of the program IQTree with 1,000 ultra-fast bootstrap (BTP) replicates.

We carried out the divergence times estimates for the complete mitogenome matrix with the program BEAST2 version 2.4.7 (^{Bouckaert et al., 2014}) using the Birth- Death model (^{Nee, 1994}) and an uncorrelated lognormal relaxed clock (^{Drummond et al., 2006}). A single partition was considered, using the GTR + I + G model of evolution. To calibrate the tree, we used the mutation rate reported by ^{Papadopoulou et al. (2010)} for the mt genome (2.4% per My; 0.012 ± 0.0012). The analysis ran for 100 million generations, sampling trees every 10,000 generations with a burn-in of 25%. The Maximum Clade Credibility Tree (MCCT) was obtained with TreeAnnotator version 1.7.4 (^{Drummond et al., 2012}).

Results

We obtained 7 complete and 9 partial mt genomes with lengths between 13,056 and 15,785 base pairs, and a GC proportion between 25.0 and 26.5% (Table 2). Of the 16 mitogenome sequences, 12 were recovered as a single contig, whereas the remaining 4 were recovered in 2 or more fragments. We did not find arrangement variation in any of the assembled mitogenomes, with all of them being similar to the Pancrustacea ground pattern (^{Boore et al., 1998}) (Fig. 2, Appendix 1).

Figure 2 Stenopelmatus mitochondrial genome composition. Protein coding genes and tRNAs are in pink color and rRNAs in orange. Linear representation of the mitochondrial genome composition is provided in Appendix 1. 

Uncorrected genetic distances for the cox1 and complete mitogenome matrices are shown in Table 3. The pairwise distances observed for the complete mt genome data set generally were slightly lower compared to the ones obtained with cox1. The NJ tree built with the cox1 data set (Fig. 3a) grouped the ingroup taxa in 9 separate clusters that were mostly congruent with the sequence divergence gap of 2%. One of the cases of incongruence was found within the putative S. sp. 1 from the Valley of Mexico, where 6 of the 8 pairwise comparisons that involved the specimen from Nepantla, Estado de Mexico (CNIN3669) were marginally higher than 2% (2.1-2.59%). Furthermore, the only sample from Apatzingán, in the state of Michoacán (S. sp. 7; CNIN3741), had a distance slightly lower than 2% with a sample of S. sp. 1 from Mexico City (CNIN3669; 1.87%) and the sample belonging to S. sp. 3 from Morelia, Michoacán (CNIN3737; 1.9%).

Figure 3 a) Neighbor-joining tree reconstructed using the cox1 barcoding locus (658 bp) showing the 9 species delimited for the S. talpa species-group based on the 2% sequence divergence criterion. b) Phylogram for the S. talpa species-group reconstructed with the maximum likelihood method using the complete mitogenome sequence. Numbers near nodes represent ultra-fast bootstrap values. 

The ML phylograms derived from the complete mt genome and protein-coding gene matrices showed identical relationships, only varying in the BTP support values of some of their clades (Fig. 3b and Appendix 2, respectively). The putative species S. sp. 7 from Apatzingán, Michoacán, in the northernmost part of the Sierra Madre del Sur (SMS) province was sister to the remaining ingroup taxa (BTP = 100), which appeared divided into 2 main clades. One of these clades (BTP = 79) had S. talpa, which is distributed in the limits between the SMO and the TMVB provinces, as sister to the taxa from the eastern part of the TMVB. The second main clade (BTP = 100) contained the putative species from the central portion of the TMVB in the Valley of Mexico and central Michoacán.

The reconstructed chronogram (Fig. 4) shows that the origin of the S. talpa species-group took place during the late Miocene to late Pliocene, 5.81-2.38 Mya, whereas the divergence between the putative undescribed species from the northernmost portion of the SMS in Apatzingán, Michoacán, and those from the TMVB and SMO provinces probably occurred during the late Pliocene 3.55-2.5 Mya. The split of the members from the central and eastern portions of the TMVB on the other hand probably occurred during the late Pliocene to early Pleistocene, between 2.94 and 1.32 Mya.

Figure 4 Bayesian ultrametric tree obtained with BEAST2 using the complete mitochondrial genome data set for the S. talpa species- group. Gray bars represent geological time periods. 

Discussion

Mitogenome features in the talpa species-group

Although the Orthoptera is one of the orders whose species are known to have among the largest nuclear insect genomes and thus are more prone to the insertion of nuclear mitochondrial paralogs (^{Song et al., 2014}), its mt gene arrangement has remained generally stable through time (^{Cameron, 2014}; ^{Song et al., 2015}). Currently, the only mt gene rearrangement that has been reported as a probable synapomorphy at the supraspecific level within the order is the trnD-trnK tRNA gene translocation within Caelifera (^{Song et al., 2015}).

Here, we have generated and assembled the complete and nearly complete mt genome of various Stenopelmatus specimens, a poorly-known group of ensiferan orthopterans mainly distributed across the Mesoamerican highlands. Previous to this study, only one mitogenome had been published for the genus and the whole family Stenopelmatidae, i.e., Ammopelmatus fuscus Haldeman, 1852 (^{Song et al., 2015}). Thus, our results provide valuable information for a better understanding of the molecular evolution of this orthopteran family and Ensifera in general.

Our results indicate that all the assembled mt genomes follow the ancestral insect mt genome structure (^{Boore et al., 1998}), thus supporting the high conservativeness of the mt genome structural organization within Orthoptera. Also, other mt genome features among the examined species, such as the conservation of the trnK-trnD pattern of the tRNAs between the protein-coding genes cox2 and atp8, appear to be similar to those described for other ensiferans (^{Fenn et al., 2008}; ^{Song et al., 2015}) as well as a high AT content reported for orthopterans in general (^{Fen et al., 2008}; ^{Sheffield et al., 2010}; ^{Song et al., 2015}). Further assembling of additional mt genomes of other stenopelmatid taxa will reveal whether the ancestral mt gene organization and other relevant mt features persist in all members of the family.

Barcoding-based species delimitation

A recent molecular phylogenetic study based on mt and 3RAD nuclear data showed the existence of an extensive, highly neglected species diversity in Stenopelmatus (^{Gutiérrez-Rodríguez et al., 2022}). Of the 34 species that were delimited in the latter study, 8 belonged to a clade named the S. talpa species-group, whose species are characterized by having a large size and a generally reddish-orange and black color pattern (although it appears to be exclusively black in S. typhlops; ^{Gutiérrez-Rodríguez et al., 2022}).

Our inclusion of additional specimens from the same and newly sampled localities with respect to the above study resulted in 9 barcoding-based species, although the delimitation of samples from Mexico City and surrounding localities decreased to only one putative species in comparison with ^{Gutierrez-Rodríguez et al.’s (2022)} study, which discriminated 3 potential evolutionary lineages. The species delineation obtained here, however, must be taken with caution, since some of the cox1 distances that involved specimens from the populations that had been regarded as a different lineage (Nepantla, Estado de Mexico) in some cases exceeded the 2% barcoding gap. Moreover, the specimen from Orizaba, Veracruz was also shown to represent a separate species. This strongly suggests the validity of S. mexicanus Saussure, 1859, a name that was recently regarded as nomen dubium and whose lectotype female was collected in Córdoba, which is less than 25 km in a straight line from Orizaba. A further molecular and morphological examination for all the aforementioned populations will help to consistently establish their species boundaries and taxon assignation.

The relationships obtained here for the members of the S. talpa species-group based on mt genome DNA sequence data were mainly similar to those recovered in a recent phylogenetic study of Stenopelmatus based on a fragment of the cox1 mtDNA gene and nuclear 3RAD data (^{Gutierrez-Rodríguez et al., 2022}). Our divergence times estimates, however, were considerably younger than those obtained by ^{Gutierrez-Rodríguez et al. (2022)}, probably due to the different source of calibration that they employed. We calibrated our molecular clock analysis taking advantage of a widely employed mutation rate that has been established for the mt genome (2.4% per My; ^{Papadopoulou et al., 2010}). In contrast, the above study employed a secondary calibration for the cox1 dataset for the time of divergence between Stenopelmatidae and Anostostomatidae Saussure, 1859 (117.4 My; ^{Song et al., 2015}) instead of using a mutation rate, since the use of single mitochondrial markers for deep phylogenies is prone to have considerably high levels of nucleotide saturation (^{Molak & Ho, 2015}). Below we thus make some inferences about the biogeography of the S. talpa species-group based on the relationships and the times of divergence obtained here, which support its recent origin and species diversification (^{Weissman et al., 2021}).

The TMVB is a 1,000 km long Neogene continental arc that is situated across central Mexico, which became an independent province from other partially overlapping arcs during the early to middle Miocene ~10 Mya, originating from the eastward subduction of the Farallon plate beneath western Mexico (^{Ferrari et al., 2012}). The TMVB is characterized by having an east-west (E-W) orientation, being divided into 4 distinct episodes: 1) early- to mid- Miocene, 2) late Miocene, 3) late Miocene-early Pliocene, and 4) late Pliocene and Pleistocene (^{Ferrari et al., 2012}; ^{Gómez-Tuena et al., 2007}). In addition, this province is divided into 4 sectors according to their age of inception: Easternmost (~16 Mya), Eastern (~19 Ma), Central (~11 Ma) and Western (~11 Ma) (^{Ferrari et al., 2012}).

According to our results, the S. talpa species-group apparently followed an east-central diversification pattern along the TMVB province. The earliest divergence within the S. talpa species-group was estimated to occur ca. 3.55-1.58 Mya and involved S. sp. 7, a species distributed between the SMS and TMVB border, followed by the clade with the species from the TMVB and adjacent areas of the SMO provinces, which started to diverge ca. 2.94-1.32 Mya. Moreover, most of the earlier divergences within the group occurred along the northern borders of the Eastern sector, while the most recent ones can be found across the Eastern and Central sectors. This is congruent with the aforementioned topographic evolution of the TMVB, which began in the early to mid-Miocene and finished in the late Pliocene to Pleistocene.

The biogeographic inferences drawn here for this apparently recently diverged group of orthopterans endemic to central Mexico will serve as a basis for future evolutionary studies of other groups of organisms with a similar geographic distribution. This work also highlights the utility that the mt genome sequence data have for phylogenetic reconstruction studies of animal taxa, and advocates for more studies that make use of this valuable, highly informative locus.