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

Kenji Yukuhiro*, Hideki Sezutsu*,1, Masanobu Itoh{dagger}, Koichi Shimizu{dagger} and Yutaka Banno{ddagger}

*Insect Genetics and Evolution Department, National Institute of Agrobiological Sciences, Tsukuba, Ibaraki, Japan;
{dagger}Department of Applied Biology, Faculty of Textile Science, Kyoto Institute of Technology, Kyoto, Japan;
{ddagger}Division of Silkworm, Institute of Genetic Resources, Faculty of Agriculture, Kyushu University, Fukuoka, Japan

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 1992Citation ). The order of these genes varies among metazoans (Boore 1999Citation ).

The mt genomes of 14 insects have been completely sequenced. These include (1) seven flies (Diptera), Drosophila yakuba (Clary and Wolstenholme 1985Citation —X03240), D. melanogaster (Lewis, Farr, and Kaguni 1995Citation —U37541), D. mauritiana, D. sechellia, D. simulans (Ballard 2000Citation —AF200830 and AF200831 for D. mauritiana, AF200832 for D. sechellia, and AF200834, AF200841, and AF200852 for D. simulans), Ceratitis capitata (Spanos et al. 2000Citation —AJ242872), and Cochliomyia himinivorax (Lessinger et al. 2000Citation —AF260826); (2) two mosquitoes (Diptera), Anopheles quadrimaculatus (Mitchell, Cockburn, and Seawright 1993Citation —L04272) and A. gambiae (Beard, Hamm, and Collins 1993Citation —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 1993CitationL06178); (5) the locust (Orthoptera), Locusta migratoria (Flook, Rowell, and Gellisson 1995CitationX80245); (6) the kissing bug (Hemiptera), Triatoma dimidiata (Dotson and Beard 2001CitationAF301594); and (7) the Wallaby louse (Phthiraptera), Heterodoxus macropus (Shao, Campbell, and Barker 2001CitationAF270939). 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, 1889Citation ), individuals in China and Far Eastern populations carry 28 chromosomes per haploid genome (Kawaguchi, 1928Citation ; Astaurov, Golisheva, and Roginskaya 1959Citation ). 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 1989Citation , p. 9.16–9.19). PCR primers were initially designed based on the sequences of Simon et al. (1994)Citation , 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-TaqTM buffer, and 1.25 U of TaKaRa Ex-TaqTM (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. 1997Citation ) and were analyzed by MEGA (Kumar, Tamura, and Nei 1993Citation ).

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 1985Citation ) which can be regarded as the gene order of the ancestral insect (Shao, Campbell, and Barker 2001Citation ). 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)Citation 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 1985Citation ; Ballard 2000Citation ); GTAA for a haplotype of D. mauritiana and two haplotypes of D. simulans; and TTAA for another haplotype of D. mauritiana (Ballard 2000Citation ). 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. mandarina cob 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. mori cob product.


View this table:
[in this window]
[in a new window]
 
Table 1 Nucleotide Differences Between Bombyx mandarina and Bombyx mori C108 Mitochondrial Genes Except for tRNA Genes

 
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.



View larger version (23K):
[in this window]
[in a new window]
 
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 1997Citation ), 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

 
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 ({chi}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)Citation 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)Citation . 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) ({chi}2df=1 = 14.34, P < 0.005), whereas no significant differences in substitution rates within Classes 1 and 2 are detected ({chi}2df=4 = 0.767, 0.9 < P and {chi}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 2000Citation ).

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 1980Citation ). If we suppose that the 0.01 D unit corresponds to 3.6 Myr for the nad5 gene according to Su et al. (1998)Citation and Osawa et al. (1999)Citation , 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 2000Citation ). 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.

Acknowledgements

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.

Footnotes

Ross Crozier, Reviewing Editor

1 Present address: Bioinfomatics group, RIKEN GSC, Yokohama, Kanagawa, Japan Back

Abbreviations: atp, ATP synthase F0 subunit; cob, cytochrome b; cox, cytochrome c oxidase subunit; nad, NADH dehydrogenase subunit; rrn, ribosomal RNA; trn, transfer RNA. Back

Keywords: mitochondrial genome Bombyx mandarina Bombyx mori genetic differentiation Back

Address for correspondence and reprints: Kenji Yukuhiro, Insect Genetics and Evolution Department, National Institute of Agrobiological Sciences, Tsukuba, Ibaraki 305-8634, Japan. kygnis{at}affrc.go.jp Back

References

    Astaurov B. L., M. D. Golisheva, I. S. Roginskaya, 1959 Chromosome complex of Ussuri geographical race of Bombyx mandarina M. with spacial reference to the problem of the origin of the domesticated silkworm, Bombyx mori L Cytology 1:327-332 [in Russian]

    Ballard J. W. O., 2000 Comparative genomics of mitochondrial DNA in members of the Drosophila melanogaster subgroup J. Mol. Evol 51:48-63[ISI][Medline]

    Beard C. B., D. M. Hamm, F. H. Collins, 1993 The mitochondrial genome of mosquito Anopheles gambiae: DNA sequence, genome organization, and comparison with mitochondrial sequences of other insects Insect Mol. Biol 2:103-104[Medline]

    Boore J. L., 1999 Animal mitochondrial genomes Nucleic Acids Res 27:1761-1780

    Clary D. O., D. R. Wolstenholme, 1985 The mitochondrial DNA molecular of Drosophila yakuba: nucleotide sequence, gene organization, and genetic code J. Mol. Evol 22:252-271[ISI][Medline]

    Crozier R. H., Y. C. Crozier, 1993 The mitochondrial genome of the honeybee Apis mellifera: complete sequence and genome organization Genetics 133:97-117[Abstract/Free Full Text]

    Dotson E. M., C. B. Beard, 2001 Sequence and organization of the mitochondrial genome of the Chagas disease vector, Triatoma dimidiata Insect Mol. Biol 10:205-215[ISI][Medline]

    Flook P. K., C. H. F. Rowell, G. Gellisson, 1995 The sequence, organization and evolution of Locusta migratoria mitochondrial genome J. Mol. Evol 41:928-941[ISI][Medline]

    Kawaguchi E., 1928 Zytologische Untersuchungen am Seidenspinner Und seinen Verwandten. I. Gametogenese von Bombyx mori L. und B. mandarina M. und ihrer Bestarde. Z. Zellforsch. Mikrosk Anat 7:154-156

    Kimura M., 1980 The simple method for estimating evolutionary rate of base substitutions thorough comparative studies of nucleotide sequences J. Mol. Evol 16:111-120[ISI][Medline]

    Kumar S., K. Tamura, M. Nei, 1993 MEGA: molecular evolutionary genetics analysis. Version 1.01 Institute of Molecular Evolutionary Genetics, The Pennsylvania State University, Philadelphia

    Lessinger A. C., A. C. Martins Junqueira, T. A. Lemos, E. L. Kemper, F. R. da Silva, A. L. Vettore, P. Arruda, A. M. Azeredo-Espin, 2000 The mitochondrial genome of the primary screwworm fly Cochliomyia hominivorax (Diptera: Calliphoridae) Insect Mol. Biol 9:521-529[ISI][Medline]

    Lewis D. L., C. L. Farr, L. S. Kaguni, 1995 Drosophila melanogaster mitochondrial DNA: completion of the nucleotide sequence and evolutionary comparisons Insect Mol. Biol 4:263-278[ISI][Medline]

    Lowe T. M., S. R. Eddy, 1997 tRNA Scan-SE: a program for improved detection of transfer RNA genes in genomic sequences Nucleic Acid Res 25:955-964[Abstract/Free Full Text]

    Mitchell S. E., A. F. Cockburn, J. A. Seawright, 1993 The mitochondrial genome of Anopheles quadrimaculatus species A: complete nucleotide sequence and gene organization Genome 36:1058-1073[ISI][Medline]

    Osawa S., Z.-H. Su, C.-G. Kim, M. Okamoto, O. Tominaga, Y. Imamura, 1999 Evolution of the carabid ground beetles Adv. Biophys 36:65-106[ISI][Medline]

    Sambrook J., F. E. Fritsch, T. Maniatis, 1989 Molecular cloning Cold Spring Harbor Laboratory Press, New York

    Sasaki C., 1889 On the affinity of our wild and domestic silkworms Ann. Zool. Jpn 2:2.

    Shao R., N. J. H. Campbell, S. Barker, 2001 Numerous gene rearrangements in the mitochondrial genome of the Wallaby louse, Heterodoxus macropus (Phthiraptera) Mol. Biol. Evol 18:858-865[Abstract/Free Full Text]

    Simon C., F. Frati, A. Beckerbach, B. Crespi, H. Liu, P. Flook, 1994 Evolution, weighting and phylogenetic utility of mitochondrial gene sequences and a compilation of the conserved polymerase chain reaction primers Ann. Entomol. Soc. Am 87:651-701[ISI]

    Spanos L. G., G. Koutroumbas, M. Kotsyfakis, C. Louis, 2000 The mitochondrial genome of the Mediterranean fruit fly, Ceratiius capitata Insect Mol. Biol 9:139-144[ISI][Medline]

    Su Z.-H., O. Tominaga, M. Okamoto, S. Osawa, 1998 Origin and diversification of hindwingless Damaster ground beetles within the Japanese island as deduced from mitochondrial ND5 gene sequences (Coleoptera: Carabidae) Mol. Biol. Evol 15:1026-1039[Abstract]

    Tayler M. F. J., S. W. McKechnie, N. Pierce, M. Kreitman, 1993 The lepidopteran mitochondrial control region: Structure and evolution Mol. Biol. Evol 10:1259-1272[Abstract]

    Thompson J. D., T. J. Gibson, F. Plewniak, F. Jeanmougin, D. G. Higgins, 1997 The ClustalX windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools Nucleic Acids Res 24:4876-4882

    Wolstenholme D. R., 1992 Animal mitochondrial DNA: structure and evolution Int. Rev. Cytol 141:173-216[ISI][Medline]

Accepted for publication March 28, 2002.