Evolutionary trend of the small subunit nuclear ribosomal DNAs (SSU rDNAs) of wood-feeding Cryptocercus cockroaches (Blattaria: Cryptocercus)

Abstract

The objective of this study is to elucidate the evolutionary rate of Cryptocercus SSU rDNAs based on the relative rate test (RRT) and comparative analysis of the primary sequences. In the RRTs, the evolutionary rate differences (dKs) of the Cryptocercus SSU rDNAs were not statistically significant (0.107 <P<0.996) in all cases, although an Asian Cryptocercus appeared to evolve more slowly than a North American Cryptocercus as indicated by the negative signs of dKs. Contrary to the RRT, the comparative analysis of SSU rDNA sequences revealed some interesting aspects in the evolutionary rates of Cryptocercus SSU rDNAs. In comparison with those of closely (Mastotermes darwiniensis and Tenodera aridifolia) or more distantly related taxa (Locusta migratoria), the SSU rDNAs of the genus Cryptocercus showed that insertions are progressing in a variable region (HV2). The SSU rDNA of Calymmaderus punctulatus species complex (Mountain Lake Biological Station, VA), a North American species, was unusually more elongated in the HV2 than those of the reference taxa or the other congeneric species.

Keywords:

• 18S rDNA
• Cryptocercus
• relative rate test
• hypervariable region

Introduction

Wood-feeding Cryptocercus cockroaches, which live in high mountainous forests, have been well known because of the peculiar characteristics of their life history and geographic distribution. Cryptocercus adults live with their offspring over an extended period of 3 years and provide their young with proctodeal feeding (Seelinger & Seelinger 1983; Nalepa 1984; Park & Choe 2003a, 2003b).

Cryptocercus species have been found only in a few forests of the Palearctic and Nearctic regions (Scudder 1862; Bey-Bienko 1950; Nalepa et al. 1997, 2001; Burnside et al. 1999; Grandcolas 1999a, 1999b, 2000; Grandcolas et al. 2001; Park et al. 2004). In Nearctic regions, they show a disjunct distribution in the USA; Cryptocercus clevelandi is distributed in the Coast and Cascade Mountains of northern California and southwestern Oregon (Nalepa et al. 1997), while Cryptocercus punctulatus species complex, also described as Cryptocercus darwini, Cryptocercus garciai, Cryptocercus punctulatus and Cryptocercus wrighti by Burnside et al. (1999), is reported to occur along the Appalachian and Allegheny Mountains (Burnside et al. 1999; Nalepa et al. 2002). In Palearctic regions, Cryptocercus has been reported in only a few forests of Southwest China and Northeast Asia including Korea, Manchuria and Siberia (Bey-Bienko 1950; Grandcolas 2000; Grandcolas et al. 2001; Nalepa et al. 2001; Park et al. 2004).

Since they are wingless and live only in rotten logs, their migration appears to be restricted to local boundary zones. Thus, the causes and dynamics that shape their present-day distribution have been the subject of vigorous discussion over the last decade, associated with their taxonomic status as well as phylogeography (Burnside et al. 1999; Nalepa & Bandi 1999; Grandcolas 1999a, 1999b; Clark et al. 2001; Hossain & Kambhampati 2001; Nalepa et al. 2002; Park et al. 2004).

In addition, to assess inter- and/or intra-specific genetic variation, or to better understand the spatiotemporal history of the geographic isolates, studies using mitochondrial DNA markers have been conducted extensively for Cryptocercus species, especially for the Nearctic cockroaches (Burnside et al. 1999; Clark et al. 2001; Hossain & Kambhampati 2001). In karyotype analysis, four variants (male 2n=37, 39, 43, 45) were found in C. punctulatus species complex in the southern Appalachian Mountains of eastern North America, and none of the variants occurs in the same locality, displaying different distribution range (Burnside et al. 1999; Nalepa et al. 2002).

According to two evolutionary scenarios of the variants suggested by Nalepa et al. (2002), the karyotype variants might have evolved independently from an ancestral population without having sequential events to one another (parallel scenario), or they might have been decreased sequentially from the ancestral karyotype (sequential scenario). As in the recent studies using karyotype and mtDNA markers (Everaerts et al. 2008), further molecular survey for these geographic isolates may provide clues to better understanding of the evolutionary trend of the genus Cryptocercus related to the molecular sequence evolution.

In this paper, we used nearly a complete sequence of small subunit nuclear ribosomal DNA (SSU rDNA or 18S rDNA) to infer the evolutionary trend of Cryptocercus. The 18S rDNA sequence has been generally known to include both very highly conserved and variable regions in their sequences. In the conserved regions high level of sequence identity is observed even among distantly related taxa (Turbeville et al. 1991; Whiting et al. 1997), while nucleotide insertions and deletions (indels) and/or nucleotide substitution are often found in the relatively fast-evolving variable regions even among closely related taxa (Hancock 1995).

These features make the utility of SSU rDNA as one of the most widely used molecular markers in multidisciplinary studies of the evolutionary trend in the nucleotide sequence and phylogenetic relationships among both closely and more distantly related taxa (Vogler et al. 1997; Whiting et al. 1997). In this study, we analysed primary structures of SSU rDNA sequences from Northeast Asian and North American Cryptocercus. Based on the comparative analyses of the sequences, we discussed the evolutionary trends between Cryptocercus of Northeast Asia and eastern North America.

Materials and methods

Sample collection

This study includes C. punctulatus species complex from the USA and two species (Cryptocercus kyebangensis and Cryptocercus relictus) from Northeast Asia. We collected C. kyebangensis from Mt. Yongmun in South Korea, C. relictus from Vladivostok in eastern Russia and C. punctulatus species complex from the Mountain Lake Biological Station (MLBS), Virginia in eastern North America. Although SSU rDNA has been known to be highly conserved at within-species level, the target gene fragments for two individuals from each locality were sequenced in order to detect potential sequence variation within the species.

DNA extraction, amplification, purification and sequencing

Genomic DNA was extracted from legs of the samples preserved in 100% ethanol using AccuPrep Genomic DNA Extraction Kit (Bioneer). Nearly complete sequences of 18S gene were amplified using polymerase chain reaction (PCR). Using a Gene®Amp PCR system 9700 (Perkin-Elmer), the polymerase chain reaction was carried out in 20 µl volumes, using 2 µl of genomic DNA as template. The PCR reactions included 10 pmol of each primer, dNTPs at 250 µM, 2.5 mM MgCl₂ and 1 unit of Taq Polymerase (Promega). Primers used for the amplification of SSU rDNA sequences are as follows: (18e) 5′-CTGGTTGATCCTGCCAGT-3′ (Hillis & Dixon 1991) and (18P) 5′-TAATGATCCTTCCGCAGGTTCACCT-3′ (Halanych et al. 1998). The temperature profile for amplifying target gene sequence was as follows: a denaturation at 94°C for 3 min and 30 cycles of 94°C for 1 min, 58°C for 90s and 72°C for 2 min followed by a final extension of 72°C for 10 min. PCR products were electrophoresed on a 1% agarose gel and purified using Cleanmix DNA Purification Kit (TA 200 CLN). Purified products were used as templates for sequencing. Sequencing reactions were performed using the ABI PRISM Dye Terminator Cycle Sequencing Ready Reaction Kit (Perkin-Elmer), according to the manufacturer's instructions, and a Gene®Amp PCR system 9700 (Perkin-Elmer). Sequences were determined by automatic sequencing on a 3730 DNA Sequencer (ABI) through 8% polyacrylamide gel. Both complementary strands were sequenced for all individuals. All of the sequences obtained in the current study are deposited at the Genbank under the accession numbers of JX091733 for C. relictus, JX091734 for C. kyebnagnesis and JX091735 for C. punctulatus species complex (MLBS).

Data analysis

For comparative sequence analysis, sequences of SSU rDNAs were aligned using Clustal X program package (Thompson et al. 1997) (Supplementary 1). As references for comparative analysis of the aligned Cryptocercus sequences, we used previously published sequences of Mastotermes darwiniensis (AF220568; Genbank accession numbers) (Isoptera) (Lo et al. 2000) and Tenodera aridifolia (AF423805) (Mantodea) (Whiting 2002) closely related to the genus Cryptocercus and that of Locusta migratoria (AF370793) (Orthoptera) (Giribet et al. 2001) distantly related to the genus. We also used a previously published sequence of C. punctulatus species complex (AF220571) (Lo et al. 2000) to compare with C. punctulatus species complex (MLBS).

The relative rate test (RRT) was conducted to estimate relative rates of nucleotide substitution of the SSU rDNAs sequences using the RRTree software (Relative Rate Test with Tree; Version 1.1) (Robinson & Huchon 2000). For the RRT, the SSU rDNA target sequences of five taxa were included as follows: C. relictus, C. kyebangensis, C. punctulatus species complex (AF220571), C. punctulatus species complex (MLBS) and Tenodera aridifolia (outgroup). The RRT has often been employed to determine whether the observed substitution rates are constant by comparing two sequences relative to the third as an outgroup (Sarich & Wilson 1973; Wu & Li 1985).

According to Dictyopteran phylogeny, cockroaches (Blattaria) and termites (Isoptera) are a sister and mantids (Mantodea) are positioned basal to the two former groups (Figure 1). Thus, for the RRT of the genus Cryptocercus, we used T. aridifolia (AF423805) (Mantodea) (Whiting 2002) as an outgroup.

Figure 1.  Phylogeny of the taxa used in this study. The phylogenetic tree was a modification of the tree that Lo et al. (2000) inferred from 18S r DNA and COII. Each species of Isoptera, Mantodea and Orthoptera was used as references for comparative analysis of nucleotide sites in the aligned sequences of the genus Cryptocercus. T. aridifolia was used as an outgroup in the relative rate test.

In the RRT, substitution rates were corrected by Kimura 2-parameter model (Kimura 1980) implemented in the RRTree software. SPSS (ver. 17) package was used for general statistical analyses of base composition frequencies.

Results

Sequence alignment

We obtained nearly complete SSU rDNA sequences for each of the two Cryptocercus individuals from South Korea, eastern Russia and eastern North America, respectively. There was no individual variation in SSU rDNA sequences for Crypotocercus collected from these localities. Thus, one SSU rDNA genotype representing each collection locality was used for further analyses. The aligned SSU rDNA sequences with reference taxa were 1772 nucleotides in length (Supplementary 1).

There were three relatively well-conserved regions (at the sites 1–140, 157–642 and 771–1772) across the aligned sequences, whereas two variable regions with many indels were found at the sites 141–156 (HV1) and 643–770 (HV2) covering 144 sites in their total length (Figure 2 and Supplementary 1).

Nucleotide composition

Base composition was homogenous across the taxa (Chi-square test; 0.77<P<1) (Table 1). In reference sequences, the chi-square values of SSU rDNAs were higher in the conserved regions (2.03) than in the variable regions (0.74), while in Cryptocercus, the values were higher in the variable regions (5.67) than in the conserved regions (0.07).

Table 1. Chi-square test of homogeneity of base frequencies across taxa.

		Statistic values
Taxa	Sites	χ ^2	Df	p
Reference and Cryptocercus	Total	4.621652	18	0.99933706
	Conserved	3.334442	18	0.99993808
	Variable (HV1 and HV2)	7.802381	18	0.98143524
Reference	Total	2.005452	6	0.91919653
	Conserved	2.030017	6	0.91691732
	Variable	0.741817	6	0.99354376
Cryptocercus	Total	0.636435	9	0.99991472
	Conserved	0.068465	9	1
	Variable (HV1 and HV2)	5.670053	9	0.77243242

In the entire regions of SSU rDNAs, the GC contents were a little higher than that of AT in both reference and Cryptocercus sequences (Mann–Whitney U; P<0.05) (Table 2). The GC contents were also a little higher than that of AT in both conserved regions and variable regions (Mann–Whitney U; P<0.05 in both cases). Average GC/AT ratio was significantly higher in the variable regions (2.79±0.34) than in the conserved region (1.11±0.03) (Mann–Whitney U; P<0.05).

Table 2. Base composition of the whole in the regions of SSU rRNA sequences.

Taxon	C	G	A	T	GC content	AT content	GC/AT ratio	No. sites
L. migratoria	0.2394	0.2824	0.2354	0.2428	0.5218	0.4782	1.0911	1768
T. aridifolia	0.2493	0.2831	0.2332	0.2344	0.5324	0.4676	1.1385	1745
M. darwinensis	0.2518	0.2896	0.2337	0.2249	0.5414	0.4586	1.1805	1716
Submean	0.2468	0.2850	0.2341	0.2340	0.5319	0.4682	1.1367	1743
	(±0.00654)	(±0.00399)	(±0.00116)	(±0.00893)	(±0.00981)	(±0.00981)	(±0.04469)	(±26.1)
C. relictus	0.2536	0.2875	0.2314	0.2275	0.5412	0.4588	1.1794	1798
C. kyebangensis	0.2521	0.2867	0.2331	0.2281	0.5388	0.4612	1.1681	1793
C. punctulatus species complex (AF220571)	0.2588	0.2854	0.2282	0.2276	0.5442	0.4558	1.1939	1731
C. punctulatus species complex (MLBS)	0.2605	0.2873	0.2256	0.2267	0.5478	0.4522	1.2113	1831
Submean	0.2563	0.2867	0.2296	0.2275	0.5430	0.4570	1.1882	1788.3
	(±0.00404)	(±0.00096)	(±0.00336)	(±0.00061)	(±0.0039)	(±0.0039)	(±0.0187)	(±41.7)
Total mean	0.2142	0.2391	0.1919	0.1897	0.5382	0.4618	1.1661	1497.2
	(±0.10302)	(±0.11666)	(±0.09238)	(±0.09261)	(±0.00866)	(±0.00866)	(±0.04)	(±713.7)
Note: The values in parentheses indicate standard deviations.

Variable regions

In the HV1, the SSU rDNAs of Cryptocercus showed a high degree of length variability by insertions and deletions (Figure 2). In comparison with those of reference taxa, SSU rDNAs of Cryptocercus were inserted by an average of 7.2 nucleotides (±3.0, n=12) in the HV1 and an average of 31.8 nucleotides (±20.2, n=12) in the HV2, respectively. One of the most interesting results is that the HV2 region of C. punctulatus species complex (MBLS) was exceptionally more elongated than that of any other congeneric members or reference taxa. It was 35.0 (±2.6, n=3, range=33–38) nucleotides longer than that of the other Cryptocercus and 58.0 (±14.5, n=3, range=44–73) longer than that of the references.

Figure 2.  HV1 (at the sites 141–156) and HV2 (at the sites 643–770) in the SSN rRNA sequences of the genus Cryptocercus aligned with those of references (L. migratoria, T. aridifolia and M. darwinensis). An unusual expansion by insertion occurs in the HV2 of the SSU rRNA of C. punctulatus species complex (MLBS). Numbers above the alignment indicate the positions of the nucleotide sites in the full sequence alignment (see Supplementary 1).

In all HV1 and HV2 (Figure 2), with the exception of the difference in gap number, there was no nucleotide substitution between Asian C. relictus and C. kyebangensis and only one substitution at the site 759 (transition; A↔T) between North American C. punctulatus species complex (AF220571) and C. punctulatus species complex (MLBS), while a total of 26 substitutions (6 sites in the HV1 and 20 sites in the HV2) were detected between the Asian and North American Cryptocercus (shown in boldface nucleotides in Figure 2).

Conserved regions

In three relatively well-conserved regions (1–140, 157–642 and 771–1772 sites), excluding the variable regions, four substitutions were found (three transitions and one transversion) between two Asian Cryptocercus species, but no substitution between two individuals of North American C. punctulatus species complex. Substitutions occurred in all 13 sites (transition at five sites and transversion at eight sites) between Asian and North American Cryptocercus.

Nucleotide sequences at 17 variable sites detected in the conserved regions can distinguish the genus Cryptocercus from the reference taxa. Nucleotide sequences at these variable sites were highly conserved in all reference taxa (no substitutions detected), but substitutions occurred only in the genus Cryptocercus (Figure 3). Seven of 17 sites were substitutions specific to all species of the genus Cryptocercus (as indicated in the boxes in Figure 3). Substitutions between Asian Cryptocercus (C. relictus and C. kyebangensis) occurred in only two sites (at sites 3 and 786), while there was no substitution between the two individuals of North American C. punctulatus species complex. Substitutions between Asian and North American Cryptocercus species occurred in an average of 14 sites (±1.4, n=2).

Figure 3.  The selected 17 sites whose nucleotides were well conserved in the SSU rDNAs of references, but substituted in those of the genus Cryptocercus. Bold nucleotides in the sites of the genus Cryptocercus SSN rRNAs indicate that they were different from those in the corresponding sites of the references SSU rRNAs. The boxed sites in the genus Cryptocercus indicate that at the sites the nucleotides were different from those of references, but same within the genus.

C. relictus was substituted in 10 of the 17 sites. In case of C. kyebangnesis substitution occurred in 12 sites. Thus, the substitution between two Asian Cryptocercus and the references occurred on average in 11 sites (±1.4, n=2), while the substitution between North American C. punctulatus species complex and the references occurred in total 13 sites.

Relative rate test

For the RRT, the data-set is composed of the entire sequences (1772 nucleotide sites) including both conserved and variable regions (HV1 and HV2). However, total nucleotide sites used in the RRT were 1675 of the 1772 sites, probably since all sites with gaps were excluded in the RRT using the RRTree software (Robinson & Huchon 2000).

Results from the RRT showed that the differences (dKs=K _ar −K _ap , K _ar −K _am , K _ak −K _ap and K _ak −K _am ) of the substitution rate (K1 and K2; substitution rates between Asian Cryptocercus and T. aridifolia and between Northeast American Cryptocercus and T. aridifolia, respectively) were all negative signs, indicating a slower substitution rate in Asian than in North American Cryptocercus (Table 3). The substitution rate of C. relictus was higher than that of C. kyebangensis since the difference of substitution rates between each of two Asian Cryptocercus and T. aridifolia (dK=K _ar −K _ak ) has a positive value. The difference of substitution rates (dK=K _ap −K _am ) between each of two individuals of C. punctulatus species complex and T. aridifolia was a negative sign, implying a slower substitution rate in C. punctulatus species complex (AF220571) than C. punctulatus species complex (MLBS).

Discussion

The Cryptocercus SSU rNDAs were not different from the general features that have ever been known from many metazoan animals, in that they have both slowly evolving conserved regions and fast evolving variable regions with relatively higher GC content (Tautz et al. 1988; Hancock 1995; Crease & Colbourne 1998; Choe et al. 1999). The comparative analysis of SSU rDNA sequences provided insights into evolutionary trends of the genus Cryptocercus. Length variations in aligned SSU rDNA sequences became much more prominent in the hypervariable regions, especially in the HV2. The unexpected expansion in the region occurred most remarkably in C. punctulatus species complex (MLBS) but not in the other Cryptocercus species. More insertions were found in the HV2 (an average of 35 nucleotides) of C. punctulatus species complex (MLBS) than that of Cryptocercus punctulatus species complex (AF220571). The high GC content in this variable region indicates that length expansion was caused mainly by insertions of G and C. With the exclusion of the insertion part, only a single transitional change was detected between HV2 of C. punctulatus species complex (MLBS) and C. punctulatus species complex (AF220571) (Figure 2: A↔T in the site 759 in HV2).

Comparative sequence analysis of the conserved regions showed a peculiar substitution patterns in a total of the 17 sites (Figure 3) where nucleotide sequences are the same all across reference taxa, but some substitutions were found between the Cryptocercus species. Among the 17 sites, seven were transitional substitutions specific to all Cryptocercus species, leading us to conjecture that these sites may be derived characters to distinguish the genus Cryptocercus from the reference taxa. In the conserved 17 sites, nucleotides were less substituted in Asian than in North American Cryptocercus, relative to reference sequences.

In the RRT, the substitution rates did not represent statistically significant differences in all cases, although the Asian Cryptocercus tend to show slower substitution rate than the North American Cryptocercus as indicated by the negative signs of dKs (Table 3). The homogeneity test of base frequencies also showed that the sequences of SSU rDNAs have evolved with the same rate of nucleotide substitution in their evolutionary history (homogeneity assumption, 0.77 <P <1 in Table 1). According to the molecular clock hypothesis (Zuckerkandl & Pauling 1965), substitutions occur at the same rate for all across the sites in closely related species. As in the RRT as well as the homogeneity test, thus, the substitution rate of Cryptocercus SSU rDNAs might not show statistically significant difference, probably because the sequences compared were too closely related to each other.

Table 3. RRT of the aligned 1675 sites excluding gaps in the entire sites of SSU rDNA sequences.

Outgroup	Lineage1	Lineage2	Number of sites	K1		K2		dK	SD_dK	Ratio dK/SD	P
T. aridifolia	C. relictus	C. kyebangensis	1675	K ar	0.03178	K ak	0.03177	0.00001	0.00164	0.00533	0.996
T. aridifolia	C. relictus	C. punctulatus species complex (AF220571)	1675	K ar	0.03178	K ap	0.03618	−0.00440	0.00346	−1.27022	0.204
T. aridifolia	C. relictus	C. punctulatus species complex (MLBS)	1675	K ar	0.03178	K am	0.03619	−0.00441	0.00352	−1.25347	0.210
T. aridifolia	C. kyebangensis	C. punctulatus species complex (AF220571)	1675	K ak	0.03177	K ap	0.03618	−0.00441	0.00362	−1.21654	0.224
T. aridifolia	C. kyebangensis	C. punctulatus species complex (MLBS)	1675	K ak	0.03177	K am	0.03619	−0.00442	0.00368	−1.20209	0.229
T. aridifolia	C. punctulatus species complex (AF220571)	C. punctulatus species complex (MLBS)	1675	K ap	0.03618	K am	0.03619	−0.00001	0.00061	−0.01785	0.986
Note: The values of K1 and K2 indicate the substitution rates between outgroup (T. aridifolia) and each species of Lineage 1 and Lineage 2, respectively.
The Kij in K1 and K2 columns indicates each substitution rate between species i (=a) and j (= r, k, p, or m)
In the Kij, the subscripts are a for T. aridifolia, r for C. relictus, k for C. kyebangensis, p for C. punctulatus species complex (AF220571) and m for C. punctulatus species complex (MLBS).
The dK indicates difference of K1 and K2 (dK=K1−K2).
Transitional and transversitional substitutions were corrected by Kimura 2-parameter.

In conclusion, nucleotides of the conserved sites were less substituted in Asian than in North American Cryptocercus, relative to reference sequences, perhaps suggesting that the Asian Cryptocercus might be more slowly evolved than the North-east American Cryptocercus. This interpretation is well in congruence with earlier arguments on the Dictyopteran phylogeny that Asian C. relictus is positioned basal to the North American Cryptocercus species (Clark et al. 2001; Lo et al. 2003). However, further examination based on different genetic markers and more extensive samplings is required to elucidate the evolutionary rate of the genus Cryptocercus.

Acknowledgements

We would like to thank Dr Christine A. Nalepa (Department of Entomology, North Carolina State University) very much for providing the specimen of North American Cryptocercus. This work was carried out with the support of ‘Forest Science & Technology Projects (Project No. S211012L030110)’ provided by Korea Forest Service.