Analysis on the Complete Mitochondrial Genome of Andraca theae (Lepidoptera: Bombycoidea).

The bombycid moth, Andraca theae (Matsumura) (Lepidoptera: Bombycoidea) is an important pest of tea in southeastern China. In the present study, the mitochondrial genome (mitogenome) of A. theae was amplified by polymerase chain reaction and sequenced. The complete mitogenome of A. theae, encoding 37 genes, was 15,737 bp in length (Genbank no. KX365419), and consisted of 13 protein-coding genes (PCGs), 22 tRNA genes, 2 ribosomal RNA genes and an adenine (A) + thymine (T)-rich region (AT-rich region). The gene order of A. theae mitogenome was typical for Lepidoptera mitogenomes. Except for cox1, which started with CGA, all other 12 PCGs started with ATN. Eleven of the 13 PCGs ended with TAA, expect for cox1 and cox2, which ended with a single T. The maximum likelihood method and the Bayesian method were used to analyze the phylogenetic relationship among 22 representative bombycoid species with a matrix consisting of the 13 PCGs of the mitogenomes of the 22 species. The topological structures of the two phylogenetic trees we constructed were almost identical, with the results indicating that the bombycid species, including A. theae, clustered into a single clade with a bootstrap value of 58% and a posterior probability of 0.98. The phylogenetic relationship among the Bombycoidea species analyzed was Lasiocampidae + (Bombycidae + (Saturniidae + Sphingidae)) which was supported by a high bootstrap value of 100% and a posterior probability of 1.00.

The moth species, Andraca theae (Matsumura), (Lepidoptera: Bombycoidea) is a serious pest of tea [Camellia sinensis (L.) Kuntze] (Theales: Theaceae) that is widely distributed throughout the tea-growing areas in southeastern China (Anhui, Hunan, Guangdong, Guizhou, and Taiwan Provinces) and Nepal (Wang et al. 2015). Tea, as mankind’s most commonly consumed drink, has become increasingly popular throughout the world. Of the many pests occurring on tea, most are rarely even reported, and little is known regarding any of the species genetic characteristics. In order to effectively manage A. theae in the field, additional basic information is needed on the species genetic features and phylogenetic position.

Lepidopteran mitochondrial DNA has characteristics typically found in the mitochondrial genome (mitogenome) of many other invertebrates, encoding 37 mitochondrial genes consisting of 13 protein-coding genes (PCGs), 22 tRNA genes, 2 ribosomal RNA (rRNA) genes and a noncoding region (AT-rich region) (Wolstenholme 1992, Lewis et al. 1995, Zhang et al. 1995, Inohira et al. 1997, Boore 1999). We were able to extract an extensive amount of information from the mitogenome of A. theae, including gene arrangement, base composition, genetic codon variation, and the secondary structures of tRNA and rRNA genes. This information could be used to understand the unique features of each individual. Mitogenome sequences are also widely used in population genetics, reconstruction of phylogenetic relationships, and evolutionary genomics and biology (Avise 1987, Ballard and Rand 2005, Cameron and Whiting 2008).

The family Bombycidae s. lat. is a relatively well-known lepidopteran taxon, belonging to the superfamily Bombycoidea. As currently defined, the family contains ± 40 genera and roughly 350 species (Lemaire and Minet 1999). Bombycid mitogenomes have rarely been studied. To date, only four species’ complete mitochondrial genome have been reported, Bombyx mandarina (Moore) (NC_003395), Bombyx huttoni Westwood (NC_026518), Bombyx mori (Linnaeus) (NC_002355), and Rondotia menciana Moore (NC_021962) (Yukuhiro et al. 2002, Peng et al. 2015, Kong and Yang. 2015).

In our study, the complete mitogenome of A. theae is sequenced and described for the first time and compared with other bombycoid species. The phylogenetic arrangement of the bombycoid species was analyzed based on mitochondrial data. The intent of this study was to contribute to understanding the phylogenetic and evolutionary relationships among the Bombycoidea.

Materials and Methods

Sampling and DNA Extraction

Larvae of A. theae were collected from a tea plantation located in Panxian Town (Liupanshui City, Guizhou Province, China; 25° 49′ N, 104° 38′ E) in October 2015 and were provided by Hui-Zhu Wang. Species identification of the specimens was based on Wang et al. (2011, 2012, 2015). Genomic DNA was extracted from fresh larvae using a Wizard Genomic DNA Purification Kit (Promega, Beijing, China) according to the manufacturer's instruction.

Polymerase Chain Reaction Amplification and Sequencing

Fourteen pairs of primers were used to amplify the entire mitogenome sequence of A. theae. The polymerase chain reaction (PCR) amplification was performed in a 25-μL reaction containing 0.2 μL of rTaq (TaKaRa Co., Dalian, China), 1-μL DNA template, 2.5-μL 10× rTaq buffer (Mg2+ free), 2.5 μL 25 mM MgCl2, 2.0-μL dNTPs, and 0.5 μL in each primer. The PCR amplifications were performed under the following cycling parameters: 5 min at 94 °C, followed by 35 cycles of 30 s at 94 °C, 30 s at 47–58 °C, and 1–2.5 min at 72 °C, with a subsequent 10-min final extension at 72 °C. PCR products were resolved by electrophoresis in a 1.0% agarose gel and were extracted using a DNA Gel Extraction Kit (Bioteke, Beijing, China). For the exact procedure used, see Ma et al. (2016).

Sequence Analysis

The fragments of A. theae mitogenome were assembled using the program Geneious version 8.1.2 (Kearse et al. 2012). When the sequencing of mitogenome was completed, it was annotated manually and automatically. The method described by Cameron and Whiting (2008) was used for the written annotation, and the online-program MITOS (Bernt et al. 2013) was used to accomplish the automated annotation. The PCGs boundaries were confirmed by the ORF finder (http://www.ncbi.nlm.nih.gov/orf/gorf.html). MEGA version 6.0 (Tamura et al. 2013) was used to determine the starting and termination codons of the PCGs. Identification of the tRNA genes was verified using the tRNAscan-SE program (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe and Eddy 1997) and the program MITOS. Unidentified tRNAs were compared with sequences from other species. The two rRNA subunits coding gene (rrnL and rrnS) were identified by a NCBI Internet BLAST search. The AT content and codon usage were calculated using Geneious version 8.1.2 (Kearse et al. 2012). The formulas of the skews of the compositions are as follows: AT skew = [A − T]/[A + T]; GC skew = [G − C]/[G + C] (Junqueira et al. 2004).

Phylogenetic Analysis

The complete mitochondrial genomes of 21 bombycoid species were obtained from GenBank (Table 1) to help with determining the relationship among Bombycoidea taxa. Twenty-one bombycoid species plus A. theae were considered as ingroups, while the geometrids Biston panterinaria (Bremer & Grey) and Phthonandria atrilineata (Butler) were used as outgroups (Kong and Yang 2015). The nucleic acid region from all 13 PCGs were excluded the start and stop codons, concatenated as a matrix and aligned using the program Geneious version 8.1.2. Gblock 0.91b with default settings was used with conserved regions of the putative amino acids (Castresana 2000). For likelihood ratio tests, jModeltest 2.1.5 and Akaike Information Criterion (Ronquist and Huelsenbeck 2003) were used to determine the best-fitting model of each PCG., We then chose two analysis approaches to construct a phylogenetic tree, the maximum likelihood (ML) and Bayesian inference (BI) analyses. The ML method was conducted using the program raxmlGUL 1.5 (https://sourceforge.net/projects/raxmlgui/) (Silvestro and Michalak 2012), in which the analysis setting chosen was “ML + rapid bootstrap” and the “number of bootstrap replicates” was 1,000. The support values of the ML tree were evaluated via a bootstrap test with 100 iterations. The program MrBayes 3.1.2 (http://morphbank.Ebc.uu.SE/mrbayes/) (Ronquist et al. 2012) was used to perform BI analysis (Darriba et al. 2012), with the MCMC analysis run for 1,000,000 generations and a burn-in series of 1,000.

Table 1.

Source and information for phylogenetic analysis

Species	Family	Subfamily	Accession number	References
Actias artemis aliena	Saturniidae	Saturniinae	KF927042	Park et al. (2016)
Actias selene	Saturniidae	Saturniinae	NC_018133	Liu et al. (2012)
Andraca theae	Bombycidae	Bombycinae	KX365419	This study
Antheraea frithi	Saturniidae	Saturniinae	NC_027071	Unpublished
Antheraea pernyi	Saturniidae	Saturniinae	NC_004622	Liu et al. (2008)
Antheraea yamamai	Saturniidae	Saturniinae	NC_012739	Kim et al. (2009)
Apatelopteryx phenax	Lasiocampidae	Lasiocampinae	KJ508055	Timmermans et al. (2014)
Attacus atlas	Saturniidae	Saturniinae	NC_021770	Chen et al. (2014)
Bombyx huttoni	Bombycidae	Bombycinae	NC_026518	Peng et al. (2015)
Bombyx mandarina	Bombycidae	Bombycinae	NC_003395	Yukuhiro et al. (2002)
Bombyx mori	Bombycidae	Bombycinae	NC_002355	Unpublished
Dendrolimus punctatus	Lasiocampidae	Pinarinae	NC_027156	Unpublished
Dendrolimus spectabilis	Lasiocampidae	Pinarinae	NC_025763	Tang et al. (2014)
Dendrolimus tabulaeformis	Lasiocampidae	Pinarinae	NC_027157	Unpublished
Eriogyna pyretorum	Saturniidae	Saturniinae	NC_012727	Jiang et al. (2009)
Manduca sexta	Sphingidae	Sphinginae	NC_010266	Cameron and Whiting (2008)
Rondotia menciana	Bombycida	Bombycinae	NC_021962	Kong and Yang (2015)
Samia canningi	Saturniidae	Saturniinae	NC_024270	Shantibala et al. (2016)
Samia cynthia cynthia	Saturniidae	Saturniinae	KC812618	Sima et al. (2013)
Samia cynthia ricini	Saturniidae	Saturniinae	NC_017869	Kim et al. (2012)
Saturnia boisduvalii	Saturniidae	Saturniinae	NC_010613	Hong et al. (2008)
Sphinx morio	Sphingidae	Sphinginae	NC_020780	Kim et al. (2013)

Results

Genome Organization and Base Composition

The complete mitochondrial genome of A. theae, which was shown to be a closed-circular molecule of 15,737 bp in length, was deposited in GenBank (NCBI) with accession number KX365419. It encoded 37 genes, including 13 PCGs (cox1-3, nad1-6, nad4l, atp6, atp8, and cytb), 22 tRNA genes, and 2 rRNA genes, and a putative AT-rich region (Table 2). Twenty-four genes were transcribed on the major strand (J-strand), with the remaining 13 genes transcribed on the minor strand (N-strand). The gene distribution in the mitochondrial genome was conserved as in other lepidopteran mitochondrial genomes (Fig. 1). We found that the mitochondrial genome of A. theae was loose, except in the AT-rich region. There were 509-bp intergenic nucleotides that were dispersed in 22 pairs of neighboring genes with their length varying from 1 to 80 bp. The longest spacer region was between trnQ and nad2 (80 bp). The overlapping nucleotides existed in eight pairs of neighboring genes ranging from 1 to 29 bp, with the longest overlap region crossing trnF and nad5 (29 bp).

Table 2.

Annotation and gene organization of the A. theae mitogenome

Gene	Position (bp)	Length (bp)	Direction	Intergenc nucleotides (intergenic nucleotides)	Anticodons	Start/stop codons	AT%
trnM	1–66	66	J		CAT		77.27
trnI	67–131	65	J	0	GAT		75.38
trnQ	129–197	69	N	−3	TTG		81.16
nad2	255–1,268	1,014	J	57		ATT/TAA	81.26
trnW	1,269–1,341	73	J	0	TCA		79.45
trnC	1,334–1,398	65	N	−8	GCA		80.00
trnY	1,408–1,473	66	N	9	GTA		74.24
cox1	1,486–3,019	1,534	J	12		CGA/T	68.38
trnL(UUR)	3,020–3,086	67	J	0	TAA		73.13
cox2	3,121–3,802	682	J	34		ATG/T	71.14
trnK	3,803–3,873	71	J	0	CTT		73.24
trnD	3,954–4,024	71	J	80	GTC		90.14
atp8	4,025–4,198	174	J	0		ATT/TAA	90.80
atp6	4,192–4,869	678	J	−7		ATG/TAA	76.04
cox3	4,882–5,670	789	J	12		ATG/TAA	71.14
trnG	5,673–5,739	67	J	2	TCC		88.06
nad3	5,737–6,093	357	J	0		ATC/TAA	79.66
trnA	6,107–6,175	69	J	13	TGC		81.16
trnR	6,253–6,315	63	J	77	TCG		77.78
trnN	6,330–6,394	65	J	14	GTT		75.38
trnS(AGN)	6,398–6,462	65	J	3	GCT		80.00
trnE	6,464–6,532	69	J	1	TTC		91.30
trnF	6,535–6,602	68	N	2	GAA		86.76
nad5	6,574–8,340	1,767	N	−29		ATT/TAA	79.40
trnH	8,341–8,404	64	N	0	GTG		82.81
nad4	8,417–9,757	1,341	N	12		ATG/TAA	77.33
nad4l	9,808–10,098	291	N	50		ATG/TAA	81.44
trnT	10,103–10,166	64	J	−8	TGT		82.81
trnP	10,167–10,229	63	N	0	TGG		80.95
nad6	10,253–10,765	513	J	23		ATA/TAA	81.85
cytb	10,769–11,920	1,152	J	3		ATG/TAA	74.08
trnS(UCN)	11,947–1,2013	67	J	26	TGA		82.09
nad1	12,035–12,973	939	N	21		ATT/TAA	74.76
trnL(CUN)	12,975–13,043	69	N	1	TAG		81.16
rrnL	13,100–14,477	1,378	N	56			83.53
trnV	14,479–14,543	65	N	1	TAC		81.51
rrnS	14,544–15,322	779	N	0			83.31
AT rich region	15,323–15,737	415	–	0			91.81

Fig. 1.

Circular map of the A. theae mitochondrial genome. Gene names are annotated using standard abbreviations; single letters are indicated based on the International Union of Pure and Applied Chemistry-International Union of Biochemistry(IUPUC-IUB) abbreviation for the corresponding amino acid.

The base compositions of the major strand of the A. theae mitogenome were as follows: A = 40.28%, T = 38.01%, C = 13.90% and G = 7.81%, with a total AT content is 78.29%, which was heavily biased toward A and T bases. AT- and GC-skews of the entire major strand of A. theae were calculated as: AT-skew= 0.029, GC-skew=  −0.281.

Protein-Coding Genes

The mitogenome of A. theae, encoded 13 PCGs, with 3,731 amino acids in total. The major strand (J-strand) contained nad2, cox1, cox2, atp8, atp6, cox3, nad3, nad6, and cytb, whereas the minor strand (N-strand) included nad5, nad4, nad4l, and nad1. Except cox1, which started with CGA, the other 12 PCGs started with ATN. The nad2, cox2, atp8, nad3, nad5, nad4, nad4l, nad6, and nad1 PCGs used ATA as the starting codon, whereas atp6, cox3, and cytb began with ATG. For the termination codon, 11 of the 13 PCGs ended with TAA, whereas cox1 and cox2 ended with a single T. Aside from the start and termination codons, the sequence of the concentrated 13 PCGs was 11,193 bp in length, the content of A + T was 75.92% and was much higher than G + C according to the summarized codon usage. The Relative Synonymous Codon Usage values in the PCGs of the A. theae mitogenome are shown in Table 3 and Figure 2.

Table 3.

Codon usage of the PCGs in A. theae

Codon (aa)	n	%	RSCU	Codon (aa)	n	%	RSCU
UUU(F)	362	9.70	1.85	UAU(Y)	160	4.29	1.74
UUC(F)	30	0.80	0.15	UAC(Y)	24	0.64	0.26
UUA(L)	382	10.24	4.2	UAA(*)	0	0	0
UUG(L)	61	1.63	0.67	UAG(*)	0	0	0
CUU(L)	38	1.02	0.42	CAU(H)	47	1.26	1.38
CUC(L)	10	0.27	0.11	CAC(H)	21	0.56	0.62
CUA(L)	50	1.34	0.55	CAA(Q)	56	1.50	1.84
CUG(L)	5	0.13	0.05	CAG(Q)	5	0.13	0.16
AUU(I)	398	10.67	1.8	AAU(N)	211	5.66	1.63
AUC(I)	44	1.18	0.2	AAC(N)	48	1.29	0.37
AUA(M)	214	5.74	1.62	AAA(K)	94	2.52	1.83
AUG(M)	51	1.37	0.38	AAG(K)	9	0.24	0.17
GUU(V)	88	2.36	2.11	GAU(D)	53	1.42	1.74
GUC(V)	8	0.21	0.19	GAC(D)	8	0.21	0.26
GUA(V)	60	1.61	1.44	GAA(E)	61	1.63	1.65
GUG(V)	11	0.29	0.26	GAG(E)	13	0.35	0.35
UCU(S)	110	2.95	2.76	UGU(C)	26	1.61	1.63
UCC(S)	20	0.54	0.5	UGC(C)	6	0.16	0.38
UCA(S)	56	1.50	1.4	UGA(W)	85	2.28	1.79
UCG(S)	6	0.16	0.15	UGG(W)	10	0.27	0.21
CCU(P)	51	1.39	1.65	CGU(R)	9	0.24	0.69
CCC(P)	43	1.15	1.39	CGC(R)	2	0.05	0.15
CCA(P)	30	0.80	0.97	CGA(R)	38	1.02	2.92
CCG(P)	0	0	0	CGG(R)	3	0.08	0.23
ACU(T)	82	2.20	2.01	AGU(S)	35	0.94	0.88
ACC(T)	23	0.62	0.56	AGC(S)	4	0.11	0.1
ACA(T)	54	1.45	1.33	AGA(S)	86	2.31	2.16
ACG(T)	4	0.11	0.1	AGG(S)	2	0.05	0.05
GCU(A)	81	2.17	2.51	GGU(G)	39	1.05	0.8
GCC(A)	13	0.35	0.4	GGC(G)	10	0.27	0.21
GCA(A)	32	0.86	0.99	GGA(G)	91	2.44	1.87
GCG(A)	3	0.08	0.09	GGG(G)	55	1.47	1.13

RSCU = Relative Synonymous Codon Usage.

Fig. 2.

Codon usage in the A. theae mitogenome. (A) Codons per thousand (CDspT) indicates the codons used in coding amino acids per thousand codons. Codon families are given on the x-axis. (B) Relative synonymous codon usage (RSCU).

Transfer RNA and Ribosomal RNA Genes

The A. theae mitochondrial genome included 22 tRNA genes ranging from 63 bp to 73 bp in size, displaying a high AT content of 80.97%. Among the 22 tRNA genes, 14 tRNAs were found in the J-strand and 8 tRNAs in the N-strand (Table 2; Fig. 1). The figure shows that all tRNAs, except trnS (AGN), trnR, trnH, and trnT, displayed the typical clover-leaf secondary structure. The trnS (AGN) lacked the dihydrouridine (DHU) arm and formed a simple loop. The secondary structures of the other three tRNA genes (trnR, trnH, and trnT) were atypical. The trnR was missing the DHU loop, whereas the TΨC loop was absent in the trnH and trnT genes (Fig. 3).

Fig. 3.

Predicted secondary structure of tRNA genes in the A. theae mitogenome. Dashes (-) indicate Watson–Crick base pairing.

The two rRNA genes, rrnL and rrnS, were 1,378 bp and 779 bp in size, respectively, with an AT content of 83.53% and 83.31%. Both of the two rRNA genes were mapped on the N-strand. The rrnL was located between trnL (CUN) and trnV, and the rrnS were situated in trnV and trnM.

Features in the AT-Rich Region

The AT-rich region of A. theae was 415 bp in length with a 91.8% AT content and was located between rrnS and trnM (Fig. 4). A conserved structure consisting of the motif “ATAGA,” which was followed by a 19-bp poly-T stretch, was observed in the downstream 18 bp of rrnS. The three microsatellites, “(TA)6,” “(TA)5,” and “(TA)8”, that were found in this region were located 159, 123, and 70 bp upstream of trnM, respectively. The poly-A stretch, which is located upstream of trnM in some lepidopteran species, was noted and was 15 bp in size.

Fig. 4.

Structure of the A + T-rich region of A. theae.

Phylogenetic Relationships

The phylogenetic relationships within the Bomycoidea, reconstructed with the concatenated nucleotides sequences of 13 PCGs of 24 mitogenomes, are shown in Figures 5 and 6. The 22 bombycids, including the 21 bombycoid species mitogenomes that were downloaded from Genbank plus A. theae, represent four families belonging to the Bombycoidea: Bombycidae, Lasiocampidae, Saturniidae, and Sphingidae. The two phylogenetic trees based on ML and BI analyses showed that the phylogenetic relationships in the Bombycoidea were Lasiocampidae + (Bombycidae + (Saturniidae + Sphingidae)), which was supported by a high bootstrap value of 100% and a posterior probability of 1.00. The bombycid species were A. theae + (R. menciana + (B. huttoni + (B. mandarina + B. mori))), which was supported as a monophyletic group by a bootstrap value of 58% and a posterior probability of 0.98.

Fig. 5.

ML phylogram constructed using 13 PCGs of mitogenomes with partitioned models. The scale bar indicates the number of substitutions per spot The values indicated at each node specify the bootstrap percentage of 1,000 replicates.

Fig. 6.

BI phylogram constructed using 13 PCGs of mitogenomes with unpartitioned models. The scale bar indicates the number of substitutions per spot The value indicated at each node specifies the bootstrap percentage of 1,000 replicates.

Discussion

The order of the gene organizations of the A. theae complete mitogenome, trnM-trnI-trnQ-nad2, was identical to the other lepidopteran mitogenomes (Flook et al. 1995). The AT content of the A. theae mitogenome (78.29%) was the same as in other bombycid species (Table 4). The AT skewness of A. theae mitogenome was 0.029, suggesting that the occurrence of A was more than T. We found that the AT skewness of A. theae was higher than in Apatelopteryx phenax, Dendrolimus spectabilis, Samia cynthia cynthia, and Sphinx morio (0.001–0.028) and lower than in Attacus atlas, Bombyx huttoni, Bombyx mori, and Rondotia menciana (0.039–0.072). However, 11 of the 22 bombycoid mitochondrial genomes were T-skewed (−0.005 to − 0.045). The active selection on nucleotide composition may be related to overcoming background mutation pressures (Meiklejohn et al. 2007).

Table 4.

Comparison of the nucleotide composition and skewness of Bombycoidea taxa

Insect species	Whole genome			PCGs codon		rrnL		rrnS		A+T rich
	Size	A+T %	AT skewness	Size	A+T %	Size	A+T %	Size	A+T %	Size	A+T %
Actias artemis aliena	15,243	78.61	−0.013	11,187	76.85	1,363	83.13	777	84.30	328	91.77
Actias selene	15,236	78.91	−0.023	11,184	77.30	1,364	83.58	762	83.99	339	87.91
Andraca theae	15,737	78.29	0.029	11,193	75.92	1,378	83.53	779	83.31	415	91.81
Antheraea frithi	15,338	80.19	−0.018	11,175	78.65	1,380	83.99	777	84.56	333	89.19
Antheraea pernyi	15,566	80.16	−0.021	11,196	78.47	1,369	83.86	775	84.13	552	90.40
Antheraea yamamai	15,338	80.29	−0.022	11,187	78.89	1,380	83.99	776	84.52	334	89.52
Apatelopteryx phenax	15,552	80.33	0.027	11,207	78.42	–	–	–	–	–	–
Attacus atlas	15,282	79.30	0.039	11,181	77.77	1,368	84.80	777	83.14	357	90.48
Bombyx huttoni	15,638	81.80	0.072	11,157	80.04	1,399	85.28	777	86.74	512	95.12
Bombyx mandarina	15,928	81.68	−0.045	11,166	79.59	1,377	84.75	783	85.95	747	95.18
Bombyx mori	15,643	81.32	0.059	11,142	79.51	1,375	84.36	783	85.57	449	95.39
Dendrolimus punctatus	15,411	79.46	0.029	11,181	77.55	1,461	84.60	779	84.98	320	92.81
Dendrolimus spectabilis	15,411	79.50	0.028	11,120	77.65	1,357	83.49	628	84.08	465	92.04
Dendrolimus tabulaeformis	15,411	79.53	0.029	11,178	77.63	1,459	84.72	778	84.70	320	92.81
Eriogyna pyretorum	15,327	80.82	−0.031	11,193	79.35	1,338	84.60	778	84.45	358	92.18
Manduca sexta	15,516	81.79	−0.005	11,154	80.24	1,391	85.26	777	85.71	324	90.35
Rondotia menciana	15,301	78.86	0.050	11,190	77.04	1,365	83.37	782	84.40	357	91.04
Samia canningi	15,384	79.88	−0.007	11,193	78.38	1,358	84.02	779	83.95	361	91.14
Samia cynthia cynthia	15,345	79.86	0.008	11,178	78.33	1,359	84.17	778	84.19	359	91.09
Samia cynthia ricini	15,384	79.78	−0.006	11,196	78.26	1,358	84.02	779	83.83	361	90.86
Saturnia boisduvalii	15,360	80.62	−0.024	11,202	79.10	1,391	84.76	774	84.11	330	91.52
Sphinx morio	15,299	81.17	0.001	11,148	79.78	1,379	84.55	773	85.25	316	92.72

The AT content in the PCGs of A. theae was 75.92%, which was the lowest value in all the bombycoid species. Twelve of the 13 PCGs started with “ATN,” whereas cox1 started with “CGA.” Having “CGA” as the starting codon may be a synapomorphic or diagnostic character in Lepidoptera (Kim et al. 2009). The incomplete stop codon of a single T was found in A. theae mitochondrial cox1 and cox2. This is a common phenomenon of mitochondrial genes, where the incomplete stop codons can be completed (TAA) by the mRNA process of polyadenylation (Anderson et al. 1981, Liu et al. 2013, Yang et al. 2013).

The longest four noncoding regions were located between trnQ-nad2 (57 bp), trnK-trnD (80 bp), trnR-trnN (77 bp), and trnL (CUN)-rrnL (56 bp). The intergenic spacer between trnQ and nad2 showed limited sequence conservation among the studied lepidopteran species (Cameron and Whiting 2008). The spacer of trnk-trnD and trnL (CUN)-rrnL had demonstrated a high AT content, and included (TA)10, (TA)9, and (TTA)2 microsatellite-like regions that were observed in the other lepidopteran mitogenomes (Cameron and Whiting 2008). There were only a few overlappings between mitochondrial genes of A. theae. The 7-bp overlapping nucleotides (ATGATAA) between atp8 and atp6, which are a common feature in all lepidopteran mitogenomes, were observed in the mitochondrial genome of A. theae. The position of maximum overlapping was not conserved in other lepidopteran species, and the overlapping was situated between many different pairs of genes. The maximum overlapping of A. theae was located from trnF to nad5.

The AT-rich region is thought to be the site of gene replication and the initiation of genome transcription (Boore 1999, Taanman 1999). In the mitochondrial genome of lepidopteran species, the motif “ATAGA” followed by poly-T was a conserved position. The motif was found in the AT-rich region of the A. theae mitogenome, and the poly-T was deemed the origin of the minor strand replication. The origin of the major strand was less conserved (Saito et al. 2005, Cameron and Whiting 2008).

The mitochondrial genome is widely used in phylogenetic analyses. After using the ML and BI methods in our study, the species of Bombycidae were found to cluster into one clade. Although the ML method value found for the Bombycidae in our study was low, similar topological structures of two separate analyses could clarified the monophyly of Bombycidae. In the study by Minet, (1994), the family Bombycidae was divided into four subfamilies, the Apatelodinae, Phiditiinae, Prismostictinae, and Bombycinae, with the genus Andraca placed in the subfamily Primostictinae (Oberthueriinae) based on morphological characters. After incorporating molecular datasets into the phylogeny, however, the classification was drastically altered. According to Regier et al. (2008), the four former bombycid subfamilies were now separated into distinct families. The former Primostictinae were now divided into two clades, one containing the tribe Oberthueriini, where the genus Andraca was placed. This clade was shown with a high bootstrap value to cluster with the Mirinidae. This outcome was consistence with Zwick et al. (2011). Although Regier et al. (2009), came to a different conclusion that grouped the two bombycid subfamilies, Phiditiinae and Prismostictinae, combining them together with the Carthaeidae and Anthelidae. It is obvious from the previous results (Minet, 1994; Regier et al. 2008, 2009; Zwick et al. 2011) that the phylogenetic relationships within the Bombycidae, especially the Bombycinae, are confusing. It is worthy of note, although, that the studies discussed above were based on nuclear genes., Our study, on the other hand, was focused on mitochondrial data, and in our study, the subfamily Oberthuerinae is placed as the sister group to the Bombycinae. However, there remains a problem with the phylogenetic relationships of families among the Bombycoidea in our study where the group “Lasiocampidae + (Bombycidae + (Saturniidae + Sphingidae)),” is consistent with some past studies (Liu et al. 2013), while differing from the findings of other previous studies where the families group as “Lasiocampidae + (Saturniidae+ (Bombycidae + Sphingidae))” (Timmermans et al. 2014). The explanation for the difference may be the result of incorporating the complete mitogenome. Significantly, the bootstrap value of the Bombycidae clade and the Saturniidae + Sphingidae clade group at a low level at around 50%. That emphasizes that the relationships in the Bombycoidea remain unsettled and will need further attention in the future.