Complete Mitochondrial Genomes and Evolutionary Insights of Two Commercially Farmed Edible Crickets (Gryllus bimaculatus and Teleogryllus mitratus) from Thailand

1. Introduction

Edible insects are increasingly recognized as sustainable alternatives to conventional livestock due to their high nutritional value, efficient feed conversion, and lower environmental footprint [1,2]. In Thailand, cricket farming has expanded rapidly and is now a major component of the edible insect industry, providing income for rural communities while supplying domestic and export markets. Two important species, the two-spotted cricket Gryllus bimaculatus and Teleogryllus mitratus, are widely reared and traded, but reliable genetic resources remain limited despite their commercial importance [3,4,5,6,7] Because farm populations can experience founder effects, selective propagation, and mixing of stocks, genetic reference data are essential for accurate species authentication, monitoring of maternal lineages, and improvement of breeding and traceability across production chains [8].

Even with their commercial importance, comprehensive genomic resources for these cricket species remain surprisingly limited. This gap in molecular data constrains our understanding of their evolutionary history, population genetics, and the genetic foundations of traits crucial for farming optimization and environmental adaptation. Previous phylogenetic studies have relied primarily on partial mitochondrial markers or nuclear genes, approaches that often yield inconsistent results or lack the resolution needed for robust evolutionary inferences [9,10]. Mitogenomes are widely used for phylogenetic inference, population structure analyses, and molecular identification due to their relatively conserved gene content, maternal inheritance, and comparatively rapid mutation rate [11,12].

The mitochondrial genome (mitogenome) is a circular DNA molecule typically ranging from 14 to 20 kb in animals and encodes 13 PCGs involved in oxidative phosphorylation (OXPHOS), 22 tRNAs, two rRNAs, and a non-coding control region that participates in replication and transcription [13,14,15]. In insects, mitogenome sequencing has also contributed to resolving taxonomic ambiguities and understanding evolutionary relationships across diverse lineages [16,17]. For edible insects and farmed species, complete mitogenomes can provide standardized markers for traceability, detection of mislabeling or contamination, and comparative analyses of mitochondrial function relevant to metabolic performance under rearing conditions [18,19,20].

Despite the growing economic significance of edible crickets, mitogenome resources remain limited for several species used in commercial farming in Southeast Asia. Previous studies have reported mitogenomes for selected Gryllidae species, but farm-derived reference genomes for Thai production systems are still scarce. Moreover, mitochondrial gene evolution, codon usage bias, and selective pressures have not been sufficiently explored in these farmed crickets, especially in relation to their metabolic genes underlying OXPHOS. As farm environments can impose distinct thermal, dietary, and density-related conditions compared to the wild, mitochondrial function may be an important component of physiological adaptation and performance, although direct evidence requires integration with ecological and nuclear-genomic data [21].

In this study, we sequenced and annotated the complete mitochondrial genomes of G. bimaculatus and T. mitratus collected from commercial farms in Thailand. We characterized genome organization, nucleotide composition, codon usage patterns, and tRNA structures. We also assessed selective pressures acting on mitochondrial PCGs using codon-based evolutionary models, and reconstructed phylogenetic relationships within Gryllidae using concatenated mitochondrial genes. By focusing on farm-derived stocks, this work provides baseline genetic references for commercially important edible crickets in Thailand and supports future efforts in traceability, genetic monitoring, and integrative breeding research.

2. Materials and Methods

2.1. Sample Collection and DNA Extraction

Cricket specimens of G. bimaculatus and T. mitratus were obtained from commercial farms in Phitsanulok and Phayao provinces, Thailand. Fresh specimens were immediately placed in labeled plastic containers and transported to the laboratory, where they were anesthetized before being stored at −20°C until processing. Total genomic DNA was extracted using the standard phenol-chloroform method [22], following established protocols for arthropod tissue. DNA quality and concentration were assessed using spectrophotometry and gel electrophoresis before submission to Macrogen, Inc. (South Korea) for library preparation and high-throughput sequencing.

2.2. Genome Sequencing and Assembly

Paired-end DNA libraries were prepared using the Illumina TruSeq™ Nano DNA Prep Kit and sequenced on an Illumina platform to generate 151 bp paired-end reads. Raw sequencing reads were evaluated for quality using FastQC v0.11.9 [23]. The adapters and low-quality bases were trimmed from raw reads using Trimmomatic v0.39 [24], with the following parameters: LEADING:3, TRAILING:3, SLIDINGWINDOW:4:20, and MINLEN:20. Post-trimming quality was reassessed to ensure successful preprocessing.

High-quality trimmed reads were assembled using NOVOplasty v4.3.1 [25], a de novo organellar genome assembler. We used partial sequences of the cytochrome c oxidase subunit I (cox1) gene from closely related cricket species as seed sequences to initiate the assembly process. Critically, k-mer length was optimized post-hoc through preliminary assembly trials: we selected k = 29 for G. bimaculatus and k = 19 for T. mitratus, reflecting species-specific sequence complexity and GC content that modulate optimal k-mer sizing. The expected genome size of 15-16 kb, informed by comparative analysis of published Gryllidae mitogenomes, further constrained the assembly search space, reducing erroneous contig formation while accelerating convergence.

Complete mitochondrial genomes were annotated using MITOS v1.1.5 [26] with default settings for invertebrate mitochondrial genetic code. Protein-coding genes, ribosomal RNA genes, transfer RNA genes, and control regions were identified and their boundaries refined manually. We validated annotations by comparing predicted sequences against the NCBI nucleotide database using BLAST [27]. The annotated mitochondrion sequences of G. bimaculatus and T. mitratus were submitted into NCBI databases. The GenBank files obtained from NCBI databases were utilized to visualize gene organization and arrangement via Organellar Genome DRAW (OGDRAW) v1.3.1 [28]. The base composition, codon usage and relative synonymous codon usage (RSCU) of protein-coding genes (PCGs) were analyzed using PhyloSuite v1.2.3 [29]. The non-synonymous/synonymous substitution ratios (Ka/Ks, ω) were calculated using HyPhy (Datamonkey) using the Single-Likelihood Ancestor Counting (SLAC) model to characterize selective pressures, while using the Branch-site Unrestricted Statistical Test for Episodic Diversification (BUSTED) model to detect gene-wide positive selection [30,31,32,33,34,35,36]. The secondary structures of tRNAs were interfered using the MITOS v1.1.1 on GALAXY server [26].

2.3. Phylogenetics Analysis

To investigate evolutionary relationships within Gryllidae, we compiled a dataset including our two sequenced mitochondrial genomes and 23 previously published Gryllidae mitogenomes retrieved from GenBank, while Ceuthophilus sp. (accession no. OR551732) and Gryllotalpa henana (accession no. NC_071757) served as outgroups based on their established phylogenetic position as sister groups to Gryllidae. (See Supplementary Table S1).

For phylogenetic reconstruction, thirteen protein-coding genes (PCGs) and two ribosomal RNA genes (12S and 16S rRNA) were extracted from each mitochondrial genome. The sequences were aligned, trimmed, and concatenated using PhyloSuite v1.2.3 [29], with individual alignments of each PCG and rRNA performed in MAFFT v7 [37]. Maximum Likelihood (ML) analyses were performed in IQ-TREE v2.2.0 [38], with model selection determined by PartitionFinder2 v2.1.1 [39] under the AIC criterion using the greedy search algorithm, and nodal support assessed with 1,000 ultrafast bootstrap replicates. Bayesian Inference (BI) was performed in MrBayes v3.2.7 [40] using two independent Markov Chain Monte Carlo (MCMC) runs of 1,000,000 generations, sampling every 100 generations and discarding the first 25% as burn-in. The resulting phylogenetic trees were visualized and edited in iTOL v7.2.1. [41].

3. Results

Assembly of high-depth sequencing reads yielded complete, circular mitochondrial genomes for both species, with coverage depths (32,391× and 63,258× for G. bimaculatus and T. mitratus, respectively) substantially exceeding thresholds required for confident base-calling and variant detection (GenBank accessions: PP230540 and PP297527). The resulting genomes, at 15,955 bp and 16,046 bp respectively, align with the expected size distribution for Gryllidae mitogenomes (typically 14–17 kb), suggesting successful capture of complete mitochondrial DNA with minimal loss of genetic material during extraction and sequencing.

As expected from the conserved nature of mitochondrial gene content across animals, both species display the canonical 37-gene architecture characteristic of insects, with composition and organization fully consistent with previously sequenced Gryllidae mitogenomes. Both mitochondrial genomes exhibit the typical insect gene organization, containing 37 genes: 13 protein-coding genes (PCGs), 22 transfer RNA (tRNA) genes, and 2 ribosomal RNA (rRNA) genes, plus a control region, consistent with the metazoan mitochondrial genetic code (Figure 1).

These two mitochondrial genomes exhibit the typical insect gene content and organization, consisting of 13 protein-coding genes (PCGs) atp6, atp8, cytb, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4l, nad5, and nad6 along with 22 transfer RNA (tRNA) genes, including two copies each of tRNA^Leu and tRNA^Ser, and two ribosomal RNA (rRNA) genes (rrnS and rrnL).

Gene organization adheres to the asymmetric strand distribution typical of insects, with both species encoding 4 PCGs and 11 tRNAs on the minority (N) strand alongside both rRNAs, while 9 PCGs and 11 tRNAs occupy the majority (J) strand. This division reflects underlying constraints related to concurrent transcription and replication, wherein genes on opposite strands experience different mutational pressures [42]. At the codon level, most PCGs initiate with standard ATN start codons; however, cox1 and nad1 deviate from this pattern, employing alternative start codons (CAA and TTG in G. bimaculatus; CGA and TTG in T. mitratus). Such non-canonical initiation codons, increasingly documented in Orthoptera [42,43], likely reflect the permissive nature of mitochondrial translation systems and may reduce ribosome stalling relative to standard AUG initiation in the AT-rich mitochondrial genome. Termination patterns showed predominant use of the complete TAA stop codon, though cox1, cox2, nad3, and nad4 in both species, plus cytb in T. mitratus, terminate with incomplete stop codons (T alone), a common feature in arthropod mitochondrial genomes [44,45] (Supplementary Table S2).

Nucleotide composition analysis revealed strong AT bias in both species, with adenine and thymine comprising 74.0% and 73.4% of total nucleotides in G. bimaculatus and T. mitratus, respectively. This AT-rich composition varied among genomic regions: protein-coding genes showed the highest AT content (73.4-74.9%), followed by rRNAs (73.5-74.9%) and tRNAs (75.5-76.0%). Strand asymmetrical analysis using AT-skew and GC-skew values revealed distinct patterns across genomic regions The complete mitochondrial genomes exhibited positive AT-skew (0.095 in G. bimaculatus, 0.082 in T. mitratus) and negative GC-skew (-0.31 and -0.28, respectively), indicating strand-specific mutational biases. Notably, PCGs displayed contrasting patterns with negative AT-skew (-0.138 to -0.141) but near-neutral GC-skew, while rRNAs showed negative AT-skew but strongly positive GC-skew (0.365-0.398). Transfer RNAs exhibited intermediate patterns with near-neutral AT-skew and moderately positive GC-skew values (Table 1).

Analysis of relative synonymous codon usage (RSCU) across all 13 PCGs revealed a strong codon bias toward A/T-ending codons in both species (Figure 2 and Supplementary Table S3). Of the 62 sense codons analyzed, both species showed a marked preference for NNU and NNA codons, with UUA (Leu2) displaying the highest RSCU values (3.97 in G. bimaculatus, 3.75 in T. mitratus) indicating that it had undergone strong selection during evolution. In contrast, the RSCU values for CUG (Leu1) which end with G and C, were relatively low. Notably, Leu2 and Ser2 consistently had the highest RSCU values in both species.

Amino acid composition analysis revealed highly conserved usage patterns between species (Figure 3 and Supplementary Table S3). Overall, both species showed highly conserved usage patterns. The most frequently used amino acid was Leucine2 (UUR codons), with counts of 438 in G. bimaculatus and 433 in T. mitratus. This was followed by Leucine1 (CUN codons) and Phenylalanine, which were also abundant (Leu1: 381 and 375; Phe: 333 and 324, respectively). In contrast, the least utilized amino acids were Cysteine (49–51) and Tryptophan (101–102). Several amino acids, including Lysine, Glycine, and Serine2, also occurred at relatively high frequencies. Minor interspecific differences were observed, with T. mitratus exhibiting slightly higher counts of Serine2 (238 vs. 226) and Glutamine (81 vs. 80), while Methionine and Proline were marginally more common in G. bimaculatus.

Differential selective pressures among mitochondrial genes have important implications for understanding metabolic evolution and identifying candidate genes for functional studies. To assess these pressures and identify genes potentially under adaptive evolution, we employed codon-based substitution rate analyses. The nonsynonymous to synonymous substitution rate ratio (Ka/Ks, ω) of mitochondrial protein-coding genes (PCGs) were exhibited that all PCGs, mean ω values were below 1, indicating that mitochondrial genes are predominantly evolving under pervasive purifying selection. Consistent with this, the BUSTED analysis detected gene-wide evidence of episodic diversifying (positive) selection in cox1, cox3, cytb, and nad5, as indicated by significant likelihood ratio tests (P < 0.05; marked by asterisks). Among the PCGs, atp8 showed a relatively elevated ω, consistent with weaker functional constraints and a faster evolutionary rate compared with other mitochondrial genes. Overall, these results demonstrate heterogeneity in selective rules among mitochondrial genes and respiratory chain complexes, with pervasive purifying selection dominating the mitochondrial genome despite episodic diversifying selection acting on specific genes. (Figure 4 and Supplementary Table S4).

All 22 mitochondrial tRNAs were successfully identified and annotated, with most displaying canonical cloverleaf secondary structures containing the acceptor stem, dihydrouridine (DHU) arm, anticodon arm, TΨC arm, and variable loop (Supplementary Figures S1 and S2). The notable exception was trnS1 (AGN), which lacks the DHU arm in both species, a characteristic feature widely observed in metazoan mitochondrial genomes that does not impair functionality (Figure 5).

Both species contain complete large and small ribosomal RNA subunits positioned in the typical arthropod arrangement: lrRNA between tRNA^Leu and tRNA^Val, and srRNA between tRNA^Val and the control region. The lrRNA genes measured 1,313 bp (G. bimaculatus) and 1,315 bp (T. mitratus), while srRNA genes spanned 762 bp and 808 bp, respectively. These size variations fall within the normal range for Gryllidae species are consistent with minor indel events in non-critical structural regions.

The phylogenetic analyses recovered well-supported clades consistent with established genera within the family Gryllidae. The phylogenetic trees generated by both ML and BI methods exhibited identical topologies (topological congruence), differing only in minor branch lengths. All three focal genera, Gryllus, Teleogryllus, and Tarbinskiellus were recovered as monophyletic groups with high statistical support (ML bootstrap ≥ 95%; BI posterior probability ≥ 0.95). The two newly sequenced species were placed firmly within their respective genera with maximum support: G. bimaculatus clustered with other Gryllus species, while T. mitratus grouped closely with T. emma and T. infernalis, confirming their evolutionary relationships. The tree was rooted using Ceuthophilus sp. and G. henana as outgroups, which were clearly separated from the ingroup taxa, providing additional confidence in the inferred topology. The congruent results from both ML and BI analyses, combined with consistently high support values, demonstrate the effectiveness of mitochondrial genome data for resolving phylogenetic relationships within edible cricket lineages (Figure 6 and Supplementary Table S5).

4. Discussion

This study presents the first complete mitochondrial genomes of G. bimaculatus and T. mitratus sourced from commercial farms in Thailand, substantially expanding genomic infrastructure for one of Southeast Asia's most economically significant edible insect industries. Beyond documenting genome sequences, our comparative analyses reveal heterogeneous evolutionary dynamics across mitochondrial genes, with implications for understanding how farmed populations may diverge genetically from wild ancestors and how mitochondrial function might be optimized through selective breeding. The overall gene content and organization of both mitogenomes conform to the conserved arthropod pattern, consistent with strong functional constraints on mitochondrial genome architecture across Gryllidae [46,47,48]. However, the presence of non-canonical start codons (e.g., CAA, CGA, TTG) and incomplete stop codons (T or TA), represents a common adaptation in insect mitochondrial genomes that likely reflects the specialized nature of mitochondrial translation systems [42,43,49,50,51].

Comparison with previously reported G. bimaculatus mitogenomes (accession number: NC_053546 and MZ440656) [47,52] confirmed overall structural conservation, with identical gene organization and comparable AT-richness across populations. However, comparing newly assembled G. bimaculatus mitogenome in present study (PP230540.1) with the published sequence MZ440656 shows a difference in the D-loop repeat array. Our assembly contains two complete 219 bp repeats in tandem followed by a 195 bp partial copy, whereas MZ440656 has one complete repeat followed by the same partial copy. The extra repeat is supported by SRA BLAST of short-read datasets DRX261895 and DRX261894 and long-read dataset DRX261898, as well as by high read depth across the region and high mitogenome coverage analysis (32,391x). In this context, the D-loop repeat polymorphism is not only a technical detail but also a potentially informative marker class for maternal-lineage tracking and intraspecific differentiation [53].

Both species exhibit pronounced AT richness, with genome-wide AT content exceeding 73–74%, consistent with broad insect patterns and Gryllidae-specific comparative work [42,54,55]. The combination of positive AT-skew and negative GC-skew indicates strand-asymmetric substitution pressures that are plausibly linked to replication and transcription associated mutational biases [56]. Importantly, this compositional landscape provides a parsimonious explanation for the observed codon-usage bias toward A/U-ending codons and the high representation of leucine and phenylalanine residues, patterns also typical of Orthopteran mitogenomes [48,57,58]. While such biases are often attributed primarily to mutational pressure, their persistence in functionally constrained OXPHOS genes underscores that selection and translational constraints may also contribute to maintaining efficient mitochondrial protein synthesis [16,17,59].

The selective pressure analyses revealed pronounced heterogeneity in evolutionary dynamics among mitochondrial protein-coding genes. Although pervasive purifying selection predominated across the mitochondrial genome, as reflected by mean ω values below one for all PCGs, several genes exhibited signatures of episodic diversifying selection. In particular, cox1, cox3, cytb, and nad5 showed significant gene-wide evidence of positive selection. This suggests that even under the overarching regime of purifying selection, certain sites in these respiratory complexes have undergone adaptive evolution, possibly in response to specific metabolic demands [60]. These genes encode key components of the oxidative phosphorylation pathway, indicating that selective fine-tuning of respiratory chain function may play an important role in mitochondrial evolution [61], suggest that the selected changes may reflect adaptation to shared environmental pressures. The intermittent nature of positive selection suggests that these pressures may be fluctuating over time rather than constant, with only certain lineages or time periods experiencing directional selection on ETC function [62]. Future work incorporating functional data on thermal tolerance, metabolic rate, and feed efficiency across natural and farmed populations could test whether positive selection signatures correlate with physiological performance differences.

In contrast, atp8 displayed a relatively elevated ω compared with other mitochondrial genes, indicating substantially weaker functional constraint and a faster evolutionary rate compared with other PCGs. This pattern, consistently observed across Orthoptera [47,63,64], likely reflects the reduced essentiality of ATP synthase subunit 8 (atp8) is a small protein (~49 amino acids in crickets) that may be partially redundant or whose sequence requirements are less stringent than larger catalytic subunits.

The elevated ω may indicate that atp8 is undergoing directional selection for altered properties (e.g., changes in membrane topology or subunit-subunit interactions) [65], though distinguishing adaptive evolution from mutation-driven relaxation of constraint requires additional functional data.

The secondary structures of mitochondrial tRNAs in both G. bimaculatus and T. mitratus conform to the classical cloverleaf model, comprising the acceptor stem, dihydrouridine (DHU) arm, anticodon arm, variable loop, and TΨC arm. The sole exception is trnS1 (AGN), which notably lacks the DHU arm a feature widely observed across metazoan mitochondrial genomes, particularly in insects [54,66,67,68]. Notably, this recurrent structural reduction is widely treated as a lineage-consistent feature of metazoan mitochondrial tRNAs rather than a loss of function, and its presence here is therefore consistent with accurate annotation and expected mitochondrial tRNA evolution [55,66].

The robust phylogenetic resolution achieved here using mitochondrial data alone represents a considerable strength for applied purposes (species authentication, population tracking) but also highlights an important limitation of mitochondrial-only phylogenetics: the single-locus nature of mtDNA means that the recovered topology, while well-supported, represents the evolutionary history of maternal lineages specifically [69,70]. In farmed populations, where animals may be crossed between diverse maternal and paternal genetic backgrounds [3,71], the mitochondrial phylogeny may not fully capture nuclear-genomic structure or the extent of admixture among populations. Nonetheless, the monophyly of Gryllus and Teleogryllus and the strong support for species-level relationships provide confidence that these genera represent genuine evolutionary lineages and not unrecognized species complexes, an important consideration given the economic value of reliable species identification.

The genomic resources presented here enable several practical innovations in cricket agriculture. First, the complete genome sequences permit design of species-specific molecular assays (real-time PCR, high-resolution melt analysis) for rapid, accurate species identification [72], a critical quality-control metric given the market premium for G. bimaculatus and T. mitratus relative to less-valued species. Second, the identification of polymorphisms in the mitochondrial D-loop (our assembly documents copy-number polymorphism in repeat arrays not present in published genomes) provides a foundation for developing population-specific markers to track maternal lineages through breeding programs, enabling better documentation of genetic management and informed selection for traits related to growth rate, feed conversion, or fecundity. Third, the genes identified under positive selection, particularly ETC genes; represent empirically informed candidate loci for future functional validation studies aimed at understanding mitochondrial contributions to physiological traits under farmer optimization. Lastly, understudied application involves quality assurance in cricket farming supply chains. Because crickets are marketed both alive and as processed products (whole-insect meals, protein powders), molecular markers enabling traceability from farm to consumer could authenticate 'origin of origin' claims and detect species substitution or contamination with wild insects or unrelated species [18,19,20]. Mitochondrial barcoding, leveraging the complete genomes presented here, offers a rapid, cost-effective pathway to implement such traceability, with particular relevance for premium markets in Thailand and the EU where food authentication and supply-chain transparency are increasingly demanded by consumers and regulators.

5. Conclusions

This study provides complete mitogenomes for G. bimaculatus and T. mitratus, characterizing the molecular adaptation of Thai farmed crickets. Our findings reveal that while purifying selection preserves OXPHOS architecture, episodic positive selection drives respiratory adaptation to farm-related stressors. These resources facilitate molecular traceability for quality assurance and provide physiological targets for selective breeding. Future research should integrate wild-population genomics and phenotypic data to develop evidence-based breeding programs, ensuring the long-term sustainability and performance of insect agriculture.