Significant Levels of Sequence Divergence and Gene Rearrangements have Occurred Between the Mitochondrial Genomes of the Wild Mulberry Silkmoth, Bombyx mandarina, and its Close Relative, the Domesticated Silkmoth, Bombyx mori

The metazoan mitochondrial (mt) genome is typically a circular, double-stranded DNA molecule between 14 and 18 kb in size. This molecule encodes 37 genes: 13 protein genes, 2 ribosomal RNA genes, and 22 tRNA genes (Wolstenholme 1992 ). The order of these genes varies among metazoans (Boore 1999 ). The mt genomes of 14 insects have been completely sequenced. These include (1) seven flies (Diptera), Drosophila yakuba (Clary and Wolstenholme 1985— X03240), D. melanogaster (Lewis, Farr, and Kaguni 1995— U37541), D. mauritiana, D. sechellia, D. simulans (Ballard 2000— AF200830 and AF200831 for D. mauritiana, AF200832 for D. sechellia, and AF200834, AF200841, and AF200852 for D. simulans), Ceratitis capitata (Spanos et al. 2000— AJ242872), and Cochliomyia himinivorax (Lessinger et al. 2000— AF260826); (2) two mosquitoes (Diptera), Anopheles quadrimaculatus (Mitchell, Cockburn, and Seawright 1993— L04272) and A. gambiae (Beard, Hamm, and Collins 1993— L20934); (3) the domesticated silkmoth (Lepidoptera), Bombyx mori, the Backokjam strain (Lee et al., unpublished data, but available in GenBank, accession number AF149768); (4) the honeybee (Hymenoptera), Apis mellifera (Crozier and Crozier 1993— L06178); (5) the locust (Orthoptera), Locusta migratoria (Flook, Rowell, and Gellisson 1995— X80245); (6) the kissing bug (Hemiptera), Triatoma dimidiata (Dotson and Beard 2001— AF301594); and (7) the Wallaby louse (Phthiraptera), Heterodoxus macropus (Shao, Campbell, and Barker 2001— AF270939). In all sampled insects other than the Wallaby louse, all mt gene rearrangements involve only tRNA genes. Bombyx mandarina is a probable wild ancestor of the domesticated silkmoth, B. mori. Bombyx mandarina includes significant variation within species, e.g., a difference in chromosome number: whereas 27 chromosomes are seen in the haploid genome of this organism living in Japan (Sasaki, 1889 ), individuals in China and Far Eastern populations carry 28 chromosomes per haploid genome (Kawaguchi, 1928 ; Astaurov, Golisheva, and Roginskaya 1959 ). Note that B. mori carries 28 chromosomes per haploid genome. In this article, the gene rearrangements and sequences are compared between the mt genomes of Japanese B. mandarina and B. mori. A B. mandarina sample was collected in Tsukuba, Ibaraki, Japan. Bombyx mori C108 strain has been maintained in the National Institute of Agrobiological Sciences. The entire mt genomes of B. mandarina and B. mori C108 strain were PCR-amplified in 13 and 11 overlapping fragments, respectively. As template in PCR, genomic DNA was prepared from a pair of silk glands of a final instar larva based on the standard extracting technique for genomic DNA (Sambrook, Fritsch, and Maniatis 1989 , p. 9.16–9.19). PCR primers were initially designed based on the sequences of Simon et al. (1994) , and some of them were modified using the mt nucleotide sequence of B. mori Backokjam strain described earlier. The PCR amplifications were carried out in a 25-μl mixture containing 10 ng of total genomic DNA as a template, 200 nM of each primer, 0.2 mM of each dNTP, 2.5 μl of TaKaRa Ex-Taq™ buffer, and 1.25 U of TaKaRa Ex-Taq™ (TaKaRa). All PCRs were performed using the following conditions: 3 min at 96°C, followed by 35 cycles of 15 s at 96°C, 15 s at 40–45°C, and 1–3 min at 68°C (1 min for each kb) followed by 5 min at 68°C. The two PCR-amplified fragments that contain the A+T–rich region and the alanine tRNA (trnA) gene, respectively, were fractionated by agarose-gel electrophoresis using SeaKem GTG agarose (BMA) and were purified using a QIAquick Gel Extraction kit (Qiagen). These were subcloned into the pGEM-T Vector (Promega) and were sequenced with dRhodamine-terminator chemistry and BigDye-terminator chemistry using two vector-specific primers and internal primers. The remaining PCR-amplified fragments were purified with the QIAquick PCR purification kit (Qiagen) and were directly sequenced using BigDye-terminator chemistry with the PCR primers and internal primers when necessary. All products were analyzed on an ABI 377 sequencer (ABI, Perkin-Elmer). Nucleotide sequences of mt protein–, ribosomal RNA–, and tRNA-coding genes of B. mandarina and B. mori C108 strain were identified through BLAST searches in NCBI and were confirmed by homology search to mt DNA sequences of D. yakuba and B. mori Backokjam strain using Genetyx Mac version 10.1 (SDC). The invertebrate mt genetic code was used for inferring amino acid sequences of the protein-coding genes. The identified sequences of the mt genomes were aligned by Clustal X (Thompson et al. 1997 ) and were analyzed by MEGA (Kumar, Tamura, and Nei 1993 ). The mt genome of the Japanese B. mandarina is a 15,928-bp-long circular DNA and putatively encodes the 37 genes described earlier. The circular mt DNA of B. mori C108 strain has 15,656 bp and the gene order identical to those of B. mandarina and B. mori Backokjam strain. The nucleotide sequence data are in the DDBJ-EMBL-GenBank databases with the accession numbers AB070263 and AB070264 for B. mandarina and B. mori C108, respectively. The common gene arrangement among the three Bombyx mt genomes is almost consistent with that of D. yakuba (Clary and Wolstenholme 1985 ) which can be regarded as the gene order of the ancestral insect (Shao, Campbell, and Barker 2001 ). The only difference in gene arrangement is a translocation of trnM, which moved to a position 5′-upstream of the trnI gene with the same orientation as in D. yakuba. Tayler et al. (1993) showed that tRNA translocation is a frequent event in this area in the evolution of insect mt genomes. The 13 protein-coding genes in mt genomes of B. mandarina and B. mori C108 strain share common initiators: six have ATA-putative initiator codons, five use ATG-putative initiation codons, and nad5 uses an ATT-putative initiator. Cox1 is an exceptional case—its ORF starts at a CGA codon for arginine, and a codon before this triplet is a TAG-stop codon. In Drosophila, a 4-bp initiation codon was proposed for each mt genome: ATAA for D. melanogaster, D. yakuba, D. sechellia, and a haplotype of D. simulans (Clary and Wolstenholme 1985 ; Ballard 2000 ); GTAA for a haplotype of D. mauritiana and two haplotypes of D. simulans; and TTAA for another haplotype of D. mauritiana (Ballard 2000 ). As these proposed 4-bp initiation codons have T and A at the second and third positions, we propose that the Bombyx cox1 translation initiation uses a novel 4-bp initiation codon, TTAG. Eleven protein-coding genes of the B. mandarina and B. mori C108 mt genomes use an identical terminator, TAA. The two remaining genes, cox1 and cox2, use an incomplete termination codon, T, which can generate a complete termination codon, TAA, by polyadenylation of the transcripts. Although the 12 protein-coding genes are identical in size, the B. mandarinacob gene is 6 bp longer than that of B. mori C108 (table 1 ). Adding two deoxyadenosines to 3′-downsteam of the AAA-Lys codon induces a frameshift to generate a 6-bp-longer transcript, which encodes two additional asparagine residues compared with the putative B. moricob product. Sequence rearrangements between the two mt genomes are found in tRNA genes and noncoding regions. Predicted secondary structures are different in tRNAs for leucine (tag), proline, histidine, valine, alanine, and glutamate between the two genomes (fig. 1 ), which are highly associated with insertions-deletions. Although 18 noncoding regions are present in both B. mandarina and B. mori C108 mt genomes, there is variation in the size of these regions between the two mt genomes. The largest noncoding region in the B. mandarina mt genome, a 747-bp sequence between the rrnS and trnM, is in a similar position to the A+T–rich region in other insects thought to contain signals for control of replication and transcription (table 1 ). Indeed, the A+T content of this region is 95.2%, significantly higher than the remainder of the B. mandarina mt genome (χ2df=1 = 95.5, P < 0.005). A tandem triplication of a 126-bp fragment is detected in this region. The first and second fragments share an identical nucleotide sequence, whereas the third one has an A-G transversion compared with the others. On the other hand, the B. mori C108 mt genome has only a 126-bp fragment, showing 99.2% sequence identity to the first and second repetitive units of B. mandarina. This rearrangement largely contributes to the difference in size between the two mt genomes. Five other noncoding regions (Regions 1–5; see Supplementary Material at MBE web site: www.molbiolevol.org) are over 40-bp long in the B. mandarina mt genome. Regions 1–4 are structurally different from those in B. mori, whereas no sequence differences are found in Region 5. Regions 1 and 2 are variable in numbers of TA-dinucleotides, similar to microsatellite sequence divergence. Regions 3 and 4 have much more complicated features: Region 3 shows not only a variable number of TA-dinucleotide repeats but also variation in the repetition number of R(T)n=2–5 (n means the number of repeats; R stands for A for B. mandarina, and G for B. mori). Region 4 of B. mandarina consists mainly of TA-dinucleotide repeats, whereas variable number of TA or TAA repeats coupled with variation in A-repetition number is seen in B. mori. Note that a conservative noncoding sequence, Region 5, contains no typical repetitive structure. Although Regions 1–4 in the Backokjam mt genome also contains repetitive structures in common with C108, the repetitive numbers of respective units are different from those of C108. The pattern of nucleotide differences between the two mt genomes varies by DNA strands. Most of the sequence changes in protein-coding genes are attributable to T-C transitions for the genes encoded by the major strand (31 A-G and 191 T-C transitions in 6,909 bp). In contrast, A-G transitions were more abundant in protein-coding genes on the minor strand (109 A-G and 26 T-C transitions in 4,296 bp). Ballard (2000) indicated a significant bias in the patterns of nucleotide substitutions between the major and minor strands of Drosophila mt genomes, which is consistent with that found in this study. On the other hand, rates of nucleotide changes between the two mt genomes vary by genome region. The rrnS gene, which is located between the rrnL gene and the A+T–rich region, is most conserved among genes or regions (table 1 ), which is also consistent with Ballard (2000) . Nucleotide sequence divergences in the rrnL gene and the A+T–rich region are also low compared with those of protein-coding genes (table 1 ). Excluding protein-coding genes smaller than 600 bp, the five genes surrounding the A+T–rich region and two rrn genes (Class 1: nad2,cox1,cox2,ATP6, and nad1) are more conserved relative to remaining genes (Class 2: cox3,nad4,nad5, and cob). The average substitution rate of Class 1 (0.0296) is significantly smaller than that of Class 2 (0.0440) (χ2df=1 = 14.34, P < 0.005), whereas no significant differences in substitution rates within Classes 1 and 2 are detected (χ2df=4 = 0.767, 0.9 < P and χ2df=3 = 3.88, 0.1 < P < 0.5, respectively). It is notable that, in Drosophila mt genomes, four disjunct regions have a significantly different nucleotide substitution process (Ballard 2000 ). Although B. mori is thought to be derived from B. mandarina, large amounts of sequence divergence between them suggest significant genetic isolation of the current Japanese population of B. mandarina from the population or populations from which B. mori was derived. For instance, the evolutionary distance (D) of the nad5 gene between the two species was estimated to be 0.0393 using Kimura's two parameter method (Kimura 1980 ). If we suppose that the 0.01 D unit corresponds to 3.6 Myr for the nad5 gene according to Su et al. (1998) and Osawa et al. (1999) , the ancestor of B. mori may have separated from that of Japanese B. mandarina about 7.1 MYA. In this study two sets of biased patterns in nucleotide substitution were confirmed, i.e., difference in preference of transitional changes between major and minor DNA strands and variation in degree of nucleotide substitution by genomic regions. Although preference of transition by DNA strands between the two mt genomes is consistent with that of Drosophila, regional variation in nucleotide substitution rate between the two mt genomes show a different feature from that seen in Drosophila mt genomes (Ballard 2000 ). These results suggest that the Bombyx mt genome is not under the common evolutionary process to Drosophila. Frequent sequence rearrangements including insertion-deletion in some tRNA genes and rearrangements driven by repetitive sequences are observed between the two mt genomes despite the stability of gene arrangement within insect mt genomes. The variable number of repetitive units in these regions may enable us to analyze the maternal lineage of B. mori in detail and to carry out population structure analysis. To address the question how B. mori has been established from B. mandarina, combined analysis of nucleotide substitutions with structural rearrangements found in mt and nuclear genes is required, because the nuclear genes have been probably subjected to mostly artificial selection. Such approaches will further contribute to understanding evolutionary dynamics of nuclear and mt genomes of lepidopteran insects. Ross Crozier, Reviewing Editor Present address: Bioinfomatics group, RIKEN GSC, Yokohama, Kanagawa, Japan Abbreviations: atp, ATP synthase F0 subunit; cob, cytochrome b; cox, cytochrome c oxidase subunit; nad, NADH dehydrogenase subunit; rrn, ribosomal RNA; trn, transfer RNA. Keywords: mitochondrial genome Bombyx mandarinaBombyx mori genetic differentiation Address for correspondence and reprints: Kenji Yukuhiro, Insect Genetics and Evolution Department, National Institute of Agrobiological Sciences, Tsukuba, Ibaraki 305-8634, Japan. kygnis@affrc.go.jp Table 1 Nucleotide Differences Between Bombyx mandarina and Bombyx mori C108 Mitochondrial Genes Except for tRNA Genes Table 1 Nucleotide Differences Between Bombyx mandarina and Bombyx mori C108 Mitochondrial Genes Except for tRNA Genes Fig. 1.—Changes of secondary structures of mt tRNAs between B. mandarina and B. mori C108 strain. tRNA secondary structures were predicted by the default search mode of tRNA-SE 1.1 (Lowe and Eddy 1997 ), specifying mt-chloroplast DNA as the source and using the invertebrate mt genetic code. Cloverleaf structures of B. mandarina tRNAs are present, whereas modified parts of those of B. mori are highlighted in boxes. Variable nucleotide sites are marked by bold-italic letters We thank two anonymous reviewers for their valuable comments. This work was supported by the Enhancement of Center of Excellence, Special Coordination Funds for Promoting Science and Technology, and the Science and Technology Agency, Japan.