Comparative and phylogenetic analysis of the mitochondrial genomes in basal hymenopterans

The Symphyta is traditionally accepted as a paraphyletic group located in a basal position of the order Hymenoptera. Herein, we conducted a comparative analysis of the mitochondrial genomes in the Symphyta by describing two newly sequenced ones, from Trichiosoma anthracinum, representing the first mitochondrial genome in family Cimbicidae, and Asiemphytus rufocephalus, from family Tenthredinidae. The sequenced lengths of these two mitochondrial genomes were 15,392 and 14,864 bp, respectively. Within the sequenced region, trnC and trnY were rearranged to the upstream of trnI-nad2 in T. anthracinum, while in A. rufocephalus all sequenced genes were arranged in the putative insect ancestral gene arrangement. Rearrangement of the tRNA genes is common in the Symphyta. The rearranged genes are mainly from trnL1 and two tRNA clusters of trnI-trnQ-trnM and trnW-trnC-trnY. The mitochondrial genomes of Symphyta show a biased usage of A and T rather than G and C. Protein-coding genes in Symphyta species show a lower evolutionary rate than those of Apocrita. The Ka/Ks ratios were all less than 1, indicating purifying selection of Symphyta species. Phylogenetic analyses supported the paraphyly and basal position of Symphyta in Hymenoptera. The well-supported phylogenetic relationship in the study is Tenthredinoidea + (Cephoidea + (Orussoidea + Apocrita)).


Results and Discussion
General features of the two newly sequenced mitochondrial genomes. We sequenced two nearly complete mitochondrial genomes from T. anthracinum (GenBank accession KT921411) and A. rufocephalus (GenBank accession KR703582). The sequenced region of the T. anthracinum mitochondrial genome was 15,392 bp long, with 13 protein-coding, two rRNA and 20 tRNA (except for trnQ and trnM) genes, and a partial A+ T-rich region. There were 619 bp intergenic nucleotides in total, in 18 locations, and the length of the intergenic spacers was 1-414 bp. The longest intergenic spacer (414 bp) was located at the start of the mitochondrial genome before trnY. Four pairs of genes overlapped each other, with a length ranging from 2 to 7 bp. Fourteen pairs of genes were directly adjacent.
The sequenced mitochondrial genome of A. rufocephalus was 14,864 bp long, with 13 protein-coding, two rRNA and 19 tRNA genes, and a partial A+ T-rich region. Part of the A+ T-rich region and three tRNA genes (trnI, trnQ and trnM) failed to be sequenced. There were 150 bp intergenic nucleotides in total, in 16 locations, and the length of the intergenic spacers was 1-27 bp. The longest intergenic spacer (27 bp) was located between cox3 and trnG. Five pairs of genes overlapped each other, with a length ranging from 1 to 8 bp. Thirteen pairs of genes were directly adjacent.
Gene rearrangement in the Symphyta. Compared with the putative ancestral gene arrangement of insects, there have been at least two rearrangement events in the mitochondrial genome of T. anthracinum, corresponding to the remote inversion of tnnY and the translocation of trnC from a location between trnW and cox1 to upstream of trnI. There was no gene rearrangement in the sequenced region of A. rufocephalus (Fig. 1). However, we could not exclude the possibility of gene rearrangement in the unsequenced tRNA cluster trnI-trnQ-trnM.
In the sequenced mitochondrial genomes of 10 Symphyta species (first 10 species in Table 1), gene rearrangement was found in seven species. No rearranged genes were found in the sequenced region of A. rufocephalus, Monocellicampa pruni, and Tenthredo tienmushana, in which the tRNA cluster trnI-trnQ-trnM is not Orussus occidentalis Figure 1. Mitochondrial genome organization of Symphyta referenced with the ancestral insect mitochondrial genomes. Genes are transcribed from left to right except for those indicated by a line. PCGs, rRNA, and A+ T-rich region genes are marked in yellow, pink, and white, respectively, while tRNA genes are marked in green and designated by the single-letter amino acid code.
detected. There was no PCG rearrangement in the Symphyta mitochondrial genomes. The rearranged tRNA genes were mainly from the tRNA clusters of trnI-trnQ-trnM and trnW-trnC-trnY, and trnL1 in Perga condei. Compared with the mitochondrial genomes of Apocrita, gene rearrangement is relatively conserved in the suborder Symphyta 3,14 ; however, it is frequent compared with most other insect orders 37-40 . Base composition. Three parameters, AT-skew, GC-skew, and A+ T content, are usually used in the investigation of the nucleotide-compositional behavior of mitochondrial genomes 29,41 . The mitochondrial genomes of T. anthracinum and A. rufocephalus both showed a very strong bias in nucleotide composition (A+ T% > G+ C%), which is typical in insect mitochondrial genomes. The A+ T content of the PCGs of all 10 Symphyta species was between 74.20% (Orussus occidentalis) and 80.30% (A. rufocephalus), while the A+ T content of T. anthracinum and A. rufocephalus was 79.30% and 80.30%, respectively ( Table 2). In the PCG sequences of Symphyta, all of the AT-skews were negative while most GC-skews were negative (except in A. rufocephalus, P. condei and Cephus pygmeus), indicating that the PCGs contained a higher percentage of T and C than A and G nucleotides (Table 3), as reported in most other insects 29 . The results of the comparative analysis of the A+ T content and the AT-and GC-skew of the PCGs within the sequenced Symphyta mitochondrial genomes are illustrated in a three-dimensional scatter-plot chart (Fig. 2). Species rom Symphyta and Apocrita were separated, except for V. bicolor, which was mixed with Symphyta. Within the Symphyta, species from each of the three superfamilies were clustered, except for T. anthracinum from Tenthredinoidea.  Table 2. Nucleotide compositions in regions of the mitochondrial genomes used in this study.
PCGs were located on the minority strand (N-strand) (Fig. 1). The entire length of the PCGs of T. anthracinum was 11,185 bp, while that of A. rufocephalus was 11,257 bp. The overall A+ T content of the 13 PCGs was 79.30% in the T. anthracinum mitochondrial genome, ranging from 72.58% (cox1) to 88.27% (atp8). In the A. rufocephalus mitochondrial genome, the A+ T content of the 13 PCGs was 80.30%, ranging from 73.45% (cox1) to 91.36% (atp8). In the T. anthracinum and A. rufocephalus mitochondrial genomes, all of the PCGs started with ATN codons. However, there are abnormal start codons in the Symphyta, including GTG (Cephus pygmeus, atp8), TTG (Cephus sareptanus, nad1), and AAA (C. pygmeus and C. sareptanus, nad2). In the T. anthracinum mitochondrial genome, two, four and seven PCGs started with ATA, ATG and ATT, respectively, while in the A. rufocephalus mitochondrial genome, four, three and six PCGs started with ATA, ATG and ATT, respectively.     The codon usage in the mitochondrial genomes of Symphyta also shows a significant bias towards A/T (Fig. 3). In the mitochondrial genomes of T. anthracinum and A. rufocephalus, Leu, Ile, Phe and Met were the most frequently used amino acids, while TTA (Leu), ATT (Ile), TTT (Phe) and ATA (Met) were the most frequent codons, as in other hymenopteran mitochondrial genomes 14,42,43 . All the three frequently used codons consisted of A and T, which may lead to the high A+ T content in the mitochondrial genome. It is obvious that the preferred codon usage is A/T in the third position rather than G/C, as almost all of the frequently used codons ended with A/T. In the mitochondrial genome of T. anthracinum, the codon Thr (ACG) was missing, while Val (GTG), Ala (GCC), Cys (TGC), Arg (CGC, CGG) and Ser (AGC) were missing in the mitochondrial genome of A. rufocephalus. The missing codons all have a high content of G/C in the third codon position.
Evolutionary rates of protein-coding genes. The rate of nonsynonymous substitutions (Ka) and synonymous substitutions (Ks), and the ratio of Ka/Ks, were calculated for PCGs of each Symphyta mitochondrial genome using A. appendiculatus as a reference sequence (Fig. 4). All of the Ka and Ks values were less than 1. Simultaneously, the ratios of Ka/Ks were all less than 1, indicating the existence of purifying selection in these species. Species of Symphyta showed lower evolutionary rates than Apocrita species, as is widely accepted 34 . Transfer RNA and ribosomal RNA genes. The anticodons of all the tRNAs of the two newly sequenced mitochondrial genomes were identical to other Symphyta species. In the mitochondrial genome of T. anthracinum, the tRNA length obtained was 1359 bp, with an A+ T content of 83.22%. The length of the tRNAs ranged from 64 bp (trnT) to 72 bp (trnH). In the 20 tRNA genes, 14 genes were located on the J-strand, while the others were encoded on the N-strand. In the mitochondrial genome of A. rufocephalus, the whole length of tRNAs obtained was 1282 bp, with a A+ T content of 83.37%. The length of the tRNAs ranged from 60 bp (trnV) to 72 bp (trnK). In the 19 tRNA genes, 12 genes were located on the J-strand, while the others were encoded on the N-strand.
In the two mitochondrial genomes, all of the tRNAs could be folded into classic clover leaf structures, except for trnV and trnS1 (AGN) in A. rufocephalus. The dihydrouridine (DHU) arm of these two tRNAs was a large loop instead of the conserved stem-and-loop structure, which is a typical feature of metazoan mitochondrial genomes 44 . All the secondary structures of the tRNA genes could be predicted by Mitos WebServer 45 . In the tRNAs of the two mitochondrial genomes, the amino acid acceptor (AA) stem and the anticodon (AC) loop were conserved as 7 bp, except for the AC-loops of trnL2 (CUN) and trnR in T. anthracinum and trnL2 (CUN) in A. rufocephalus, which were 9 bp. The length of tRNA usually depends on the size of the variable loop and the D-loo 46 . The DHU arm was 3-4 bp and the AC arm was 4-5 bp, while the TΨ C arm varied from 3-5 bp. The variable loops were less consistent, ranging from 4-10 bp. In the mitochondrial genomes, some base pairs were not the classic A-U and C-G, based on the secondary structure. Five mismatched base pairs were found in the tRNAs of each of T. anthracinum and A. rufocephalus. Among the five mismatched base pairs of T. anthracinum, four of them were U-U pairs located in the AA stem, while the other was an A-C pair located in the TΨ C stem (Table 4). However, in A. rufocephalus, there were three U-U pairs located in the AA stem and one U-U pair in the AC stem, and the other was an A-A pair located in the TΨ C stem (Table 4).
In the mitochondrial genome of T. anthracinum, the rrnL was 1351 bp long with an A+ T content of 84.46%, while the rrnS was 800 bp long with an A+ T content of 83.50%. In the mitochondrial genome of A. rufocephalus, the rrnL was 1339 bp long with an A+ T content of 85.06%, while the rrnS was 803 bp long with an A+ T content of 83.44%. In the 10 Symphyta species, the length of rrnS ranged from 728 (P. condei) to 1014 bp (C. sareptanus), while the length of rrnL ranged from 1336 (O. occidentalis) to 1386 bp (C. cinctus).

Intergenic spacers and overlapping regions. We identified 20 intergenic spacers in the T. anthracinum
mitochondrial genome, and 17 intergenic spacers in the A. rufocephalus mitochondrial genome. In the T. anthracinum mitochondrial genome, the longest intergenic spacer was located at the start of the mitochondrial genome before trnY and the second longest region was the incomplete A+ T-rich region with a length of 100 bp, while the other 18 regions ranged from 1 to 47 bp. In the mitochondrial genome of A. rufocephalus, the longest intergenic spacer was also the incomplete A+ T-rich region, with a length of 55 bp, and the other 16 regions ranged from 1 to 27 bp. In both species, we did not sequence the full length of the A+ T-rich region. The content of these two incomplete A+ T-rich regions was 78% and 85.45%, respectively.
There was an "ATTATAA" motif between nad4 and nad4l in the N-strand of the T. anthracinum mitochondrial genome, but this intergenic region did not exist in the A. rufocephalus mitochondrial genome, because of the direct adjacency of nad4 and nad4l. There was an overlap sequence of "ATGATAA" present between atp8 and atp6. This characteristic 7 bp overlap of the two PCGs located on the J-strand appears in most species of Symphyta and other insects.
Phylogenetic relationships. We reconstructed the phylogeny among the Symphyta (Fig. 5). Before reconstructing the phylogenetic tree, the saturation of different partitions generated from PartitionFinder was tested (Fig. 6). The results indicated that partitions 6, 9, 10 and 11 identified by PartitionFinder had experienced substitution saturation. From the results of the data partitioning, we found that these four partitions contained all the third positions of PCGs and two rRNA genes. We also tested the saturation of all partitions and different codon  positions in the PCGs (Fig. 6). The results showed that substitution saturation occurred in the third position. Hence, we reconstructed the phylogenetic relationship using four data matrices, i.e. 11,931 sites in the P123 matrix (containing the three codon positions of PCGs), 16,548 sites in the P123R matrix (containing the three codon positions of PCGs, two rRNA genes and 20 tRNA genes), 9483 sites in the P12T matrix (excluding partitions 6, 9, 10 and 11) and 3977 sites in the AA matrix (containing the amino acid sequences). For Bayesian and ML analyses, all matrices of nucleotide sequences (P123, P123R and P12T) generated the same topology with high posterior probabilities. However, the matrix of the amino acid sequences generated from both BI and ML analyses showed a different topology to that of the nucleotide sequences, in which Cephoidea forms a sister group with Orussoidea, and then these two groups form a sister group to Apocrita. These topologies have not been reported in previous studies 9,20,21 . However, the relationship among the Tenthredinoidea was stable no matter which matrix and method was used. We present the BI and ML results based on three matrices of the nucleotide sequences in Fig. 6.
Our analyses support the paraphyly and basal position of Symphyta in Hymenoptera and are congruent with traditional views, while Orussoidea forms a sister group with the suborder Apocrita 9,16,34 . There are three strongly supported lineages in the Symphyta, Tenthredinoidea, Cephoidea, and Orussoidea. Our analyses support a phylogenetic relationship of Tenthredinoidea + (Cephoidea + (Orussoidea + Apocrita)) in the Symphyta, as in most other analyses 9,20,21 . The analyses support the monophyly of Tenthredinoidea. Within the Tenthredinoidea, Tenthredinidae forms a sister relationship with Cimibicidae, while the family Pergidae forms a sister lineage to Cimibicidae + Tenthredinidae. Within the Tenthredinidae, the subfamilies Tenthredininae and Allantinae form a sister lineage and are then a sister to the Nematinae family, as revealed in Malm and Nyman 9 , and this is also well defined by morphology 20,21 .

Materials and Methods
Sample preparation and DNA extraction. A specimen of T. anthracinum was collected from Hanmi in Tibet, China, while a specimen of A. rufocephalus was collected from Songshan Forest Park in Yanqing County, Beijing, China. All the specimens were stored at − 80 °C in 100% ethanol prior to DNA extraction. Total genomic DNA was extracted from legs of the specimens using the DNeasy tissue kit (Qiagen Hilden, Germany), following the manufacturer's protocols. The voucher specimens are kept in the Evolutionary Biology Laboratory of Zhejiang University, China.
Primer design, PCR amplification and sequencing. First, we used a set of universal primers for the insect mitochondrial genome 2,8 for amplification and sequencing of partial gene segments. Then we designed specific primers based on the sequenced segments to amplify regions to bridge the gaps between different segments. PCR was carried out using Takara LA Taq (Takara Biomedical, Japan) with the following conditions: initial denaturation at 96 °C for 3 min and then 40 cycles at 95 °C for 30 s, annealing at 45-53 °C for 30 s, extension at 60 °C for 1-2 min, and final extension at 60 °C for 10 min. PCR components were added following the Takara LA Taq protocols. The PCR products were directly sequenced by TSINGKE Company (Beijing, China) using a primer-walking strategy from both strands.
Mitochondrial genome annotation. The tRNA genes were initially identified using the tRNAscan-SE search server, with a Cove cutoff score of 5, and Mitos WebServer. The tRNAscan-SE search server set the parameters that the source was Mito/Chloroplast and the genetic code was the Invertebrate Mito genetic code, while Mitos WebServer set the parameters that the genetic code was the Invertebrate Mito genetic code. The tRNAs that could not be found using these two approaches were confirmed by sequence alignment with their homologs in related species. Protein-coding genes were identified by Blast searches in GenBank, according to other published hymenopteran mitochondrial genomes. The rRNA genes and control region were identified by the boundary  Figure 6. Phylogenetic relationships of the Symphyta inferred from nucleotide sequences of the mitochondrial genome using three matrices (P123, P123R and P12T), referring to the Bayesian/Maximum Likelihood methods. The numbers separated by "/" indicate the posterior probability and bootstrap values of the corresponding nodes (P123 using BI/P123R using BI/P12T using BI/P123 using ML/P123R using ML/P12T using ML). "*" indicates that the node was fully supported by all six inferences.
Scientific RepoRts | 6:20972 | DOI: 10.1038/srep20972 of the tRNA genes and by comparing them with other insect mitochondrial genomes. Secondary structures of tRNAs were predicted using the tRNAscan-SE search server (Lowe and Eddy, 1997) and Mitos WebServer 45 , and others were predicted manually if they could not be found using the tRNAscan-SE search server.
Comparative analysis of the mitochondrial genomes from Symphyta. We compared the mitochondrial genomes of 10 species from the Symphyta, including the two newly sequenced ones in our study.
Features of gene arrangement, base composition, codon usage and base substitution of PCGs were analyzed. Since there were several tRNA genes not obtained in some species, we analyzed the base composition using the PCGs only.
The base composition was calculated by Mega6 47 . The AT and GC asymmetries, called the AT-skew and GC-skew, were calculated according to the formulas suggested by Hassanin et al. 41 : AT-skew = (A − T%)/ (A + T%) and GC-skew = (G − C%)/(G + C%). The intergenic spacers and overlapping regions between genes were counted manually. The relative synonymous codon usage (RSCU) of all protein-coding genes was calculated by CodonW (written by John Peden, University of Nottingham, UK). We used the software package DnaSP 4.0 48 to calculate the number of synonymous substitutions per synonymous site (Ks) and the number of nonsynonymous substitutions per nonsynonymous site (Ka) for each Symphyta mitochondrial genome, using that of Ascaloptynx appendiculatus from Neuroptera as a reference sequence.

Phylogenetic analysis. Taxa selection.
To investigate the phylogenetic relationships within the Symphyta, we used all eight species of Symphyta for which the mitochondrial genomes had previously been published, and the two newly sequenced mitochondrial genomes in this study ( Table 1). The 10 species can be divided into three superfamilies, Tenthredinoidea, Cephoidea and Orussoidea. There were six species belonging to Tenthredinoidea, three species belonging to Cephoidea, and one species belonging to Orussoidea. Four mitochondrial genomes from the Apocrita were used to construct the phylogenetic tree, because of the close relationship between Apocrita and Symphyta. The four insects represented four families, four superfamilies and both parts of Apocrita. Vanhornia eucnemidarum (Vanhorniidae) 49 and Diadegma semiclausum (Ichneumonidae) 43 represented Terebrantia, while Apis mellifera (Apidae) 50 and Vespa bicolor (Vespidae) 51 represented Aculeata (Table 1). A. appendiculatus 52 from the order Neuroptera was set as an outgroup because of the close relationship between Hymenoptera and Neuroptera.
Sequence alignment, data partition and substitution model selection. We used MAFFT version 7.205 for the alignment of protein-coding and RNA genes, which implements consistency-based algorithms 53 . We used the G-INS-I and Q-INS-I algorithms in MAFFT 54 for protein-coding and RNA alignment, respectively. The alignment of the nucleotide sequences was guided by the amino acid sequence alignment using the Perl script TranslatorX version 1.1 55 .
Data partitioning, and the ability to apply specific models to different partitions, is ideal for analyzing mitochondrial genomes 3 . We used PartitionFinder version 1.1.1 56 to simultaneously confirm partition schemes and choose substitution models for the matrix. The search models for DNA and amino acid sequences were set to be "mrbayes" and "all protein", respectively. The greedy algorithm was used, with branch lengths estimated, to search for the best-fit partitioning model. The branch lengths were set as linked.
Phylogenetic inference. We constructed the phylogenetic relationships among the Symphyta with the Bayesian inference method (BI) using Mrbayes v3.2.5 57 , and the Maximum Likelihood (ML) method using RAxML v8.0.0 58 , based on the nucleotide sequences of the 13 protein-coding genes, 20 tRNAs (except trnI and trnQ, as they were missing in more than half of the 14 species) and two rRNAs. In the BI, we used the GTR+ I+ G/ GTR+ G/HKY+ I+ G/HKY+ G model for nucleotide sequences and the MtArt+ I+ G+ F/MtArt+ I+ G model for amino acid sequences (Table 5). Four simultaneous Markov chains were run for 10 million generations, with tree sampling occurring every 1000 generations and a burn-in of 25% trees. In the ML analyses, we used the GTRGAMMA and PROTCATMTART models for nucleotide and amino acid partitions, respectively. For each analysis, 200 runs were conducted to find the highest-likelihood tree, followed by analysis of 1000 bootstrap replicates.