The Complete Nucleotide Sequence of the Mitochondrial Genome of Bactrocera minax (Diptera: Tephritidae)

The complete 16,043 bp mitochondrial genome (mitogenome) of Bactrocera minax (Diptera: Tephritidae) has been sequenced. The genome encodes 37 genes usually found in insect mitogenomes. The mitogenome information for B. minax was compared to the homologous sequences of Bactrocera oleae, Bactrocera tryoni, Bactrocera philippinensis, Bactrocera carambolae, Bactrocera papayae, Bactrocera dorsalis, Bactrocera correcta, Bactrocera cucurbitae and Ceratitis capitata. The analysis indicated the structure and organization are typical of, and similar to, the nine closely related species mentioned above, although it contains the lowest genome-wide A+T content (67.3%). Four short intergenic spacers with a high degree of conservation among the nine tephritid species mentioned above and B. minax were observed, which also have clear counterparts in the control regions (CRs). Correlation analysis among these ten tephritid species revealed close positive correlation between the A+T content of zero-fold degenerate sites (P0FD), the ratio of nucleotide substitution frequency at P0FD sites to all degenerate sites (zero-fold degenerate sites, two-fold degenerate sites and four-fold degenerate sites) and amino acid sequence distance (ASD) were found. Further, significant positive correlation was observed between the A+T content of four-fold degenerate sites (P4FD) and the ratio of nucleotide substitution frequency at P4FD sites to all degenerate sites; however, we found significant negative correlation between ASD and the A+T content of P4FD, and the ratio of nucleotide substitution frequency at P4FD sites to all degenerate sites. A higher nucleotide substitution frequency at non-synonymous sites compared to synonymous sites was observed in nad4, the first time that has been observed in an insect mitogenome. A poly(T) stretch at the 5′ end of the CR followed by a [TA(A)]n-like stretch was also found. In addition, a highly conserved G+A-rich sequence block was observed in front of the poly(T) stretch among the ten tephritid species and two tandem repeats were present in the CR.


Introduction
The family Tephritidae, generally known as ''true'' fruit flies, includes 471 genera and 4257 species distributed throughout the temperate and tropical areas of the world. Many species are of critical importance to man either as pests of fruit and vegetable crops or as beneficial species for the control of weeds [1]. The fruit fly Bactrocera minax Enderlein (Diptera: Tephritidae), generally known as the Chinese citrus fruit fly, has been a serious pest of commercial citrus crops in China for more than half a century [2]. This species has been recorded in southern China, India (West Bengal and Sikkim) and Bhutan [2,3] wild and cultivated citrus species [4]. Some hosts are endemic to southern China and the eastern Himalayan region [5] but B. minax has been reported on the kumquat Fortunella crassifolia [6] and the boxthorn Lycium chinense [2].
B. minax was first collected from India and Sikkim and designated B. minax Enderlein [1]. Drew [3] provided a detailed description and illustration of the B. minax type specimens collected in 1920 and assigned the species to the genus Bactrocera (Polistomimetes). White and Wang [7] designated a lectotype of B. minax and assigned the species to the Bactrocera (Tetradacus); in addition, they indicated that Bactrocera citri Chen, collected from China in 1940, should be placed in synonymy with B. minax.
A wide variety of questions about the biology and phylogeny of B. minax have been addressed with the aid of molecular tools. These studies could have used two main sources of genetic data; namely, nuclear sequence data and, most frequently, mitochondrial sequence data. Insect mitochondrial DNA (mtDNA) usually occurs as a double-stranded closed circular molecule, ranging in size from 14-20 kb and generally encoding 13 protein-coding genes (PCGs), two ribosomal RNAs (rRNAs) and 22 transfer RNA (tRNAs), which is conserved across bilaterian metazoans with only a few exceptions (e.g. loss of a small number of genes in some derived groups) [8]. The molecule contains at least one sequence of variable length known as the A+T-rich region or control region Table 1. Summary of primers used for complete mtgenome of B. minax amplification. Note: Lowercase ''r'' behind some primer names represents these primers were designed on basis of Simon et al. (1994); Lowercase ''r'' behind some primer names represents these primers were designed by us. doi:10.1371/journal.pone.0100558.t001 (CR), which contains initiation sites for transcription and replication [9] and ranges in size from tens to several thousands of base pairs [10][11][12][13]. As the results of highly conservative gene structures among phyla, maternal inheritance, high copy number and relatively fast evolution rates compared to nuclear DNA [14], mitochondrial genome (mitogenome) sequences have been regarded as useful molecular markers in studies focusing on comparative and evolutionary genomics, molecular evolution, phylogenetics, phylogeography and population genetics [15]. Many complete or nearly complete mitogenomes have been sequenced and comparative analyses at the genus or species level have used multiple complete mitochondrial genes instead of one or partial genes, including molecular systematics [16][17][18][19][20], population genetics/phylogeography [16], diagnostics [21], molecular evolutionary studies [13,22,23], the frequency and type of gene rearrangements [24,25] and the evolution of genome size [26]. To date, more than 500 insect mitogenomes have been sequenced from all orders, including 77 dipterans in 24 families, and are available in Genbank. In this study, we sequenced the complete sequence of the mitogenome of B. minax (Diptera: Tephritidae).
Genbank contains information for only ten Tephritidae species; Bactrocera oleae, Bactrocera tryoni, Bactrocera philippinensis, Bactrocera    (Table 1). Amplification was done in a thermocycler (Eppendorf Mastercycler 5333) in 50 ml reactions containing 5 ml of 25 mM MgCl 2 , 5 ml of 106PCR Buffer (Mg 2+ free), 8 ml of a dNTP mixture (2.5 mM each), 3 ml of 10 mM each primer, 0.5 ml of 5 U/ml Taq polymerase (Takara Biomedical, Japan) and 2 ml of a 1/10 dilution of the DNA extract. Amplification conditions were: 59 of pre-PCR denaturation at 94uC followed by 34 cycles of 30 s at 94uC, 1 min at 40-58uC (depending on the primer pair) and 2 min at 72uC. The F21 fragment ( Fig. 1) was amplified using LA Taq (Takara Biomedical, Japan) and a cycle consisting of a pre-PCR denaturation at 96uC for 2 min followed by 30 cycles of 10 s at 98uC and 2 min at 58uC with a final elongation step of 10 min at 72uC. PCR products were separated by electrophoresis and purified using a QIAquick Gel Extraction Kit (QIAGEN). PCR products were sequenced directly on both strands using amplification and additional ad hoc primers as needed. Individual sequences were combined in a consensus contig using DNAStar package software (DNAStar Inc.).

Sequence analysis and gene annotation
Genes encoded on the B. minax mitogenome were located initially by comparison to homologous full-length insect mitochondrial sequences using DNAStar. Nucleotide sequences of PCGs were translated using the invertebrate mtDNA genetic code. tRNA genes were identified initially using tRNAscan-SE Search Server version 1.21 (available online at http://lowelab.ucsc.edu/ tRNAscan-SE/) [32] and refined using tRNAscan-SE and RNAshapes [33]. The presence and secondary structures of tRNA genes that could not be located by tRNAscan-SE owing to variant morphology were annotated manually by comparison to the sequences of other insect tRNAs [34][35][36][37]. Codon usage analysis and relative synonymous codon usage (RSCU) in PCGs were calculated using CodonW version 1.4.2 (John Peden, available at http://codonw.sourceforge.net/index.html) [38]. Potential secondary structure folds of non-coding sequences and sequences in the CR were calculated with the DNA mfold web server using default settings (http://mfold.bioinfo.rpi.edu/cgi-bin/dna-form1. cgi) [39]. The presence of tandem repeats in the CR was investigated using the Tandem Repeats Finder available online (http://tandem.bu.edu/trf/trf.html) [40]. The A+T content and nucleotide substitution frequency at synonymous sites and nonsynonymous sites (the number of synonymous substitutions per site and the number of non-synonymous substitutions per site) were calculated on the basis of the data using MEGA 4.0 [41]. The correlation analysis was done by the bivariate method using SPSS version 13 (SPSS Inc., Chicago, IL). The overall average amino acid distance among each of the PCGs from ten tephritid species (B. minax, B. oleae, B. tryoni B. dorsalis B. philippinensis, B. carambolae, B. papayae, B. correcta, B. cucurbitae and C. capitata) were calculated by the method of Poisson distances by MEGA 4.0 [41]. The complete B. minax mtDNA sequence was deposited in Genbank under accession no. HM776033.

Genome organization
The mitochondrial genome of B. minax is a closed circular molecule of 16043 bp; hence, it is longer than the other nine tephritid mitogenomes available (range 15,815 bp in B. oleae to 15,980 bp in C. capitata) but is still well within the range of other insect mitogenomes (14,503 bp in Rhopalomyia pomum [42] to 19517 bp in Drosophila melanogaster [11]). The gene content is typical of metazoan mitogenomes, with 13 PCGs (cox1-3, cob, nad1-6, nad4l, atp6 and atp8), 22 tRNAs and two genes for ribosomal RNA subunits (rrnS and rrnL). A long uninterrupted non-coding region of 1141 bp, likely homologous to the insect A+T-rich region, is present between rrnS and trnI, corresponding to position 14,903 to 16,043 in the annotated sequence. The gene order in the   B. minax mitogenome corresponds to the typical and plesiomorphic state hypothesized for the Pancrustacea, and is shared with all tephritids analyzed to date (Fig. 1). Genes in the B. minax mitogenome overlap by a total of 43 bp, distributed in 12 segments from 1 to 17 bp long and are separated by a total of 178 bp dispersed in 16 intergenic spacers from 2 to 42 bp (without taking the tRNA-like sequence into account; Table 2). Despite its relatively large size, the B. minax mitogenome has more overlapping sequences between genes compared to those of other tephritids; genes overlap by a total of 35
Considering the two strands separately, the PCGs on the Majority strand (J-strand, nine PCGs are located on this strand) (61.5%) have a lower A+T content compared to the Minority strand (N-strand, the other four PCGs are located on this strand) (68.9%). Furthermore, PCGs encoded on the J-strand have a comparable content of A (31.0%) and T (30.5%), whereas PCGs on the N-strand show a strong bias for T content (46.3%) compared to A content (22.6% A). The above situation has been observed in the other tephritid mitogenomes available (data not shown) and in other insects [34][35][36][37][43][44][45][46][47][48][49][50]. However, tRNAs on the two opposite strands have nearly equal A+T contents, which has been found in the other nine tephritid species. For three PCG codon positions, the third codon positions have significantly higher A+T content than the first and second codon positions owing to genetic code degeneracies. In particular, T in each codon position of PCGs on the N-strand is over-represented. With exception of the second codon position over-representing T, however, the first and third codon positions of PCGs show a preponderance of A on  the J-strand and T on the N-strand, which is similar to many insect mitogenomes [34][35][36][37][43][44][45][46][47][48][49][50] (Table 3).
The base compositional bias for A+T in PCGs is reflected in the relative synonymous codon usage statistics of the B. minax mitogenome (Table 4). With the exception of amino acid His, codons with A or T in the third codon position are generally strongly over-represented compared to codons terminating with either G or C. The ratio of G+C-rich (Pro, Ala, Arg and Gly) codons to A+T-rich codons (Phe, Ile, Met, Tyr, Asn and Lys) in B. minax PCGs was 0.44, which is higher compared to the other nine tephritids B. dorsalis With the exception of first codon positions, G is underrepresented compared to C in coding genes on the J-strand (PCGs, tRNAs, CR and intergenic nucleotides), while the G content is higher compared to C in coding genes on the N-strand (PCGs, tRNAs and rRNAs). This base compositional bias is in line with the general trend in the mitogenome toward a lower G content [51].
Base compositional heterogeneity and among-site rate variation (ASRV) are known to affect phylogenetic inference, resulting in the identification of incorrect phylogenetic relationships [52]. The easiest solution is simply to avoid non-stationary genes [53] but most earlier studies used relatively intuitive mitogenome data partitioning schemes, including by gene type (PCG, rRNA and tRNA), by gene, by codon position, by codon and gene, or by the strand on which the coding gene is located [15]. Inevitably, different intuitive partitioning schemes can each result in strong conflicting topologies, especially at deeper phylogenetic levels [25,54,55]. Therefore, selection of stationary, reversible compositional homogeneous is vital for reliable phylogenetic inference [52,56].
Many earlier studies were focused on the A+T content of different genes or regions to investigate the base compositional heterogeneity and among-site rate variation ASRV [57]. For mitogenomes, composition bias of A+T content was verified in most earlier studies; e.g. A+T content was usually over-represented in non-coding regions [58] and the third codon position generally had stronger A+T composition bias compared to the other two codon positions [59] etc.. We asked how variability between PCGs is related to underlying A+T content and its distribution across synonymous and non-synonymous sites.
In this study, the A+T content of zero-fold sites (P 0FD ), two-fold (P 2FD ) and four-fold degenerate sites (P 4FD ) was determined for each of the PCGs from ten tephritid species (B. minax, B. oleae, B. tryoni B. dorsalis B. philippinensis, B. carambolae, B. papayae, B. correcta, B. cucurbitae and C. capitata) (Fig. 2). Nucleotide substitution frequency was calculated in P 0FD , P 2FD and P 4FD for each of the PCGs among five tephritid species (Fig. 3). After analyzing the correlation between A+T content and nucleotide substitution frequency for each of the PCGs, we found a significant positive correlation between A+T content percentage of zero-fold degenerate sites (AT 0F ) and nucleotide substitution frequency at P 0FD (r = 0.735, P = 0.004) as well as between A+T content percentage of four-fold degenerate sites (AT 4F ) and nucleotide substitution frequency at P 4FD (r = 0.864, P = 0.000) ( Table 5). Correlation analysis indicated there is a significant positive correlation between AT 0F and ASD (r = 0.752, P = 0.003), ASD and the nucleotide Table 5. A+T content percentage and nucleotide substitution frequency at 0-fold degenerate sites (P 0FD ), 2-fold degenerate sites (P 2FD ) and 4-fold degenerate sites (P 4FD ) (the number of substitutions per P 0FD , P 2FD  substitution number of zero-fold degenerate sites/the nucleotide substitution number of all degenerate sites (R 0F/all ) (r = 0.983, P = 0.000), AT 0F and R 0F/all (r = 0.760, P = 0.003) ( Table 6). Interestingly, the significant positive correlation was observed between AT 4F and the nucleotide substitution number of four-fold degenerate sites/the nucleotide substitution number of all degenerate sites (R 4F/all ) (r = 0.809, P = 0.001); however, there was significant negative correlation between AT 4F and ASD (r = 2 0.828, P = 0.000), between R 4F/all and ASD (r = 20.970, P = 0.000) ( Table 6). On the basis of the above results, we can hypothesize divergence at the amino acid level of less well conserved PCGs is due to higher A+T at P 0FD in those genes and/ or lower A+T at P 4FD . On the basis of this result, when we choose which PCGs are used to analyze phylogenic relationships for different evolutionary time scales, the A+T content of P 0FD and/or P 4FD of PCGs could be useful to judge the homogenesis of PCGs. Nucleotide substitution is considered to be a reflection of evolution at the molecular level. Many earlier studies indicated the substitution was directional bias across different genes in the mitogenome [15]. Some researchers have proposed variation of A+T% among taxa is associated with directional mutation pressure and has a phylogenetic component [57,60,61]. In this study, with the exception of nad4, all PCGs had significantly lower variation of A+T content among the ten tephritid species at P 0FD compared to both P 2FD and P 4FD sites. We observed that, with the exception of nad4, P 0FD sites had lower nucleotide substitution frequency compared to both P 2FD and P 4FD sites (Fig. 3). The P 0FD of nad4 had a higher nucleotide substitution frequency (0.783) compared to both P 2FD (0.172) and P 4FD (0.004), and the R 0F/all was 0.936. As a result of functional constraints, the number of nucleotide substitution per non-synonymous site is usually lower than that per synonymous site [62]. In this study, a higher nucleotide substitution frequency at P 0FD of nad4 indicates the non-synonymous nucleotide substitution frequency was higher compared to the synonymous sites for this gene. Higher number of nucleotide substitution per non-synonymous site has been observed at the variable-region genes of immunoglobulins [63] and some genes of the histocompatibility complex [64] but this is the first reported occurrence in the mitogenome.

Protein-coding genes
With the exception of cox1 and nad3, all protein coding genes start with an ATN codon, with ATG used in cox2, atp6, cox3, nad4, nad4l, nad6 and cob, ATT in nad2, atp8 and nad5 and ATA in nad1. Genes for cox1 and nad3 used TCG and GTC as initiation codons, respectively. The initiation codon for cox1 was TCG(S) in B. minax, which was observed in other Diptera species [54]. GTC being the initiation codon for nad3 was a new observation in tephritids, but it is common in other insects [65].
With the exception of nad3, nad5 and nad1, all PCGs are terminated by complete stop codons: TAG is used for nad2, atp6 and cob, TAA is used for cox2, atp8, cox3, nad4, nad4l and nad6 and TA is used for cox1. The remaining genes, nad3, nad5 and nad1, are terminated by incomplete stop codons ''T''.
4. Transfer RNA genes, ribosomal RNA genes and tRNAlike structure All of 22 tRNA genes typical of metazoan mitogenomes were identified in the B. minax mitogenome, and the predicted structures are shown in Fig. 4. All tRNAs display a typical clover-leaf secondary structure, except for trnS (AGN) , where the DHU arm appears to be replaced by seven unpaired nucleotides, a feature typical of other animal mitochondria [66]. Nuclear magnetic resonance analysis of the tertiary structure of nematode trnS (AGN)) suggested such aberrant tRNA can fit the ribosome by adjusting its structural conformation and function in a way similar to that of usual tRNAs in the ribosome [67].
Like most insect tRNAs, all B. minax tRNAs have a length of 7 bp for the anticodon loop, 7 bp for the acceptor stem and 5 bp for anticodon stem. Most of the size variability in the B. minax tRNA genes originated from length variation in the DHU arms (loop size 4-9 bp, stem size 3-4 bp) and the TYC arms (loop size 2-9 bp, stem size 3-5 bp); in addition, trnA and trnH contained U-U mismatches. trnS (UCN) encodes an A-C mismatch, trnH encodes  an A-G mismatch and trnRtrnR has a U-C mismatch in the acceptor stem. Additionally, trnV contains a U-U mismatch in the TYC stem. Anticodon sequences were the same as in B. dorsalis, B. oleae, B. tryoni and C. capitata, which are considered common for other insects, including Gryllotalpa orientalis [68], Philaenus spumarius [35], Phthonandria atrilineata [50] and Artogeia melete [36].
On the basis of the sequence similarity of B. dorsalis, the two genes coding for the small and the large ribosomal subunits were located in the B. minax mitogenome between trnL (CUN) and trnV and between trnV and the CR region. The length of B. minax rrnS and rrnL was 782 bp and 1333 bp, respectively, similar to B. dorsalis, B. oleae and C. capitata.

Intergenic spacers
In B. minax, the two longest intergenic spacers were 42 bp between trnC and trnY and 28 bp between trnR and trnN. In B. dorsalis, the second longest intergenic spacer was 45 bp between trnC and trnY. In B. tryoni, the second longest intergenic spacer was 33 bp between trnR and trnN and the third longest intergenic spacer was 30 bp between trnC and trnY. In B. oleae, the longest intergenic spacer was 28 bp between trnR and trnN. In B. minax, however, only a 10 bp intergenic spacer was observed between trnQ and trnM, which is shorter compared to 66 bp in B. dorsalis, 71 bp in B. tryoni and 47 bp in C. capitata at the same location. Yu et al. [48] reported the 45 bp intergenic spacer located between trnC and trnY in B. dorsalis had a clear counterpart in the CR with the first 33 of 45 bp matching. These counterparts were predicted to form a small internal stem and a long stem structure pairing with the partially complementary sequence in the CR. A similar phenomenon was observed in the B. tryoni mitogenome, where both the second longest (33 bp between trnR and trnN) and the third longest intergenic spacer (30 bp between trnC and trnY) have clear counterparts (32 out of 33 bases and 25 out of 30 bases, respectively) on the N-strand of the CR. These two intergenic spacers have highly significant similarity and their counterparts were located in the same position of the CR. We asked whether the 42 bp intergenic spacer located between trnC and trnY in B. minax had these features. The first 15/42 bp of the spacer have a clear counterpart in the CR at positions 15,670-15,684. The 42 bp of intergenic spacer was predicted to form two stem-loop secondary structures with 4 bp loops and one with a 3 bp stem and the other with a 4 bp stem. The first 15 of the 42 bp formed one of the two structures; a 4 bp stem with a 4 bp loop and a 3 bp flanking sequence. The counterpart in the CR also formed a long stem structure with the neighboring sequence. Yu et al. [48] compared the 33 bp counterpart in the CR from B. dorsalis with the B. oleae CR and found 25 of the 33 bp were identical. Surprisingly, of the original 33 bases present in the B. minax CR, 23 were identical. Therefore, the results obtained in this study support the hypothesis that the secondary structures of the counterparts in both the intergenic spacer and the CR might have a major role in recombination [48,69].
Additionally, all four intergenic spacers have clear counterparts in the CR of the ten tephritid species (data not shown) but these intergenic spacers cannot form the secondary structures (even though some can be predicted to form stem-loop structures with 2-3 bp stems). Some earlier studies focused on longer intergenic spacers with potential secondary structure and tried to find original sequences and structures in the CR [48]. Even among the close tephritid species, however, these longer intergenic spacers had significantly different features, including sequence, length and location. Cameron et al. [70] suggested the possibility that stemloop structures instead of tRNAs in the 39 end of PCGs enhance the rearrangement. Two of four small intergenic spaces locate the 39 end of PCGs without forming stem-loop structures. These results might explain why no rearrangement was found in tephritid species. This is the first report of shorter intergenic spacers with highly conserved sequences and locations among four tephritid species, which should attract more attention to the shorter intergenic spacers, even though the functions of these are not clear.

CR
The CR has a high A+T content among the mitochondrial genes of both vertebrates and invertebrates, and the initiation of replication is one of the most interesting features of this region [8]. Zhang and Hewitt [71] proposed conserved structural features on the basis of comparison of the CRs of one dipteran and two orthopteran species. These features include: (1) a poly(T) stretch at the 59 end of the CR; (2) a [TA(A)] n -like stretch after the poly(T) stretch; (3) a highly conserved stem-loop structure; (4) a stem-loop structure with a highly conserved flanking sequence of a TATA consensus at the 59 end and a G(A) n T consensus at the 39 end; and (5) a G+A-rich sequence downstream of the secondary structure. The B. minax CR was found to have three of the five features proposed by Zhang and Hewitt [71].
The CR from four tephritid species, including B. minax, presented a conspicuous poly(T) stretch at the 59 end. This sequence stretch has been found to be conserved within hymenoptera [49]. Further, the poly(T) stretch has been observed to be followed by a [TA(A)] n -like stretch (Fig. 5). Our results suggest that this poly(T) region might be involved in the control of transcription and/or replication, or have some other unknown functions [10]. Additionally, a highly conserved G+A-rich sequence block was found in front of the poly(T) stretch among the four tephritid species and these sequences can be predicted to form secondary structures with a stem-loop. The highly conserved G+A-rich sequence with a poly(T) stretch nearby has been found in other dipteran and orthopteran species [71].
In the B. minax CR, more than ten sequences have the potential to form stem-loop structures with perfect matches and loops of variable size. In addition, several other stem-loop structures with some mismatch in the stems can be predicted. However, obvious stem-loop structures with conserved flanking sequences were not found in the CR of these ten tephritid species. In addition, The B. minax CR does not contain any tRNA-like sequence, but contains two tandem repeats ranging in size from 33 to 45 bp. The sequence TATTAATTTTATTAAA occurred twice and the sequence CCTTTTAAATTTTCC occurred three times. The two repeats were located at positions from 15,325 to 15,357 and from 15,858 to 15,903, respectively. For other tephritid species, we found one tandem repeat in the CR of B. doraslis, B. correcta, B. curcubitae and C. capitata, two in B. philippinensis and B. carambolae, three in B. oleae and B. papaya but none in B. tryoni.