Comparative analysis of mitogenomes among three species of grasshoppers (Orthoptera: Acridoidea: Gomphocerinae) and their phylogenetic implications

Sample collection and DNA extraction

The samples used in this study are shown in Table 1. Specimens were preserved in 95% alcohol and then transferred to 4 °C for cryopreservation. We performed an accurate identification under stereomicroscope based on its morphological characteristics by Fauna Sinica (Yin & Xia, 2003). Using the Animal Tissues/Cells Genomic DNA Extraction Kit (Solarbio, Beijing, China) and according to the manufacturer’s instructions, total genomic DNA was extracted from the hind femora muscle. The quality of the total DNA was checked with a 1% agarose gel and the concentration was measured with a Nanodrop 2000 spectrophotometer. DNAs were stored at −20 °C for long term storage and further molecular analyses.

DNA sequencing and assembly

Complete mitochondrial genomes were sequenced using the Illumina Novaseq 6000 platform with 151 bp paired end reads at Personalbio, Shanghai, China. Using the invertebrate genetic code in Genious 8.1.3 (Kearse et al., 2012), by using the closely related known grasshoppers as reference sequences, complete mitochondrial genomes of Dasyhippus peipingensis, Myrmeleotettix palpalis, and Aeropedellus prominemarginis were assembled and aligned by using ClustalW (Larkin et al., 2007). Complete mitochondrial genome sequences were manually proofread in Genious 8.1.3 to check the accuracy of the assembly.

Mitogenome assembly and annotation

The mitochondrial genome sequences were annotated using Genious 8.1.3 and online website GeSeq (https://chlorobox.mpimp-golm.mpg.de/geseq.html). The mitochondrial genome map was visualized by Proksee (https://proksee.ca/). The tRNA Scan-SE v. 1.21(Lowe & Eddy, 1997) was used to identify tRNAs and confirm their secondary structure, and the secondary structure was visualized using Adobe lllustrator 2020. The nucleotide composition, base composition skew, codon usage and relative synonymous codon usage (RSCU) of protein-coding genes were analyzed by PhyloSuite v1.2.2 (Zhang et al., 2019). The skew values were calculated using the formulae: GC skew = (G − C)/(G + C) and AT skew = (A − T)/(A + T) (Perna & Kocher, 1995). The “ggplot2” package in the RStudio was used to visualize relative synonymous codon usage (RSCU). Nucleotide diversity (Pi) and sliding window analysis (sliding window: 100 bp, step size: 25 bp) of 13 PCGs among 10 Gomphocerinae species was performed using DnaSP 5.0 (Librado & Rozas, 2009). Non-synonymous substitutions (Ka) and synonymous substitutions (Ks) were calculated for total PCGs by KaKs Calculator Toolbox 2.0 (Wang et al., 2009), and the evolutionary rate of PCGs was assessed by Ka/Ks value.

Phylogenetic analyses and divergence time estimate

The complete mitochondrial genomes of 34 species from 13 subfamilies within family Acrididae were selected as ingroup, and seven species within Pamphagidae, Pyrgomorphoidea, and Tetrigoidea were chosen as outgroups (Table 2). All mitochondrial genomes (except the three sequenced in this study) were obtained from GenBank (accession numbers given in Table 2). The heterogeneity of nucleotide divergence of two matrixes were analyzed by AliGROOVE 1.5 (Kück et al., 2014). The nucleotide sequences of the 13 protein-coding genes were analyzed. Substitution saturation of PCGs based on Xia’s test implemented in DAMBE (Xia, 2017). The 13 PCGs were extracted by PhyloSuite v1.2.2 and aligned using MAFFT v7.313 (Katoh & Standley, 2013). Intergenic gaps were removed by Gblocks v 0.91b (Castresana, 2000), and the 13 protein-coding genes were concatenated in PhyloSuite v1.2.2.

The best-fit models of each gene were selected by ModelFinder (Lanfear et al., 2017), maximum likelihood (ML) and Bayesian inference (BI) were selected based on AICc, the results are presented in Tables S1 and S2. ML and BI analyses were performed in IQ-TREE (Nguyen et al., 2015) and MrBays v3.2.7 (Ronquist & Huelsenbeck, 2003), respectively. The ML phylogenetic analyses were using “Ultrafast” algorithm (bootstrap number = 5,000). In BI phylogenetic analyses, Markov chain Monte Carlo (MCMC) run for 10,000,000 generations, sampling trees every 10,000 generations with the 25% burn in. The data were imported into Tracer v.1.7.2 (Rambaut et al., 2018), ESS >200 proving that the data were converged. The consensus trees were displayed and managed visually by iTOL (https://itol.embl.de) (accessed February 2023).

Divergence time in Aricidiae was estimated using the 13 PCGs with relaxed molecular clock model in BEAST 1.10.4 (Drummond et al., 2012). Coalescent: Constant Size model was used for the prior tree, ModelFinder was used to find the best model GTR+G+F4. Divergence time tree nodes from Chang et al. (2020) were used for calibration (separation time 115 Mya for Tetrigoidea; 71 Mya for Chrotogonidae; 56 Mya for Pamphagidae; 35 Mya for Catantopinae; 33 Mya for Oedipodinae). The Markov chain was run 100,000,000 generations, sampling every 10,000 generations, 25% was burn in. The stability of the results was verified by Tracer v1.7.2 with most parameters having more than 200 effective sample size (ESS) values. ChiPlot (https://www.chiplot.online/tvbot.html) was used to visualize maximum clade credibility tree with 95% highest probability density (95% HPD).

Genome content and organization

We sequenced and annotated the whole mitochondrial genomes which performed visual editing. The complete mitogenome sequence of Dasyhippus peipingensis, Myrmeleotettix palpalis, and Aeropedellus prominemarginis were 15,628, 15,621 and 15,629 bp, respectively. These mitogenomes showed typical insect mitogenome structure, which composed of circular double-stranded DNA molecules (Fig. 1). Each mitogenome includes 13 protein-coding genes (PCGs), 22 tRNA genes, two rRNA genes and an A+T rich region (control region). There were 23 genes (including nine PCGs and 14 tRNAs) are encoded on majority-strand (J-strand) and 14 genes (including four PCGs, eight tRNAs and two rRNAs) are transcribed from the minority-stand (N-strand) (Table 3).

The nucleotide compositions of the three mitochondrial genomes revealed a distinct A/T bias: 75.7% (Dasyhippus peipingensis), 75.2% (Myrmeleotettix palpalis), and 75.0% (Aeropedellus prominemarginis). All mitochondrial genomes were positive for A+T skew and negative for GC skew (Table 4). The complete mitochondrial genomes and PCGs of three grasshopper species had A+T contents higher than 64% at different compositional sites and locations. The control regions (A+T-rich region) of the mitochondrial genomes were all located between tRNA-Ile and rrnS, with sizes of 729 bp (Dasyhippus peipingensis), 728 bp (Myrmeleotettix palpalis), and 735 bp (Aeropedellus prominemarginis), and the A+T content was >83%, which were also referred to as AT-rich regions. The structures and nucleotide compositions of the three species are generally consistent with the mitochondrial genome structure of the Acrididae (Zhang et al., 2023; Zheng et al., 2021), indicating that these mitochondrial genome’s structure is highly conserved (Wei et al., 2010).

The total lengths of the 13 PCGs of Dasyhippus peipingensis, Myrmeleotettix palpalis and Aeropedellus prominemarginis are 11,190, 11,193 and 11,190 bp, respectively (Table 4), accounting for 71.6%, 71.7%, and 71.6% of the whole mitochondrial genome, respectively. The size of the PCGs ranges from 162 bp (ATP8) to 1,719 bp (ND5). Among the 13 PCGs, 9 PCGs (ATP6, ATP8, COX1, COX2, COX3, CYTB, ND2, ND3, and ND6) are encoded on the J-strand, while 4 PCGs (ND1, ND4, ND4L, and ND5) are encoded on the N-strand. The third codon position has the highest A+T content, while the second codon position has the lowest A+T content. All the initiation codons in the mitogenomes of the three species were ATN, with ATG being the most frequently used, termination codons were TAN, and TAA was the most frequently used termination codon.

To indicate the frequency of codon usage, the relative synonymous codon usage (RSCU) values of the three mitochondrial genomes were visualized (Fig. 2; Table S3). Comparative analysis showed that the synonymous codon preferences were highly conserved among three mitochondrial genomes. The most frequently used codons are TTT, TTA, ATT, and ATA, therefore, Phe, Leu (UUR), IIe, and Met are most frequently used amino acids, accounting for 7.86%, 9.67%, 9.31%, and 5.81% of total, respectively. In addition, RSCU analysis also showed a bias towards using more A/T at the third codon position rather than G/C. Similarly, the frequency of codon usage indicates the preference of nucleotide A/T in three species.

As other insects, the mitochondrial genome of three species contain 22 tRNA genes with lengths ranging from 62–71 bp (Table 3), the total length of tRNAs is ranging from 1,475 bp (Dasyhippus peipingensis) to 1,476 bp (Myrmeleotettix palpalis and Aeropedellus prominemarginis). The A+T content of tRNAs is 73.2%, 72.6%, and 72.7% for three mitochondrial genomes, with positive AT skew and GC skew. Most tRNAs could be folded into the typical clover leaf secondary structure, except that tRNA-Ser (AGN) lacked a dihydrouridine (DHU) arm and formed a simple loop (Fig. 3). The secondary structure of tRNAs is usually conserved in the amino acid acceptor arm and anticodon loop, while DHU and TψC are more variable. In addition to the classic base pairs A-U and C-G, there are also noncanonical base pairings (G-U and A-C) and mismatched base pairs (A-A and A-G) distributed throughout the tRNA arms, with G-U noncanonical base pairs being the most abundant.

The two ribosomal RNA genes are encoded on the N-strand among three grasshoppers (Table 3), rrnL is located between tRNA-Leu (CUN) and tRNA-Val, while rrnS is flanked by tRNA-Val and A+T rich regions. The rrnL of Dasyhippus peipingensis, Myrmeleotettix palpalis and Aeropedellus prominemarginis are 1,319, 1,316 and 1,312 bp in length, contains the A+T content ranging from 77% to 78%. The rrnS is 843bp in Dasyhippus peipingensis and Myrmeleotettix palpalis, 844 bp in Aeropedellus prominemarginis, with AT content ranging from 73.6% to 75.5%. Therefore, there were no significant differences in rRNAs among three species. Both rrnL and rrnS exhibit negative AT-skew and positive GC-skew in three mitogenomes.

Nucleotide diversity analysis can identify regions with large nucleotide divergence, which is useful for designing species-specific markers in groups within taxa where morphological identification is difficult and taxonomic boundaries are blurred (Ma et al., 2019; Xie et al., 2011; Yuan et al., 2022). The nucleotide diversity (Pi values) of 10 species was analyzed by sliding window analysis (Fig. 4). In the mitochondrial genomes of the 10 gomphocerine species, nucleotide diversity is highly variable. Nucleotide diversity ranges from 0.08 to 0.132, with higher nucleotide diversity in genes ND2, ND6, and ATP8, which are 0.132, 0.121, and 0.115, respectively. In contrast, the nucleotide diversity of COX1, ND1, and ND4L is lower, with 0.094, 0.084, and 0.080, respectively. This indicates that COX1, ND1, and ND4L are relatively conserved genes.

Ka/Ks indicates the ratio between the non-synonymous substitution rate (Ka) and the synonymous substitution rate (Ks) of two protein-coding genes, which can be used as an important marker to estimate the evolutionary rate. We calculated the Ka/Ks values of the mitochondrial genomes among 10 species of gomphocerine (Fig. 5; Table S4). The Ka/Ks values of all PCGs were less than one, indicating that these genes evolved under purifying selection and were evolutionarily conserved in the mitochondrial genome. ND6 had the highest Ka/Ks value, followed by ATP8 and ND5 and COX1 had minimum Ka/Ks value (Ka/Ks = 0.058) and low evolutionary rate, indicating that the COX1 gene had strong purifying selection and evolutionary conservation, which could be used as an important marker to identify relatedness among species, therefore, a partial fragment of COX1 is often used as DNA barcodes for inferring species phylogenetic relationships (Hebert, Ratnasingham & deWaard, 2003). In contrast, ND6 had the highest Ka/Ks value (Ka/Ks = 0.293), showing a faster evolutionary rate with less selection pressure in PCGs, which undergone relatively weak purifying selection. It can be used to assess intraspecific relationships and is more suitable as a potential molecular marker in population genetics (Yuan et al., 2021; Pu et al., 2022; Chen et al., 2021).

Phylogenetic analysis

The heterogeneity sequence divergence of the two matrixes PCG123 and PCG12 was assessed (Fig. 6), both indicated that the mitochondrial genomes of gomphocerine species showed lower heterogeneity (the similar scores for pairwise sequence comparisons were the lowest). Furthermore, among the species in Gomphocerinae, Orinhippus tibetanus shows higher heterogeneity than others. Substitution saturation of PCGs of 41 sequences were tested, and Xia’s analysis showed Iss < Iss.c and p < 0.05, revealed that base substitutions had not saturated and phylogenetic analysis could be performed.

We investigated the phylogenetic position of species among Gomphocerinae within Acrididae (Fig. 7). Tree topologies were consistent from both BI and ML analysis with high bootstrap values (BS) and bayesian posterior probability values (PP) in most clades. In Gomphocerinae, except for Orinhippus tibetanus, all other species showed significant monophyly (PP = 100, BS = 1). Orinhippus tibetanus and the species of Oedipodinae are recovered as sister groups, which is greatly supported (PPs = 1, BSs = 100). The results further support the reclassification of Orinhippus tibetanus as a member of the Oedipodinae using molecular systematics by Gao, Chen & Jiang (2017). In addition, by means of species descriptions of the genus and type species of the Orinhippus, previous research has classified the genus as a member of Locustinae (=Oedipodinae) (Uvarov, 1921; Bei-Bienko & Mishchenko, 1951), our results corroborate this ancient view from a molecular point of view. Therefore, we believe that antennal structure is of lesser importance as a higher order taxonomic character for grasshoppers. Moreover, this study clarifies the phylogenetic status of the genus Aeropedellus in Gomphocerinae, which is recovered a sister group with the genus Dasyhippus.

Phylogenetic analysis showed that the filiform antennae tend to be ancestral, while the condensed and expanded antennal flagellum which ends gradually to form a clubbed shape, appears to be more evolutionarily advanced (Fig. 8). This may first occur in micro-club-shape antennae of plateau species Pacris xizangensis, and the subsequently differentiated genera Myrmeleotettix and Aeropedellus show similar antennal morphology, in which the male antennal ends are slightly expanded at 5–7 segments, forming antennal ends that are twice as wide as long. Starting from Dasyhippus, the male antennal ends are very expanded at 7-8 segments, forming antennal ends that are significantly wider than long. Similar antennal morphology occurs in Gomphocerus and Gomphocerippus. The results indicated a possible antennal evolutionary trend in Acrididae, in which the filiform antennae are more ancient and gradually evolved into clubbed shape antennae.

Convergent evolution refers to independent lineages evolving similar phenotypes under similar selective pressures (Fraser & Whiting, 2019), but the phenomenon of convergent evolution is not easy to identify in evolution. Both Orinhippus tibetanus and Pacris xizangensis are distributed in Tibet, with very similar altitudes and environmental factors in their habitat, both species are classified in Gomphocerinae. However, based on phylogenetic analysis using mitochondrial genomes, Orinhippus tibetanus and Pacris xizangensis are on two independent clades. This study speculated that Orinhippus tibetanus and Pacris xizangensis may have been subjected to similar environmental selection pressures that formed similar antennal morphology convergently.

The phylogenetic relationships among Melanoplinae, Catantopinae, and Oxyinae are unclear in this tree, which may be the result of incomplete sampling. However, since they are not the focus of this research, no further discussion has been conducted. The low confidence level may be due to insufficient sampling, and in further studies, a wide range of sampling and multiple methods may be used to explore the phylogenetic relationships of these subfamily units in Acrididae. Further understanding of the triggering factors for evolution and the convergence of ecological forms in entire tree of life may help clarify genomic constraints and historical contingencies that have led to convergent evolution.

Divergence time estimation

The divergence time estimates (95% HPD) of each species based on the topology recovered from BEAST analysis were exhibited in Fig. 9. The divergence time tree indicated that Orinhippus diverged first, early in the Oligocene 33 Mya (32–35 Mya, 95% HPD), while the divergence events of the other clubbed antennae species occurred in late Miocene 19 Mya (15–23 Mya, 95% HPD). The divergence between micro-clubbed-shape and clubbed-shape antennae occurred between Aeropedellus and Dasyhippus, and before divergence of Gomphocerus and Gomphocerippus, which was about 8–11 Mya.

According to the BEAST analysis, the formation of club shape antennae occurred approximately 18 Mya in Pacris, after which there were two independent evolutionary events, resulting in extreme enlargement of the antennal ends (Gomphocerus and Dasyhippus) and slight enlargement (Myrmeleotettix and Aeropedellus). The clubbed shape antennae are an important taxonomic characteristic at the family level, which is occurred only in species of the Gomphocerinae (Yin, 1982; Yin & Xia, 2003). However, the divergence times indicate that the earliest clade of the gomphocerine grasshoppers, Pacris xizangensis, diverged at approximately 18 Mya (15–23 Mya, 95% HPD). Compared with the divergence times of subfamilies such as Oedipodinae, Acridinae, and Melanoplinae, the clubbed antennae grasshoppers have a relatively brief divergence history, may not to reach the family category, which confirms the current classification system that places them in Gomphocerinae under the tribe Gomphocerini (Otte, 1981). The results of the divergence time estimates reconfirmed that antennal morphology should be given less weight as a taxonomic character for grasshoppers higher-level classification. It also suggests that species within the subfamily Gomphocerinae are not monophyletic.

Limitations in selecting gene fragments, lateral gene transfer during early evolution, or ancestral polymorphisms may result in phylogenetic incongruence between morphological and molecular data during speciation events. Feng et al. (2022) suggested that constructing a species relationship tree based on only partial genes and phenotypes may not be reliable, and genome-wide data is the gold standard for reconstructing the evolutionary history of species. However, the genome of grasshoppers is significantly larger than other insects (Alfsnes, Leinaas & Hessen, 2017; Husemann et al., 2020), whole genome sequencing is expensive and it is difficult to analyze the data. Therefore, it is often more practical to select more conservative molecular markers to explore the true evolutionary relationships among species in Acrididae. In addition, incomplete lineage sorting and convergent evolution may also cause contradictions between morphology and molecular data. The taxonomic category of clubbed antennae grasshoppers in Acrididae requires deeper investigation using larger scale sampling.