The Complete Mitochondrial Genome of Coptotermes 'suzhouensis' (syn. Coptotermes formosanus) (Isoptera: Rhinotermitidae) and Molecular Phylogeny Analysis.

Coptotermes suzhouensis (Isoptera: Rhinotermitidae) is a significant subterranean termite pest of wooden structures and is widely distributed in southeastern China. The complete mitochondrial DNA sequence of C. suzhouensis was analyzed in this study. The mitogenome was a circular molecule of 15,764 bp in length, which contained 13 protein-coding genes (PCGs), 22 transfer RNA genes, two ribosomal RNA genes, and an A+T-rich region with a gene arrangement typical of Isoptera mitogenomes. All PCGs were initiated by ATN codons and terminated by complete termination codons (TAA), except COX2, ND5, and Cytb, which ended with an incomplete termination codon T. All tRNAs displayed a typical clover-leaf structure, except for tRNASer(AGN), which did not contain the stem-loop structure in the DHU arm. The A+T content (69.23%) of the A+T-rich region (949 bp) was higher than that of the entire mitogenome (65.60%), and two different sets of repeat units (A+B) were distributed in this region. Comparison of complete mitogenome sequences with those of Coptotermes formosanus indicated that the two taxa have very high genetic similarity. Forty-one representative termite species were used to construct phylogenetic trees by maximum likelihood, maximum parsimony, and Bayesian inference methods. The phylogenetic analyses also strongly supported (BPP, MLBP, and MPBP = 100%) that all C. suzhouensis and C. formosanus samples gathered into one clade with genetic distances between 0.000 and 0.002. This study provides molecular evidence for a more robust phylogenetic position of C. suzhouensis and inferrs that C. suzhouensis was the synonymy of C. formosanus.

Introduction

Termites (Isoptera) (Eggleton et al. 2007) comprise more than 3,000 species in approximately 283 genera (Inward et al. 2007a, Krishna et al. 2013, Cheng 2014). Of all termite species recognized globally, only 183 are significant pest species known to damage buildings (Krishna et al. 2013, Chouvenc et al. 2015), of which the genus Coptotermes contains the largest number (18) of pest species (Rust and Su 2012). Due to their destructive effect on wooden structures and their essential role in decomposition and nutrient recycling, research on the biological characteristics, classification, identification, and phylogenetic relationships of termites has received increasing attention (Wolstenholme 1992, Brody et al. 2010, Hausberger et al. 2011, Korb et al. 2015, Chouvenc et al. 2016, Bourguignon et al. 2016, Rocha et al. 2017). In China, 4 families, 44 genera, and 473 species of termites have been recorded (Huang et al. 2000, Cheng 2014). The Rhinotermitidae are widely distributed in China, with approximately 200 described species in seven genera, including representatives of major pest genera Coptotermes, Reticulitermes, and Heterotermes (Huang et al. 2000, Cheng 2012). The difficulty in identifying termites is recognized by many termite classification experts, especially in the genus Coptotermes (Emerson 1952, Watson and Gay 1991, Huang et al. 2000, Chouvenc et al. 2016). Taxonomic limits between species in genus Coptotermes, are poorly established, which has resulted in different names being used for the same species (Chouvenc et al. 2016). This has created additional difficulties in the identification of termites, which were already been complicated by the limited morphological features available for identifications. The number of species in the genus Coptotermes is frequently revised, as synonyms of other species are recognised (Engel and Krishna 2004, Krishna et al. 2013, Cameron 2014a). Most species of Coptotermes are considered to be similar to each other, and their identification is thus very difficult. Kambhampati (2000) also thought that many synonyms were found in the nomenclature of termites based on morphological classifications. The classification and identification of Coptotermes have been revised based on molecular data (Cai and Chen 1964, Huang et al. 2000, Xu et al. 2009, Cheng 2012). To date, 22 species of Coptotermes have been recorded in China (Krishna et al. 2013). Coptotermes suzhouensis (Isoptera: Rhinotermitidae) is very similar to Coptotermes formosanus (Isoptera: Rhinotermitidae), a major invasive species, in morphology. C. suzhouensis is widely distributed in the regions of southeastern China and is a pest of wooden buildings (Xia and He 1986, Huang et al. 2000, Li 2002, Cheng 2012). There has always been controversy regarding the relationship between the two species (C. suzhouensis and C. formosanus).

In the study, we determined the complete mitogenome sequence of C. suzhouensis and compared it to available mitogenomes of other termites using phylogenetic analyses. In doing so, we provide robust molecular evidence for the taxonomic status of C. suzhouensis and reveal the evolutionary relationships of C. suzhouensis and C. formosanus.

Materials and Methods

Collection and Storage

Specimens of C. suzhouensis were collected from colonies of termites living in wooden buildings in Feixi County in Hefei, Anhui Province, China (31°42′31″N, 117°10′38″E). The specimens were preserved in 95% ethanol and stored at −20°C until DNA extraction.

Morphological Identification of Termites

Soldier specimens of all the populations collected were identified by the Hefei Termite Control Institute based on their morphological characteristics (Xia and He 1986, Huang et al. 2000).

DNA Isolation, Polymerase Chain Reaction, and Sequencing

Whole genomic DNA was extracted from the heads of 20 soldiers using One-tube General Sample DNAup Kit (Sangon Biotech, Shanghai, China) for polymerase chain reaction (PCR), following manufacturer’s instructions. Mitochondrial genome sequences of the termites related to C. suzhouensis were downloaded from the NCBI database. ClustalX ver1.8 was used for alignment, and Premier Primer 5 software was applied to design five pairs of primers from conserved regions (Table 1). Fragments were amplified by long distance PCR using Takara LA Taq (Takara Bio, Japan) with the following cycling conditions: an initial denaturation for 3 min at 94°C, followed by 35 cycles of denaturation at 94°C (30 s), annealing at 52−58°C (30 s; Table 1), elongation at 72°C (10 min); and a final extension period 72°C (10 min). PCR products were checked using 1% agarose gels and were sequenced at MAP Biotech Company (Shanghai, China) with ABI3730 Genetic Analyzer.

Table 1.

The primers for PCR in this study

Primers	Sequences (5ʹ–3ʹ)	Positions (5ʹ–3ʹ)	Size of PCR product (bp)	Annealing temperatures (°C)
Primer-1-F	TATCGCCATACCATCACTACGACTCCTA	3,340–3,367	5,154	55
Primer-1-R	TGCTCCCCCTTCTCTTAATCTTCTCGGT	8,493-8,466
Primer-2-F	GAACCAAAGCAGACACAGGAGTAGGAGC	7,481–7,508	4,742	58
Primer-2-R	TGGGCTTCGTGCTTTGGCTCAGACTATC	12,222-12,195
Primer-3-F	AGAAACCAACTCCGATTCCCCCTCAGCA	11,988–12,015	7,427	58
Primer-3-R	GTCGTCCTGGTGTGGCGTCTGTTTTTAC	3,650-3,623
Primer-COX3-F	ATTCCACCAATACGACAACAGCCTA	4,872–4,896	582	52
Primer-COX3-R	GAGAAGTGTAGGGCTGCTTGTCGTA	5,453-5,429
Primer-Cytb-F	GACATCAATACCGCATTTTCCAGAG	10,638–10,662	850	52
Primer-Cytb-R	GTCGTGCTCCGATTCAGGTAAGTAG	11,487-11,463

Sequence Assembly, Annotation, and Analysis

Sequences were assembled and aligned with complete mitogenomes from C. formosanus in Sequencher 4.1.4. The location, size, and coding direction of each gene, including 13 protein coding genes (PCGs), 22 tRNA genes, and two rRNA genes were determined with DOGMA (Wyman et al. 2004). Secondary structures of tRNA were predicted using MITOS (Bernt et al. 2013) and tRNAscan-SE Search Server v.2.0 online (Lowe and Chan 2016). Amino acid composition and coding region for each PCG was identified with ORF Finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html). Nucleotide composition statistics (excluding stop codons) and relative synonymous codon usage (RSCU) of 13 PCGs were calculated with MEGA 5.1 (Tamura et al. 2014). Composition skew analysis was calculated according to formulas AT skew = [A−T]/[A+T] and GC skew = [G−C]/[G+C], respectively (Perna and Kocher 1995).

The Reconstruction of Phylogenetic Trees

Along with the C. suzhouensis mitochondrial genome, 41 termite mitogenomes and two outgroup species (Shelfordella lateralis and Periplaneta australasiae: Blattodea) were used in phylogenetic analysis. Sampling and sequence availability used in this study are summarized in Table 2. Concatenated amino acid sequences from the 13 PCGs was used in phylogenetic analysis, with maximum likelihood (RAxML7.03; Stamatakis et al. 2008), maximum parsimony (PAUP 4.0b10; Swofford 2002), and Bayesian inference methods (Mr Bayes v.3.1.2) (Huelsenbeck and Ronquist 2001). Modeltest ver3.06 (Posada and Crandall 1998) was used to infer best-fitting model for Bayesian inference (BI) and maximum likelihood (ML) analysis. The GTR + I + G model was selected based on the Akaike information criterion (Akaike 1974). Bayesian inference was performed with the following settings: four MCMC chains (one cold chain and three hot chains) for 10,000,000 generations until the average standard deviation of split frequencies reached a value less than 0.01. Bayesian posterior probabilities (BPP) were calculated from the sample points after the MCMC algorithm had started to converge. In ML and MP analyses, a heuristic search with 100 random addition replicates was applied. BPP values were mapped onto the tree, and nodal support for ML&MP was assessed using nonparametric bootstrapping (Felsenstein 1985) in PAUP for the MP analyses (MPBP) and RAxML (Stamatakis et al. 2008) for ML (MLBP) using 1,000 pseudoreplicates each.

Table 2.

Detailed information of the termites species analyzed in this study

Family	Species	Size (bp)	Accession number	References
Rhinotermitidae	Reticulitermes aculabialis	16,475	KP334994.1	Wang et al. (2015)
	Reticulitermes chinensis	15,925	KM216388.1	Chen et al. (2016)
	Reticulitermes flavipes	16,565	EF206314.1	Cameron and Whiting (2007)
	Reticulitermes hageni	16,590	EF206320.1	Cameron and Whiting (2007)
	Reticulitermes santonensis	16,567	EF206315.1	Cameron and Whiting (2007)
	Reticulitermes virginicus	16,513	EF206318.1	Cameron and Whiting (2007)
	Reticulitermes labralis	15,914	KU877221.1	Wang et al. (2016)
	Reticulitermes speratus kyushuensis	15,898	KY484910.1	Lee et al. (2017)
	Reticulitermes sp.	14,907	KU925239.1	Bourguignon et al. (2017)
	Reticulitermes grassei	14,910	KU925237.1	Bourguignon et al. (2017)
	Reticulitermes flaviceps	16,485	KX712090.1	Chen et al. (2016)
	Heterotermes sp.	16,370	JX144936.1	Cameron et al. (2012)
	Heterotermes cf.occiduus	14,919	KU925230.1	Bourguignon et al. (2017)
	Heterotermes cf. paradoxus	14,904	KU925225.1	Bourguignon et al. (2017)
	Heterotermes nr. tenuis	14,944	KU925228.1	Bourguignon et al. (2017)
	Heterotermes platycephalus	14,919	KU925231.1	Bourguignon et al. (2017)
	Heterotermes tenuis	14,940	KU925233.1	Bourguignon et al. (2017)
	Heterotermes validus	14,922	KU925235.1	Bourguignon et al. (2017)
	Coptotermes acinaciformis raffrayi	14,897	KU925196.1	Bourguignon et al. 2017
	Coptotermes acinaciforis	14,896	KU925198.1	Bourguignon et al. (2017)
	Coptotermes amanii	14,894	KU925200.1	Bourguignon et al. (2017)
	Coptotermes lacteus	16,326	JX144934.1	Cameron et al. (2012)
	Coptotermes frenchi	14,912	KU925204.1	Bourguignon et al. (2017)
	Coptotermes gestroi	14,919	KU925205.1	Bourguignon et al. (2017)
	Coptotermes heimi	14,908	KU925206.1	Bourguignon et al. (2017)
	Coptotermes kalshoveni	14,889	KU925210.1	Bourguignon et al. (2017)
	Coptotermes michaelseni	14,900	KU925212.1	Bourguignon et al. 2017
	Coptotermes remotus	14,742	KU925213.1	Bourguignon et al. 2017
	Coptotermes sjoestedti	14,899	KU925217.1	Bourguignon et al. 2017
	Coptotermes sepangensis	14,905	KU925215.1	Bourguignon et al. (2017)
	Coptotermes travians	14,892	KU925222.1	Bourguignon et al. (2017)
	Coptotermes testaceus	15,752	KR872938.1	Li et al. (2016)
	Coptotermes formosanus	14,908	KU925203.1	Bourguignon et al. (2017)
	Coptotermes formosanus	16,324	AB626147.1	Tokuda et al. (2012)
	Coptotermes formosanus	16,326	AB626146.1	Tokuda et al. (2012)
	Coptotermes formosanus	1,326	AB626145.1	Tokuda et al. (2012)
	Coptotermes suzhouensis	15,764	MG000963	present study
Termitidae	Macrotermes barneyi	15,940	JX050221.1	Wei et al. (2012)
	Acanthotermes acanthothorax	15,231	KP026280.1	Bourguignon et al. (2015)
	Ancistrotermes pakistanicus	15,299	KP026267.1	Bourguignon et al. (2015)
	Macrotermes natalensis	16,325	KM405637.1	Meng et al. (2016)
	Macrotermes subhyalinus	16,351	JX144937.1	Cameron et al. (2012)
Serritermitidae	Serritermes serrifer	14,783	KP026264.1	Bourguignon et al. (2015)
	Glossotermes occulatus	14,791	KP026291.1	Bourguignon et al. (2015)
Kalotermitidae	Glyptotermes satsumensis	15,611	KP026257.1	Bourguignon et al. (2015)
	Cryptotermes secundus	15,695	KP026283.1	Bourguignon et al. (2015)
Mastotermitidae	Mastotermes darwiniensis	15,487	JX144929.1	Cameron et al. (2012)
Outgroup	Shelfordella lateralis	15,601	KU684413.1	Cheng et al. (2016)
	Periplaneta australasiae	15,605	KX640825.1	Ma et al. (2017)

Results

Sequencing and Organization of Mitochondrial Genome

The mitochondrial genome of C. suzhouensis is a typical circular DNA molecule of 15,764 bp in length (GenBank: MG000963). The mitochondrial DNA (mtDNA) consisted of 13 PCGs (ATP6, ATP8, COX1-3, ND1-6, ND4L, and Cytb), 2 rRNA genes (srRNA and lrRNA), 22 tRNA genes, and a noncoding A+T-rich region (Fig. 1, Table 4). Nucleotide composition was A/T (65.60%) biased and composed as follows: A = 6,891 (43.71%), T = 3,450 (21.89%), G = 1,854 (11.76%), C = 3,569 (22.64%), which is commonly observed in termites’ mitogenomes (Bourguignon et al. 2015; Table 3). In total, 23 genes (9 PCGs and 14 tRNAs) were located on the majority strand (H-strand) and the others (4 PCGs, 8 tRNAs, and 2 rRNAs) were located on the minority strand (L-strand; Fig. 1). The order and the orientation of the genes were identical to those previously reported from other Coptotermes species, which retain the ancestral insect arrangement.

Fig. 1.

Circular map of the mitogenome of C. suzhouensis. Genes encoded on the H-strand (clockwise orientation) are colored in the outside. Genes encoded on the L-strand (anticlockwise orientation) are colored in the inside.

Table 4.

Annotation and gene organization of the C. suzhouensis mitogenome

Gene	Coding strand	Region	Size (bp)	Intergenic nucleotide	Overlapping nucleotide	Anticodon	Start codon	Stop codon
tRNAIle	H	1–66	66			GAT (30–32)
tRNAGln	L	64–133	70		3	TTG (101–103)
tRNAMet	H	155–223	69	21		CAT (187–189)
ND2	H	224–1,261	1,038				ATG	TAA
tRNATrp	H	1,267–1,334	68	5		TCA (1,298–1,300)
tRNACys	L	1,327–1,395	69		8	GCA (1,363–1,365)
tRNATyr	L	1,406–1,481	76	10		GTA (1,442–1,444)
COX1	H	1,483–3,030	1,548	1			ATT	TAA
tRNALeu(UUR)	H	3,042–3,107	66	11		TAA (3,071–3,073)
COX2	H	3,116–3,800	685	8			ATG	T
tRNALys	H	3,801–3,871	71			CTT (3,831–3,833)
tRNAAsp	H	3,871–3,935	65		1	GTC (3,902–3,904)
ATP8	H	3,936–4,094	159				ATA	TAA
ATP6	H	4,088–4,768	681		7		ATG	TAA
COX3	H	4,768–5,556	789		1		ATG	TAA
tRNAGly	H	5,563–5,630	68	6		TCC (5,597–5,599)
ND3	H	5,631–5,984	354				ATA	TAA
tRNAAla	H	5,994–6,056	63	9		TGC (6,023–6,025)
tRNAArg	H	6,062–6,130	69	5		TCG (6,091–6,093)
tRNAAsn	H	6,136–6,203	68	5		GTT (6,167–6,169)
tRNASer(AGN)	H	6,201–6,272	72		3	GCT (6,228–6,230)
tRNAGlu	H	6,270–6,333	64		3	TTC(6,299–6,301)
tRNAPhe	L	6,345–6,412	68	11		GAA (6,376–6,378)
ND5	L	6,413–8,138	1,726				ATA	T
tRNAHis	L	8,142–8,206	65	3		GTG (8,173–8,175)
ND4	L	8,220–9,554	1,335	13			ATG	TAA
ND4L	L	9,548–9,835	288		7		ATG	TAA
tRNAThr	H	9,839–9,903	65	3		TGT (9,869–9,871)
tRNAPro	L	9,903–9,971	69		1	TGG (9,938–9,940)
ND6	H	9,973–10,464	492	1			ATA	TAA
Cytb	H	10,464–11,595	1,132		1		ATG	T
tRNASer(UCN)	H	11,596–11,668	73			TGA (11,629–11,631)
ND1	L	11,688–12,626	939	19			ATG	TAA
tRNALeu(CUN)	L	12,633–12,699	67	6		TAG (12,668–12,670)
lrRNA(16S)	L	12,700–14,019	1,320
tRNAVal	L	14,021–14,087	67	1		TAC (14,056–14,058)
srRNA(12S)	L	14,089–14,815	727	1
A+T-rich	nc	14,816–15,764	949
Repeat A1	nc	14,913–14,978	66
Repeat A2	nc	15,102–15,167	66
Repeat B1	nc	15,156–15,717	562
Repeat B2	nc	15,718–15,755	38

Table 3.

Nucleotides composition of the C. suzhouensis mitochondrion in different regions

Feature	Proportion of nucleotides
	A%	T%	G%	C%	A+T%	AT-skew	GC-skew	Siz e(bp)
Whole genome	43.71	21.89	11.76	22.64	65.60	0.33	−0.32	15,764
PCGs	43.52	21.06	12.05	23.37	64.58	0.35	−0.32	11,166
tRNA genes	39.05	26.84	13.82	20.29	65.89	0.19	−0.19	1,498
srRNA genes	45.25	20.63	10.45	23.66	65.89	0.37	−0.39	727
lrRNA gene	48.56	22.05	8.56	20.83	70.61	0.38	−0.42	1,320
A+T-rich region	43.94	25.29	11.70	19.07	69.23	0.27	−0.24	949

The mitogenome of C. suzhouensis harbored a total of 139 bp of intergenic spacers, made up of 19 individual regions ranging in size from 1 to 21 bp. There were 35 base pairs of overlapping regions total, at 10 intergenic positions, which ranged in size from 1 to 8 bp (Table 4).

Protein-Coding Genes

The total length of the 13 PCGs was 11,166 bp, representing 70.83% of the entire mitochondrial genome. PCGs used ATN as initiation codon, all had ATG as start codon, except for ATP8, ND3, ND5, ND6 (ATA), and COX1 (ATT). Ten PCGs used the standard stop codon TAA, whereas COX2, ND5, and Cytb genes used a single T nucleotide. Codon usage of the PCGs exhibited an AT bias with an A+T composition of 64.58% (Table 3). It was found that the relative synonymous codon usage values of NNU/NNA codons was essentially greater than of NNC/NNG, indicating higher U+A bias at third condons. Our analysis showed that UUU (Phe), CUA (Leu), AUA (Met), and AUU (Ile) were the most frequently used codons, accounting for 19.32% of all the codons (Table 5).

Table 5.

Codon usage in 13 PCGs of the C. suzhouensis mitochondrial genome

Codon (aa)	Count	%	RSCU	Codon (aa)	Count	%	RSCU	Codon (aa)	Count	%	RSCU	Codon (aa)	Count	%	RSCU
UUU(F)	201	5.42	1.24	UCU(S)	96	2.59	2.12	UAU(Y)	102	2.75	1.29	UGU(C)	38	1.02	1.55
UUC(F)	122	3.29	0.76	UCC(S)	19	0.51	0.42	UAC(Y)	56	1.51	0.71	UGC(C)	11	0.3	0.45
UUA(L)	103	2.78	1.16	UCA(S)	102	2.75	2.25	UAA(*)	0	0	0	UGA(W)	68	1.83	1.33
UUG(L)	146	3.93	1.64	UCG(S)	10	0.27	0.22	UAG(*)	0	0	0	UGG(W)	34	0.92	0.67
CUU(L)	61	1.64	0.68	CCU(P)	45	1.21	1.22	CAU(H)	22	0.59	0.61	CGU(R)	21	0.57	1.4
CUC(L)	17	0.46	0.19	CCC(P)	17	0.46	0.46	CAC(H)	50	1.35	1.39	CGC(R)	0	0	0
CUA(L)	189	5.09	2.12	CCA(P)	82	2.21	2.23	CAA(Q)	53	1.43	1.51	CGA(R)	35	0.94	2.33
CUG(L)	19	0.51	0.21	CCG(P)	3	0.08	0.08	CAG(Q)	17	0.46	0.49	CGG(R)	4	0.11	0.27
AUU(I)	155	4.18	1.11	ACU(T)	46	1.24	0.76	AAU(N)	65	1.75	0.87	AGU(S)	32	0.86	0.71
AUC(I)	124	3.34	0.89	ACC(T)	44	1.19	0.73	AAC(N)	85	2.29	1.13	AGC(S)	6	0.16	0.13
AUA(M)	172	4.63	1.42	ACA(T)	141	3.8	2.33	AAA(K)	60	1.62	1.52	AGA(S)	79	2.13	1.75
AUG(M)	70	1.89	0.58	ACG(T)	11	0.3	0.18	AAG(K)	19	0.51	0.48	AGG(S)	18	0.49	0.4
GUU(V)	114	3.07	1.82	GCU(A)	50	1.35	1.09	GAU(D)	36	0.97	1.04	GGU(G)	91	2.45	1.45
GUC(V)	18	0.49	0.29	GCC(A)	27	0.73	0.59	GAC(D)	33	0.89	0.96	GGC(G)	14	0.38	0.22
GUA(V)	96	2.59	1.53	GCA(A)	94	2.53	2.04	GAA(E)	63	1.7	1.47	GGA(G)	117	3.15	1.86
GUG(V)	23	0.62	0.37	GCG(A)	13	0.35	0.28	GAG(E)	23	0.62	0.53	GGG(G)	29	0.78	0.46

Transfer RNA and Ribosomal RNA Genes

The mtDNA contained 22 tRNA genes, ranging in size from 63 bp (tRNA) to 76 bp (tRNA) in length, and were AT biased (65.89%). The secondary structures of tRNAs were predicted by MITOS online (Bernt et al. 2013) and tRNAscan-SE 2.0 online (Lowe and Chan 2016). The genes that could not be detected by the softwares were determined through comparison with published mitogenomes of Coptotermes (Cameron et al. 2012, Tokuda et al. 2012). All tRNAs could be folded into typical cloverleaf secondary structures, with the exception of the tRNA, which was lacking a stable stem-loop structure in the DHU arm, as observed in most insects (Cameron 2014b; Fig. 2) .

Fig. 2.

Secondary structures for 22 tRNAs of C. suzhouensis mitogenome predicted by the the MITOS and tRNAscan-SE 2.0 online.

As typically observed in other insect mitogenomes, two rRNA genes (srRNA and lrRNA) were found on the L strand of the mitogenome, which were located between tRNA and tRNA, tRNA and CR region, respectively. The rRNAs of C. suzhouensis were 1,320 bp for lrRNA and 727 bp for srRNA in length, and the A+T content of the two genes was 70.61 and 65.89%, respectively.

Control Region

The 949-bp control region of C. suzhouensis was located between srRNA and tRNA with an A+T content of 69.23%, which was higher than that of the complete mitogenome (65.60%). There were two different sets of repeat units in the CR zone (A+B repeats). The A repeats contained two identical units A1 and A2 (66 bp). The B repeats consisted of one complete unit B1 (562 bp) and a partial unit B2 (38 bp).

Homology Analysis of Mitochondrial Sequences of C. suzhouensis

To explore the phylogenetic potential of the determined sequence, we performed multiple alignment of the mitogenomes determined for C. suzhouensis and Rhinotermitidae. The nucleic acid similarity rate between these taxa was found to rang from 85 to 99%. Further, C. suzhouensis shared the highest homology with C. formosanus, with nucleic acid similarity of more than 99%, whereas the deduced amino acids similarity of individual PCGs ranged from 99.72 to 100%. The base composition of Rhinotermitidae mitogenomes showed a high degree of similarity with A+T biased (61.77−66.33%).

Differences in mitochondrial sequences between C. suzhouensis and C. formosanus were shown in Table 6. Comparison of complete mitogenome sequences with the three populations of C. formosanus Shiraki (AB626145.1−AB626147.1) showed five to seven site differences and a long fragment deletion (B1 repeat region of 562 bp; 3.47% sequence divergence; Table 6). Variable sites between C. suzhouensis and C. formosanus were located in four PCGs (COX1, COX2, Cytb, and ND1), tRNA, lrRNA (16S), and Repeat A1 region; variation in COX1 and COX2 was synonymous. However, two sites in Cytb and ND1 were nonsynomyous. The A to G (C. suzhouensis) at position 10,722 led to serine (AGA) to glycine (GGA) transformation, whereas the C to T at 11,978 led to valine (GTT) to isoleucine (ATT) transformation, respectively. In contrast, the differences between C. suzhouensis and C. formosanus (KU925203) were greater, with 31 differences, of which 13 were synonymous, and 5 were nonsynomyous. Comparative analyses suggested that the complete mtDNA of C. suzhouensis and C. formosanus are highly similar, consistant with them being the same species.

Table 6.

Differences in the nucleotides of C. suzhouensis and C. formosanus mitochondrial genomes

Position	Gene	C. suzhouensis	C. formosanus	C. formosanus	C. formosanus	C. formosanus
		(MG000963) 15,764 bp	(AB626145.1) 16,326 bp	(AB626146.1) 16,326 bp	(AB626147.1) 16,324 bp	(KU925203.1) 14,908 bp
482	ND2	T(Phe)	T	T	T	C(Leu)
646		A(Met)	A	A	A	G(Met)
1,578	COX1	T(Leu)	T	T	T	C(Leu)
1,990		A(Met)	A	A	A	G(Val)
2,061		G(Leu)	G	G	G	A(Leu)
2,334		T(Phe)	C(Phe)	C(Phe)	C(Phe)	C(Phe)
3,056	tRNA Leu(UUR)	C	C	C	C	T
3,063		A	A	A	A	G
3,337	COX2	A(Val)	A	A	A	G(Val)
3,397		G(Thr)	A(Thr)	A(Thr)	A(Thr)	A(Thr)
3,412		T(Ala)	T	T	T	A(Ala)
3,824	tRNA Lys	A	A	C	A	A
4,223	ATP6	A(Met)	A	A	A	T(Leu)
4,348		C(Phe)	C	C	C	T(Phe)
4,465		T(His)	T	T	T	C(His)
5,580	tRNA Gly	A	A	A	—	—
5,581		A	A	A	—	—
5,928	ND3	T(Leu)	T	T	T	C(Leu)
5,981		A(Lys)	A	A	A	C(Asn)
6,020	tRNA Ala	A	A	A	A	G
7,598	ND5	T(Met)	T	T	T	C(Val)
7,719		C(Met)	C	C	C	T(Met)
7,866		C(Ser)	C	C	C	A(Ser)
8,523	ND4	C(Trp)	C	C	C	T(Trp)
10,722	Cytb	G(Gly)	G	A(Ser)	G	G
11,687–11,688	intergenic region (tRNA Ser(UCN) — ND1)	—	—	—	—	TTAC
11,978	ND1	T(Ile)	C(Val)	C(Val)	T	T
13,040	lrRNA(16S)	C	T	T	T	T
13,660		C	C	C	T	T
14,514	srRNA(12S)	A	A	A	A	—
14,590		A	A	A	A	T
14,664		T	T	T	T	C
14,765		C	C	C	C	T
14,963	CR	T	C	C	T	—
14,816–15,764		(14,816–15,764)	(14,816–16,326)	(14,816–16,326)	(14,816–16,326)	—
		949bp	1,511bp	1,511bp	1,511bp	—
		(A1,A2,B1,B3)	(A1,A2.B1,B2,B3)	(A1.,A2,B1,B2,B3)	(A1,A2,B1,B2,B3)	—

Nonsynonymous substitutions are indicated with the bold font; corresponding amino acids are shown in parentheses. Positions are relevant to MG000963. Deletions are indicated as —.

Phylogenetic Relationships

Phylogenetic trees were built from 13 PCGs using three different methods (ML, MP, and BI). The topological structure of the phylogenetic tree was largely consistent, and the classification of the family was clear. The analyzed species (41 termites) were divided into five major clades with the basic framework: (Mastotermitidae + (Kalotermitidae + (Serritermitidae + (Termitidae + Rhinotermitidae)))) (Fig. 3). In the phylogenetic tree, Coptotermes, Heterotermes, and Reticulitermes formed a monophyletic group, and the relationship between the three genera was Reticulitermes + (Heterotermes + Coptotermes), which is consistent with morphological data, and findings of previous molecular studies (Bourguignon et al. 2015, Bourguignon et al. 2017, Inward et al. 2007b, Lo et al. 2004).

Fig. 3.

Phylogenetic trees inferred with the amino acid sequences of 13 PCGs of the mitogenome of 41 termites species. S. lateralis (KU684413.1) and P. australasiae (KX640825.1) were used as outgroups. Numbers above the branches represent Bayesian posterior probabilities and bootstrap branch support for Maximum likelihood and Maximum parsimony, respectively. Nodes, which all support rates were 100%, were marked with an asterisk.

Phylogenetic analyses suggested that C. suzhouensis and C. formosanus were clustered in one branch with strong support (BPP = 100%, MLBP = 100%, and MPBP = 100%), and these two groups formed a sister group to (C. kalshoveni + (C. remotus + C. sepangensis)). The genetic distances between C. suzhouensis and several C. formosanus samples was 0.000 (AB626145.1−AB626147.1), while the genetic distance between C. suzhouensis and C. formosanus was 0.002 (KU925203.1), which indicated that C. suzhouensis shared considerably close evolutionary relationships with C. formosanus.

Discussion

Termites play an important role in nutrient cycling and decomposition but are also often considered pest species as they may damage wooden buildings. Identifying termite species, especially those in the genus Coptotermes, is very difficult, and the taxonomic validity of many named Coptotermes species remains unclear (Chouvenc et al. 2015). mtDNA has been extensively used as an informative molecular marker for diverse evolutionary studies of animals, including in molecular evolution, phylogenetics, and population genetics (Gissi et al. 2008, Cameron 2014b, Qin et al. 2015). Thus, molecular tools may aid in the identification of termite species and resolves the relationships between C. suzhouensis and C. formosanus.

The present study analyzed the complete mtDNA sequence of C. suzhouensis. The gene arrangement of the mitogenome was identical to that of Coptotermes, as well as consistent with the ancestral arrangement of the insect mitogenome, indicating that it has been conserved during the evolution of these insects (Wei et al. 2010). Genetic analysis indicated interspecies genetic distance within Rhinotermitidae was 0.029−0.186. In Coptotermes, the difference among individuals within species (0.000−0.019) was lower than that among species (0.029−0.116), whereas the average genetic distance between C. suzhouensis and C. formosanus was 0.000−0.002, indicating that the two taxa share considerably close evolutionary relationships with each other (consistent with the findings of phylogenetic analyses; Fig. 3, Supp Table 1 [online only]). Our phylogenetic analyses support the viewpoint that C. suzhouensis was the synonymy of C. formosanus.

In Rhinotermitidae, the initiation and termination of 13 PCGs were essentially consistent. The majority of PCGs utilized canonical start codons (ATN) and stop codons (TAA, TAG, TA, or T). For C. suzhouensis, the stop codons of 13 PCGs were TAA, except for 3 PCGs (COX2, ND5, and Cytb), which were terminated with the partial stop codon T. Incomplete stop codons are common features of insect mitogenomes (Dietrich and Brune 2014, Hervé and Brune 2017), and it has been proposed that the complete termination codons are generated by the post-transcriptional polyadenylation (Ojala et al. 1981).

The A+T rich region, known for the initiation of replication, is located between srRNA and tRNA in Rhinotermitidae. And, the A+T content (69.23%) of the control region (949 bp) in C. suzhouensis was higher than that of the entire mitogenome (65.60%), which is consistent with other Rhinotermitidae mitogenomes. The control region has the highly variability in the base composition and structure within Rhinotermitidae. In termites, the complicated double repeat units were first found in Reticulitermes, consisting of short (type-A) repeats adjacent to the srRNA end and long (type-B) repeats adjacent to the tRNA end (Supp Table 2 [online only]). Cameron (2012) considered that the repeat units are involved in the replication-mediated processes and are a synapomorphic feature of Neotermitoidea, secondarily lost A-type repeats in higher Termitidae (non-macrotermitine termitids). In C. suzhouensis, the repeat units consist of two full A units, one full B, and one partial B unit (A-A-B-Bp). In addition, the second A and full B units overlap 12 bp. Compared to C. formosanus (A-A-B-B-Bp), a deletion of one complete B repeat was found in the CR zone of C. suzhouensis, and the phenomenon of repeat unit reduction was also observed in other termites [Reticulitermes virginicus (Isoptera: Rhinotermitidae), Cameron et al. 2007] . Comparative analysis of the control region within Rhinotermitidae revealed two conserved sequences ‘AATCCTAAACTTATCT’ (located at B1 repeat region: from 15,223 to 15,238) and ‘AGATAAGTTTAGGATT’ (located at B1 repeat region: from 15,254 to 15,269) in the control regions of C. suzhouensis (Supp Fig. 1 [online only]) that formed the hairpin loop structure. The hairpin loop known as RGC (rare genomic change) is considered to have a high degree of sequence conservation in all termite species and speculated to be the origin of replication for the mitogenome (Saito et al. 2005). Further, a poly-T stretch (located at B1 repeat region: from 15,238 to 15,245) was found in the loop of hairpin structure, spanning 8-bp long, which also existed in other Coptotermes species. No typical microsatellite-like regions were found in the AT-rich region of C. suzhouensis, which are commonly found in other insect species, but absent from all reported Rhinotermitidae species identified thus far (Beier et al. 2017).

In termites, species descriptions have historically relied upon morphological characters of the soldiers and alates, which has contributed to synonyms within Coptotermes (more than 40 junior synonyms) (Krishna et al. 2013). A recent analysis using molecular data has provided some insight into Coptotermes phylogenetics and radiation (Bourguignon et al. 2017). For example, Coptotermes havilandi (Isoptera: Rhinotermitidae) and Coptotermes gestroi (Isoptera: Rhinotermitidae), Coptotermes elisae (Desneux) (Isoptera: Rhinotermitidae), and Coptotermes curvignathus (Isoptera: Rhinotermitidae), as well as C. gestroi and Coptotermes vastator (Isoptera: Rhinotermitidae), were each revealed to be synonymous (Kirton and Brown 2003, Wang 2004, Kirton et al. 2005, Bengkeok et al. 2007, Lee et al. 2015) . In this study, we described the complete mitogenome sequence of C. suzhouensis and compared it to the available mitogenomes of other termites using phylogenetic analyses. Our results suggest that C. suzhouensis and C. formosanus are likely to be synonyms.

Click here for additional data file.

Click here for additional data file.