Mitogenomics of Chinch Bugs from China and Implications for Its Coevolutionary Relationship with Grasses

Simple Summary Blissidae is a group with high species richness in Lygaeoidea, and most of them live in the leaf sheaths of Poaceae plants. Here, 10 new mitogenomes from 10 species of eight genera from Blissidae were sequenced and analyzed. Gene rearrangement is only found in Pirkimeru japonicus (PiGXBS1), which is formed as the duplication of tRNA-H. Coupled with published data, phylogenetic analyses and divergence time were performed in Blissidae. The divergence within Blissidae began about 56 million years ago (Ma), in which the genus level divergence was concentrated at 30–51 Ma, slightly later than the diversification of Poaceae. The consistency of the divergence time between Blissidae and Poaceae might hint at the coevolutionary relationship between them. Our study provides a valuable resource for understanding insect–host relationships. Abstract Blissidae (the Chinch bug) is a group with high species richness in Lygaeoidea, but there are only a few descriptions of mitochondrial genomes available. We obtained mitogenomes from 10 species of eight genera from Blissidae through second-generation sequencing technology. The length of the mitochondrial genome (excluding the control region) is between 14643 and 14385 bp; the content of AT is between 74.1% and 77.9%. The sequence of the evolution rate of protein coding genes was as follows: ND5 > ATP8 > ND6 > ND2 > ND4 > ND4L > ND1 > ATP6 > ND3 > COIII > COII > CYTB > COI. The mitogenomic structure of Blissidae is highly conservative. Gene rearrangement is only found in Pirkimeru japonicus (PiGXBS1), which is formed as the duplication of tRNA-H. The intergenic spacer between ND4 and tRNA-H, which form an obvious stem-and-loop structure, was found in all samples in this study. The phylogenetic trees generated by BI and ML indicated that Blissidae can be divided into three major clades: Clade A (only included Macropes); Clade B ((Pirkimerus + Bochrus) + Iphicrates); and Clade C ((Ischnodemus + Capodemus) + (Cavelerius + Dimorphopterus)). The divergence within the Blissidae began at about 56 Ma. At the genus level, the divergence was concentrated at 30–51 Ma, slightly later than the diversification of Poaceae. The consistency of divergence time between Blissidae and Poaceae might hint at the coevolutionary relationship between them, but further molecular and biological evidence is still needed to prove it.


Introduction
Mitogenomes are widely used in the study of molecular evolution, phylogeny, population genetics, and phylogeography [1][2][3][4]. The insect mitogenome is a double-stranded, closed circular DNA molecule. It is generally between 14 and 20 kb in length, mainly influenced by the length of the control region (CR) [5,6]. Insect mitogenomes have conserved genetic composition and arrangement similar to that of Drosophila yakuba (Diptera: Drosophilidae) [7,8]. The information of the mitochondrial gene structure has been considered as the key signal of evolutionary biology [2,[9][10][11][12].
Lygaeoidea is the second largest superfamily within Pentatomomorpha, currently comprising 14 families. In the published mitogenomes of Lygaeoidea, the gene structure is consistent with that of Drosophila yakuba [2,10,13,14]. As a group with relatively high species richness in Lygaeoidea, Blissidae has more than 420 species belonging to 55 genera, which are widely distributed in tropical and subtropical areas (Lygaeoidea Species File: http://lygaeoidea.speciesfile.org/Common/basic/Taxa.aspx?TaxonNameID=1208147, accessed on 1 June 2022). To [15,16]. Almost all chinch bugs are sap-feeding and live in leaf sheaths, with their bodies showing significant morphological specialization and extreme flattening with the adaptation to the living environment [17][18][19][20][21]. Poaceae species are the main host plants of Blissidae. In previous study, through the comparison of the host plants utilized by chinch bugs and the plant phylogenies, it was hypothesized that the Blissidae species might be radiated and diversified in the Poaceae [22]. However, research on the molecular biology of chinch bugs is very scarce, and research on its internal phylogenetic relationship has never been carried out. The study of the internal phylogeny and evolutionary history of Blissidae will deepen our understanding of the family, and also provide important molecular data from species for the discussion of its relationship with host plants.
In this study, we obtained 10 species of eight genera from Blissidae (Table S1) and analyzed the characteristics of the mitogenome of Blissidae. At the same time, the phylogeny and divergence time of Blissidae was preliminarily studied and its relationship with the host plants was discussed.

Sample Collection
We obtained 10 species of Blissidae in China (Table S1). All specimens were preserved in 100% ethanol immediately in the field and stored at −20 • C before DNA extraction. The specimens were deposited in the Insect Molecular Systematics Lab, Institute of Entomology, College of Life Sciences, Nankai University, Tianjin, China (NKUM). The specimens were identified based on their morphological characteristics.

DNA Extraction, Sequencing, and Assembling
Genomic DNA was extracted from the mid-leg using a DNeasy Blood & Tissue Kit (QIAGEN). Sequencing was performed by Novogene (Tianjin, China) with an insert size of 250 bp and a pair-end 150 bp sequencing strategy on the Illumina platform. The mitogenomes were assembled in MitoZ [23] and IDBA-master [24] methods for assembly and mutual verification.

Annotation of Mitogenome
We use the Mitos Web Server (http://mitos.bioinf.uni-leipzig.de/index.py/, accessed on 1 June 2022) [25] to identify the boundaries of tRNA genes and the secondary structures of tRNAs. The start and stop codons of the protein-coding genes (PCGs) are determined by ORF Finder using invertebrate mitochondrial genetic codes, which are implemented by the NCBI website (https://www.ncbi.nlm.nih.gov/orffinder/, accessed on 1 June 2022). The rRNA boundaries are predicted by comparison with homologous regions of the other published lygaeoidea mitogenomes in GenBank.

Phylogenetic Inference and Divergence Time Estimation
We carried out phylogenetic analyses based on the mitochondrial dataset of 18 samples, including 10 samples of Blissidae obtained in this study, and eight outgroups were selected (four of them belonging to Lygaeoidea, and the other four belonging to Pyrrhocoridea) (Tables S1 and S2). Alignments of all PCGs with other mitogenomes were performed based on their amino acid sequences using MUSCLE implemented in MEGA X [26]. The rRNAs were aligned using MAFFT v7 with the option -G-INS-I [27]. All individual genes were Insects 2022, 13, 643 3 of 8 then concatenated as a combined matrix. The PCG12RT matrix, including all 13 PCGs with third codon positions removed, two rRNA and all tRNA genes, was used for the final phylogenetic analysis.
The phylogenetic analyses were conducted utilizing Bayesian inference (BI) and maximum likelihood (ML). PartitionFinder 2 [28] was used to determine the optimal partitioning strategies and the best-fit nucleotide substitution model. The results of PartitionFinder for MrBayes v3.2.5 [29] and IQ-TREE v2 [30] revealed the best-fit nucleotide substitution model in Table S3. Under the BI method, we used Phylobayes-MPI v.1.5a [31] and MrBayes 3.2.5 [29]. Regarding MrBayes analysis, a Markov chain Monte Carlo analysis (four chains) was run for 10,000,000 generations, and samples were recorded every 1000 generations. The first 25% of the samples were discarded as burn-in, and the remaining samples were used to summarize the Bayesian posterior probabilities (BPP). In the Phylobayes analysis, two independent Markov chain Monte Carlo chains were run after the removal of constant sites from the alignment and were stopped after the two runs had satisfactorily converged (maxdiff < 0.1). A consensus tree was computed from the remaining trees combined from two runs after the initial 25% trees of each run were discarded as burn-in. ML analyses of molecular dataset were conducted with IQ-TREE v2 [30], and the node support values were assessed by bootstrap resampling (BP) calculated using 1000 replicates. Finally, FigTree v 1.3.1 [32] was used to visualize the phylogenetic tree.
MCMCTree from the Phylogenetic Analysis by Maximum Likelihood (PAML) 4.9 package [33] was used to estimate divergence times in Blissidae with a relaxed molecular clock based on the ML tree of PCG12RT. The earliest fossil of Blissidae, Eoblissus gallicus (58.3-53.7 Ma), and the fossil of Berytidae, Metacanthus serratus (33.9-28.4 Ma), were used in the calculation. The main parameters during operation are as follows: model = 4, rgene_gamma = 2, 20. After a 20,000 burn-in, the MCMCTree program was run for 1,000,000 MCMC steps and sampled every 100. The robustness of the MCMCTree results was checked by comparing the consistency of at least two independent runs, with all parameters at least 200 for the effective sample sizes.

Mitogenome Comparison and Rearrangement
In the 10 mitogenomes obtained in this study from Blissidae, the length of the mitochondrial genome (excluding the control region) is between 14643 and 14385 bp; the content of AT is between 74.1 and 77.9% (Figures 1 and 2). The nucleotide skew statistics for the whole genome obviously showed AT-skew and CG-skew; the AT content of rRNA was significantly higher than that of PCGs and tRNA (Table 1). For PCGs, the ratios of the nonsynonymous nucleotide changes (Ka) versus the synonymous nucleotide changes (Ks) were all below 1, indicating that they evolved according to the purifying selection [34]. The sequence of the evolution rate of the protein-coding genes was as follows: (Table S4).
The subregion of the intergenic spacer between the tRNA-H and the ND4 was found in all samples in our study ( Figure 3A). Through structural prediction, we found that the subregion of the intergenic spacer formed a stem-and-loop structure. The mitogenome structure was consistent and identical to the putative ancestral arrangement of insects [35] among all of the species we obtained (expect the sample Pirkimeru japonicus), showing the highly conservative mitogenome structure of Blissidae. It is noteworthy that gene rearrangement of tRNA-H replication was found in Pirkimeru japonicus (PiGXBS1) (Figures 1, 2 and 3B). The second tRNA-H was located between tRNA-M and ATP8, and the distance between the two tRNA-H was large (Figures 2 and 3B). Due to the consistent neck ring structure at the front ends of the two tRNA-H, we inferred that the rearrangement was due to the recent recombination caused by the stem-and-loop structure [10,36,37].

Phylogenetic Analyses and Divergence Time Estimation
The phylogenetic trees generated by BI and ML were highly consistent. The 10 species could be divided into three major clades: Clade A (only including Macropes); Clad B ((Pirkimerus + Bochrus) + Iphicrates); and Clade C ((Ischnodemus + Capodemus) + (Cavelerius + Dimorphopterus)) ( Figure 4). The monophyly of Macropes was supported by the phylogenetic results. The divergence within Blissidae began at about 56 Ma, and the divergence at the genus level was concentrated at 30-51 Ma. Previous research showed that Poaceae plants were the main hosts of the species included in the study, including Bambusoideae, Pooideae, Oryzoideae, Panicoideae, Chloridoideae, and Arundinoideae [22]. Recent studies revealed that the diversification within these subfamilies of Poaceae was about 66-54 Ma [38], slightly earlier than the diversification of species in Blissidae. The species of Blissidae live in leaf sheaths and are highly dependent on the host plants; we speculated that Blissidae species might have evolved with the continuous radiation of the host after the diversification of Poaceae. However, at present, the speculation about the co-evolution between Blissidae and Poaceae was very preliminary through the consistency of time in this study, and we still need more subsequent evidence from biology and molecular science to prove it. Meanwhile, Clade B contained three species, with bamboo as the host plant, indicating that its monophyly might be involved with host adaptation. Moreover, in combination with base composition analysis, we found that there were significant differences in AT skew between different clades, and its biological significance needs to be further discussed with more genomic data.
At present, there are still some problems in the reconstruction of the phylogeny of Blissidae, especially with the lack of known Blissidae molecular data. Therefore, more molecular data need to be acquired to provide a more comprehensive view of the evolution of Blissidae and Lygaeoidea. At the same time, we should pay more attention to the relevant host information of Blissidae species, and further reveal the evolutionary history between them through more molecular and biological evidence in the future. The divergence within Blissidae began at about 56 Ma, and the divergence at the genus level was concentrated at 30-51 Ma. Previous research showed that Poaceae plants were the main hosts of the species included in the study, including Bambusoideae, Pooideae, Oryzoideae, Panicoideae, Chloridoideae, and Arundinoideae [22]. Recent studies revealed that the diversification within these subfamilies of Poaceae was about 66-54 Ma [38], slightly earlier than the diversification of species in Blissidae. The species of Blissidae live in leaf sheaths and are highly dependent on the host plants; we speculated that Blissidae species might have evolved with the continuous radiation of the host after the diversification of Poaceae. However, at present, the speculation about the co-evolution between Blissidae and Poaceae was very preliminary through the consistency of time in this study, and we still need more subsequent evidence from biology and molecular science to prove it. Meanwhile, Clade B contained three species, with bamboo as the host plant, indicating that its monophyly might be involved with host adaptation. Moreover, in combination with base composition analysis, we found that there were significant differences in AT skew between different clades, and its biological significance needs to be further discussed with more genomic data.
At present, there are still some problems in the reconstruction of the phylogeny of Blissidae, especially with the lack of known Blissidae molecular data. Therefore, more molecular data need to be acquired to provide a more comprehensive view of the evolution of Blissidae and Lygaeoidea. At the same time, we should pay more attention to the relevant host information of Blissidae species, and further reveal the evolutionary history between them through more molecular and biological evidence in the future.

Conclusions
In this study, 10 new mitogenomes of eight genera from Blissidae were sequenced and analyzed. Coupled with published data, phylogenetic analyses and divergence time were performed in Blissidae. A gene rearrangement with the duplication of tRNA-H was identified within Pirkimeru japonicus, which was the first instance of mitochondrial gene rearrangement discovered in Lygaeoidea. The divergence within Blissidae began at about 56 Ma, and the divergence at the genus level was concentrated at 51-30 Ma, slightly Insects 2022, 13, 643 7 of 8 later than that of the host plant Poaceae. The consistency of the divergence time between Blissidae and Poaceae might hint at the coevolutionary relationship between them, but further molecular and biological evidence is still needed to prove it.
Supplementary Materials: The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/insects13070643/s1, Table S1: Experimental sample information of Blissidae; Table S2: Public data used in phylogenetic analysis; Table S3: Best-fit models of sequence evolution and partitioning schemes selected by PartitionFinder for phylogenetic reconstructions; Table S4: The nonsynonymous nucleotide changes (Ka) and the synonymous nucleotide changes (Ks) of PCGs.