Comparative analysis of codon usage patterns and phylogenetic implications of five mitochondrial genomes of the genus Japanagallia Ishihara, 1955 (Hemiptera, Cicadellidae, Megophthalminae)

Sample preparation and DNA extraction

Collection information for the specimens used in this study was provided in Table S2. Live specimens were collected by sweep nets and immediately preserved in 100% ethanol and stored at −20 °C until identification and DNA extraction. Samples were identified by their morphological characteristics (Rakitov, 1998; Dietrich, 2005; Viraktamath, 2011). Total genomic DNA was extracted from the whole body of adult males using a DNeasy©Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocols. Genomic DNA were stored at −20 °C, and male external genitalia were preserved in glycerol at room temperature. Voucher adult specimens with male genitalia and DNA sample have been deposited at the Institute of Entomology, Guizhou University, Guiyang, China (GUGC).

Sequence assembly and annotation

Genomes of the five species were sequenced by next-generation sequencing using Illumina HiSeq 4000 (Berry Genomic, Beijing, China). Each mitogenome was assembled using Geneious Primer (v.2019.2.1) (Kearse et al., 2012) using J. spinosa (GenBank NC_035685) (Wang et al., 2017) as the reference.

The assembled mitogenomes were annotated using MITOS web server with the invertebrate genetic code (Bernt et al., 2013), compared with mitogenomes of other leafhoppers, and identified by BLAST searches on NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi) to confirm that the sequences were correct (Altschul et al., 1997). The 13 protein-coding genes (PCGs) boundaries were identified by ORF Finder of Geneious Prime applying the invertebrate mitochondrial genetic code. Abnormal start and stop codons were determined through comparisons with other insects. The locations and secondary structures of 22 tRNA genes were estimated using tRNAscan-SE version 1.21 (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe & Eddy, 1997) and ARWEN version 1.2 (Laslett & Canbäck, 2008). The location of 16S rRNA was determined according to the location of trnL2 and trnV. The target gene was compared with the existing mitochondrial genome sequence by BLAST to find out the location of 12S rRNA and control region.

Mitochondrial genomic composition

The compositional properties of the 13 protein-coding genes in the mitochondrial genome of five Japanagallia species, including all nucleotides composition (A, T, G and C), nucleotides composition at the third position of codons (A3, C3, T3 and G3), all GC and AT contents, the average of nucleotide G and C present at first and second positions of codon (GC12), GC contents at the first, second, third position (GC1, GC2 and GC3) in percentages, were calculated by MEGA6.06 software (Tamura et al., 2013). Meanwhile, in order to understand the bias of nucleotide, we calculated skew value (GC and AT) (Tillier & Collins, 2000) using the following formula:

GC skew = $\frac{\mathrm{G}-\mathrm{C}}{\mathrm{G}+\mathrm{C}}$ and AT skew = $\frac{\mathrm{A}-\mathrm{T}}{\mathrm{A}+\mathrm{T}}.$

Relative synonymous codon usage

The relative synonymous codon usage (RSCU) definition refers to the relative probability between synonymous codons encoding corresponding amino acids for a specific codon (Sharp, Tuohy & Mosurski, 1986). If the value of RSCU > 1.6, the codon is over-represented, and if RSCU < 0.6, the codon is under-represented (Wong et al., 2010). When the RSCU value greater than 1 or less than 1 represent two types of codons, namely high bias (>1) and low bias (<1). However, when RSCU = 1, we consider that the codon usage lacks bias (Sharp & Li, 1986). The RSCU values of each codon of 13 protein-coding genes in the mitochondrial genome were calculated in the software MEGA6.06 (Tamura et al., 2013). Based on the relative synonymous codon values of 13 protein-coding genes of five newly sequenced mitogenomes, a heat map was drawn through TBtools v1.120 (Chen et al., 2020)

Neutrality plot

To understand the effects of mutational pressure and natural selection on codon usage bias in species of five Japanagallia species, the neutrality plot can be drawn with GC12 as the ordinate and GC3 (average of GC1 and GC2) as the abscissa. In the plot, as the slope gets closer to or equal to one, it indicates that mutational pressure plays a decisive role in codon usage. Whereas when the slope is far from one, it means that codon use is mainly influenced by other forces, such as natural selection (Sueoka, 1988; Sueoka, 1999).

Parity rule two plot

The parity rule two (PR2) plot is drawn with AT-bias value [A3/(A3 + T3)] as the y-axis and GC-bias value [G3/(G3 + C3)] as the x-axis to analyze the role of natural selection and mutational pressure on codon selection (Sueoka, 1995). The central point 0.5 of the graph shows that the influence of natural selection and mutational pressure are equal (Yang et al., 2015).

Effective number of codons

The effective number of codons (ENC) is often widely used to evaluate whether synonymous codons are used equally, with values ranging from 20 to 61 (Wright, 1990). A lower ENC value means that higher codon usage bias, and vice versa. Typically, an ENC value 35 reflects a strong preference for codons (Wright, 1990; Jiang et al., 2008). The ENC value 20 indicates that each amino acid is encoded with only one codon, while the ENC value 61 indicates that each codon can be used equally when encoding amino acids. The ENC values were estimated in MEGA6.06 (Tamura et al., 2013).

Correspondence analysis

The correspondence analysis (COA) is a multivariate statistical technique used to determine the major trends of codon usage variation among genes (Clarke & Greenacre, 1985; Greenacre, 1894; James & McCulloch, 1990). To understand trends in codon usage in axis1 and aix2 (Perrière & Thioulouse, 2002; Shields & Sharp, 1987), we performed COA analysis based on the RSCU values for 13 protein-coding genes in the mitogenomic using Pest software.

Grand average of hydropathy

The grand average of hydropathy (GRAVY) score was calculated by the sum of the products of the frequency of each amino acid and the corresponding hydropathy index of each amino acid (Kyte & Doolittle, 1982). The positive GRAVY value indicates that the protein is essentially hydrophobic, while the negative value represents that the protein is hydrophilic in nature. GRAVY is calculated using Galaxy (Afgan et al., 2018).

Phylogenetic analysis

The phylogenetic analysis was performed using whole mitogenomes of the 81 leafhopper species and five treehopper species as the ingroup, Cervaphis quercus (NC_024926) and Cacopsylla coccinea (NC_027087) were selected as members of the outgroups. Taxonomic information and mitochondrial accession numbers for each species are listed in Table S2. Three datasets were concatenated for phylogenetic analysis: (1) amino acid sequences of 13 PCGs with 3,594 amino acids (AA); (2) 13 PCGs and 2 rRNA genes (12S and 16S rRNA) with 12,281 nucleotides (PCG-rRNA); (3) the first and second codon positions of the 13 PCGs and 2 rRNA genes (12S and 16S rRNA) with 8,686 nucleotides (PCG12-rRNA). Each of the PCG (without stop codons) was aligned individually based on the invertebrate mitochondrial genetic code using the MAFFT algorithm in the Translator X online server (http://translatorx.co.uk/) (Abascal, Zardoya & Telford, 2010). Two rRNA genes were aligned using MAFFT v7 (https://mafft.cbrc.jp/alignment/server/index.html) with the Q-INS-I strategy (Katoh, Rozewicki & Yamada, 2019), and the poorly aligned positions and divergent regions were removed using the Gblocks 0.91b component in PhyloSuite v1.2.2 (Talavera & Castresana, 2007; Zhang et al., 2020). Next, the clipped sequences were concatenated into three different datasets using MEGA 6.06 (Tamura et al., 2013), the position of each gene fragment was recorded, and the resulting 15 alignments were manually checked and corrected.

The phylogenetic analyses were reconstructed using Bayesian Inference (BI) and Maximum Likelihood (ML) methods. The three concatenated datasets obtained were converted into “nex” and “phy” formats using Mesquite v2.75. The optimal partition scheme for each dataset and the best model for each partition was determined in PartitionFinder version 2.1.1 (Miller, Pfeiffer & Schwartz, 2010) using the Akaike information criterion correct (AICc), and using the “greedy” algorithm with branch lengths estimated an “linked” a greedy search algorithm (Lanfear et al., 2017). Maximum likelihood (ML) phylogenetic trees were constructed with IQ-TREE v1.6.3 (Nguyen et al., 2015) using an ultrafast bootstrap (UFB) approximation approach for 10,000 replicates (Nguyen et al., 2015). Bayesian analysis used defaulted settings by simulating four independent runs for 100 million generations and sampling once every 1,000 generations in MrBayes 3.2.6 (Huelsenbeck & Ronquist, 2001; Nylander et al., 2004), after the average standard deviation of split frequencies fell below 0.01, the initial 25% of samples were discarded as burn-in, and the remaining trees were used to generate a consensus tree and calculate the posterior probability (PP) of each branch. The beautification of phylogenetic trees using FigTree v1.4.4 and Adobe Illustrator CS6 software.

Mitogenome organization and composition analysis

The size of the complete mitochondrial genomes of the five Japanagallia species were 15,365 bp in J. curvipenis, 15,575 bp in J. malaise, 15,533 bp in J. multispina, 15,717 bp in J. turriformis, and 15,396 bp in Japanagallia sp., respectively. Figure 1 shows the circular maps of five Japanagallia species mitogenomes. Like the mitochondrial genome of most leafhoppers, each of the newly sequenced mitogenomes contained 37 coding genes (13 PCGs, 2 rRNAs genes, 22 tRNAs genes), and a control region.

At the same time, in order to accurately understand the codon bias, we analyzed the 13 protein-coding genes of each genome of five Japanagallia species. The nucleotide base frequencies (A, T, G, C and GC) as well as with frequencies at its third codon position (A3, T3, G3 and C3) were calculated (Table S3). Figure 2 depicted the average nucleotide compositions of all the genes. We found that the average percentage of nucleotide T was the most abundant base followed by nucleotide A, C and G, the trend of the nucleotide frequencies was T>A>C>G. Likewise, the average percentage of nucleotide composition at the third codon position revealed that the average percentage of T3 was the highest followed by A3, C3 and G3 (T3>A3>C3>G3).

Moreover, the average percentage of GC contents at different codon positions were also analyzed (Fig. 3), and the results showed that the percentage of GC2 content was higher than the content GC1 and GC3, and the content GC3 was the lowest in the five Japanagallia species, the overall percentage of GC content were less than 50%. Therefore, this compositional study suggested that the A, T, G, and C nucleotides are unequally distributed, and the genes mostly preferred T or A nucleotide, while the C or G were less preferred.

Relative synonymous codon usage analysis

We analyzed the RSCU value of each codon corresponding to amino acid in the CDS of 13 mitochondrial genes from the five Japanagallia species. By analyzing the overall RSCU (Table 1) of the CDS of mitochondrial 13 protein-coding genes, we observed that T ending codons were mostly preferred followed by A ending ones and thus the frequency of A or T was higher than that of G and C ending codons (Table S4). At the same time, we performed hierarchical clustering of the heatmap (Fig. 4). The color and degree of intensity clearly indicated the RSCU value and it varies from blue (low value of RSCU) to red (high value of RSCU). Orange and red color indicate over-represented codons (RSCU > 1.6), while blue under-represented codons (RSCU < 0.6). Figure 4 and Table S4 show that the codons UUA encoding Leucine amino acid (RSCU > 3) was favored by nature in all the 13 mitochondrial genes of five Japanagallia species.

Neutrality plot analysis

In our analysis, the neutrality plot results showed that the regression coefficients of GC12 versus GC3 in all 13 mitochondrial genes were less than 0.5 (ATP6: 0.076, ATP8: 0.017, COX1: 0.073, COX2: 0.211, COX3: −0.063, CYTB: 0.114, ND1: 0.235, ND2: 0.181, ND3: 0.032, ND4: 0.025, ND4L: 0.052, ND5: 0.040, ND6: 0.231) (Fig. 5), indicating that mutation pressure effect accounted only 1.69%–23.46%. Therefore, the above results suggested that natural selection might have played a dominant role, whereas mutation pressure and other factors might have played a lesser role in shaping CUB (Uddin & Chakraborty, 2020).

Parity plot analysis

The overall GC bias [G3/(G3 + C3)] and AT bias [A3/(A3 + T3)] was 0.441 and 0.514, separately. We plotted [G3/(G3 + C3)] along the x-axis and [A3/(A3 + T3)] along the y-axis of the graphical plot for the 13 mitochondrial genes of five Japanagallia species (Fig. 6). According to the parity plot, the central coordinate (0.5, 0.5) indicates that mutation pressure is the only force of codon bias, that is, the frequencies of nucleotides A and T are equal to the frequencies of C and G at the third position (Kawabe & Miyashita, 2003). In contrast, the degree of deviation from the center that is both mutations and selection will be involved in codon usage bias, with G and C unequal in the use of A and T bases (Sueoka, 1995). In our present study, overall the values of GC bias < 0.5 and AT bias > 0.5 (Table 2) can predict the preference of pyrimidine (U/C) over purine (A/C) at the third codon position (Zhang et al., 2018). Furthermore, we observed an uneven distribution of genes in the four regions of the PR2 plane, confirming the role of mutational pressure and natural selection in the use of mitochondrial gene codons in five Japanagallia species.

Meanwhile, our results showed that AT and GC are not used equally frequently at the third codon in 13 mitochondrial genes, suggesting that mutational pressure and natural selection may affect the codon preference in these genes.

Analysis the relationship of ENC and compositional attributes

The average ENC value for mitochondrial genes of the Japanagallia in our analysis was 34.564 (ATP6), 39.0844 (ATP8), 38.3864 (COX1), 36.4992 (COX2), 37.3108 (COX3), 36.1256 (CYTB), 31.1794 (ND1), 36.1856 (ND2), 34.291 (ND3), 32.482 (ND4), 31.2502 (ND4L), 33.8502 (ND5) and 32.1766 (ND6), which represented relatively high codon bias in these genes. At the same time, in order to better understand the impact of translational selection or mutation pressure on the codon usage of mitochondrial genes of the Japanagallia, we performed correlation analysis between ENC, the overall composition of nucleotide and its third position of the codon (Table 3). We found that in some mitochondrial genes, there was a significant positive correlation between homogenous nucleotides, while there was a significant negative correlation between most heterogeneous nucleotides, suggesting that the base composition bias of Japanagallia was affected by mutation pressure (Wong et al., 2010). In addition, there was a significant positive correlation between the ENC of COX2, ND3 and ND4 genes and GC3 (P < 0.05), indicating that the codon preference of these genes may be affected by GC content.

Correspondence analysis of codon bias

In order to investigate the trends of codon usage variation in five Japanagallia species for the 13 mitochondrial protein-coding genes, we performed COA using the RSCU value of codons. The distribution of axes in 13 mitochondrial genes is shown in Fig. 7. In our analysis, all zero row codons and two stop codons (UAA, UAG) in the matrix were deleted. We observed that the first principal axis accounted for 33.63% (ATP6), 47.67% (ATP8), 36.49% (COX1), 32.93% (COX2), 35.40% (COX3), 39.00% (CYTB), 41.17% (ND1), 33.9% (ND2), 36.09% (ND3), 36.56% (ND4), 43.15% (ND4L), 39.16% (ND5), 38.02% (ND6) of all variations in the gene set. The second axis accounted for 28.99% (ATP6), 40.47% (ATP8), 28.15% (COX1), 27.86% (COX2), 24.81% (COX3), 24.28% (CYTB), 24.10% (ND1), 30.60% (ND2), 27.66% (ND3), 29.96% (ND4), 33.69% (ND4L), 27.82% (ND5), 33.62% (ND6) of all variations within the gene get.

In addition, by observation we found that codons were mostly located near the axes and concentrated in the center of the graph, indicating that codon bias of these genes may be related to the base composition of mutation bias supporting the view of Butt, Nasrullah & Tong (2014). However, it was also observed that some gene codons were discretely distributed, indicating that the codon usage of 13 mitochondrial genes might be affected by natural selection and other factors (Wei et al., 2014).

Grand average of hydropathy analysis

We calculated the protein properties and amino acid composition of 13 protein-coding genes in the mitochondrial genome of five Japanagallia species. The grand average of hydrophilicity (GRAVY) of the 13 mitochondrial protein-coding genes (Table 4) was calculated to be positive, suggesting that the proteins of the Japanagallia are hydrophobic and may be to maintain their biological functions. At the same time, the overall frequency of the usage amino acid in the mitochondrial genes (Fig. 8) showed that the contents of Leucine (Leu), Isoleucine (Ile), Methionine (Met), Serine (Ser) and Phenylalanine (Phe) were higher than those of other amino acids, and Arginine (Arg) had the lowest content.

Phylogenetic analysis

The phylogenetic relationships were analyzed using ML and BI to reconstruct among the 81 species of the 12 subfamilies of leafhoppers, five species of treehoppers, and two outgroup species based on the following three datasets (AA, PCGRNA, PCG12RNA), and six phylogenetic trees were obtained (ML-AA, BI-AA, ML-PCGRNA, BI-PCGRNA, ML-PCG12RNA and BI-PCG12RNA) (Figs. 9 and 10; Figs. S1–S5).

In the six trees, most nodes received high nodal support values, and a few nodes only obtain medium or low support, and the monophyly of each subfamily was generally well supported, while the relationship among subfamilies were discrepant across analyses. At the same time, all phylogenetic trees showed the treehoppers clustered into one group and a sister group of Megophthalminae with high support, which supported the previous view that leafhoppers originated from the Cicadellinae, and further confirm that Cicadellidae is a paraphyletic group. These results are consistent with many previous studies based on mitogenomes (Dietrich et al., 2001, Dietrich et al., 2017; Zhao & Liang, 2016; Du, Dai & Dietrich, 2017; Duet al., 2017; Du, Dietrich & Dai, 2019; Skinner et al., 2020; Jiang et al., 2021; Wang, Wang & Dai, 2021). Within Megophthalminae, all phylogenetic trees showed that Japanagallia and Durgade were sister groups, which was consistent with previous reports based on molecular data (Wang et al., 2017).