Taxon sampling

Samples were collected from rotting logs at 28 localities in China (Fig. 1), from three main regions: Hengduan Mountains, Qin-Daba Mountains and Manchuria (Table 1). Collections were limited to single-day-trip trekking in one valley, meaning that deep valleys were difficult to reach and usually not sampled. Specimens were preserved in analytical pure ethanol and stored at −80°C until processing. All specimens are deposited in the institute of Entomology, College of Plant Protection, South-west University, Chongqing, China.

Fig. 1. 

Map of sampled localities for Cryptocercus. Numbers for sampling localities match those in Table 1 (red circles denote samples collected by our team, blue circles denote data obtained from GenBank). The map was generated by  www.simplemappr.net using GPS coordinates.

Table 1. 

Samples used in species delimitation: number of location, sample collection localities, abbreviation of location, and GenBank accession numbers

All localities are in China unless otherwise indicated. The letter ‘n’ after the voucher abbreviation means the sample is a nymph and the letter ‘f’ means the sample is female. n.a., not available

Morphological analyses

The terminology for morphological structures follows McKittrick (1964), Li et al. (2013) and Wang et al. (2015). We identified the species based on five standard characters of female genitalia (anterior margin of Tergite VII, shape of the spermatheca, shape of the basivalvulae, shape and coloration of the laterosternal shelf, length of the paraprocts), which are stable and can be used to distinguish species to a certain extent (Grandcolas 2000; Grandcolas et al. 2001, 2005; Aldrich et al. 2004; Wang et al. 2015). We use abbreviations to indicate the following structures: paraprocts (pp.), first to third valve (v.I, v.II, v.III), paratergites (pt.), anterior arch (a.a.), basivalvula (bsv.), laterosternite IX (ltst.IX), laterosternal shelf (ltst.sh), vestibular sclerite (vst.s.), intersternal fold (inst.f.), central apodeme (c.a.) and spermatheca (sp.). Preparation of genitalia followed the standard procedure, using 10% NaOH to clear internal tissues, then all genitalia were mounted on microscope slides using a Motic K400 stereomicroscope. Photographs of the female genitalia were made using a Leica M205A microscope with a Leica DFC camera. Adults were photographed with a digital camera (Canon 50D) with the aid of the Helicon Focus software.

Karyotype analysis

For karyotype analysis, three individuals from three different logs within each population were examined. Mitotic chromosomes from the testes of males were examined following Luykx (1983). Karyotypes are reported as the diploid complement.

DNA extraction and sequencing

Genomic DNA was extracted from legs by using the TIANamp Genomic DNA Kit (DP304, TIANGEN). Primers are provided in  Table S1 (IS17003_AC.pdf), available as supplementary material to this paper. All the reactions were carried out in volumes of 25 μL, containing 14.25 μL of ultrapure water, 2.5 μL of 10 × buffer (Mg²⁺ free), 2 μL of MgCl₂ (25 mM), 2 μL of dNTP mixture, 1 μL of each primer, 0.25 μL of Taq polymerase, and 2 μL of DNA template. The following steps were performed on a programmable thermal cycler. The amplification protocol setting used were 94°C for 5 min; followed by 35 cycles at 94°C for 45s, 48°C for 45s, and 72°C for 45s; and final extension at 72°C for 10 min. PCR reactions were separated by electrophoresis on a 1% agarose gel. All DNA purification and sequencing was carried out by BGI Tech (Beijing, China) using the aforementioned primers. All sequences were deposited in GenBank, at the National Center for Biotechnology Information (Table 1).

Sequence analysis

A total of 135 cytochrome c oxidase subunit I (COI) sequences were analysed, including: 86 sequences from this study; five sequences representing four species from North America and Korea downloaded from GenBank; and 44 sequences representing termite outgroups (Table S2). To provide a comparison with the results from COI, we also analysed two sequences of 28S rRNA for each species (74 total sequences) (Table 1). COI sequences were aligned using MUSCLE 3.8 (Edgar 2004) and adjusted visually and manually corrected after translation into amino acid sequences; 28S rRNA sequences were inspected visually and manually. Intraspecific and interspecific genetic divergence values were quantified based on the Kimura 2-parameter (K2P) distance model (Kimura 1980), using MEGA 6.0.6 (Tamura et al. 2013), with 1000 bootstrap replicates.

To explore phylogenetic relationship among these closely related species, phylogenetic trees were constructed on two different datasets (single COI dataset and the concatenated dataset) using Maximum Likelihood (ML) and Bayesian Inference (BI). Bayesian analysis was implemented in MrBayes 3.2 (Ronquist et al. 2012), assigning site-specific models. The COI dataset was divided into three partitions by codon position (pos1–3), and PartitionFinder 1.1.1 (Lanfear et al. 2012) was used to determine the best-fitting models for each partition. The best-fitting models were as follows: COI_pos1, TVM+G; COI_pos2, TrNef+G; COI_pos3, TrN+I+G; 28S rRNA, TrN+I+G. Two independent sets of Markov chains were run, each with one cold and three heated chains for 1 × 10⁷ generations, and every 1000th generation was sampled. Convergence was inferred when a standard deviation of split frequencies <0.01 was completed. Sump and sumt burninfrac was set to 25% and contype was set to allcompat. Maximum likelihood was implemented in RAxML 7.3.0 (Stamatakis et al. 2008); as the software does not allow different substitution models for different partitions we used the GTRGAMMA model for analyses. Bootstrap values were performed for 1000 replicates.

All COI sequences were analysed using ABGD and GMYC. ABGD analysis was performed using a web interface ( http://wwwabi.snv.jussieu.fr/public/abgd/, accessed 4 May 2017); the default parameters were used except for the relative gap width, which was set as 1.0. The GMYC method requires a fully resolved ultrametric tree as input for the analysis. We constructed a Bayesian inference tree in BEAST 1.8.1 (Drummond and Rambaut 2007) using the best-fitting models from PartitionFinder. Rate variation was modelled among branches using an uncorrelated log-normal relaxed clock model with the mean clock rate fixed to 1.0. A UPGMA starting tree was used in the analysis, and the constant size coalescent was used as a prior on divergence times. We performed two replicate MCMC runs, and sampled every 5000 steps over a total of 50 million generations. A maximum clade credibility tree was obtained using Tree Annotator within the BEAST software package with a burn-in of 1000 trees. We then applied the single-threshold GMYC method to the ultrametric gene tree generated by BEAST using the SPLITS package (Ezard et al. 2009) in R software (R Core Team 2013). The groups delimited were compared with a one-species null model using a likelihood ratio test.

Karyotypes and taxonomy of novel Cryptocercus species

Through examination of the morphology and chromosome numbers of cockroaches from previously unsampled areas, we identified nine new candidate Cryptocercus species: C. changbaiensis, sp. nov.; C. weixiensis, sp. nov.; C. luanchuanensis, sp. nov. (2n = 39) (Fig. 2H); C. chengkouensis, sp. nov. (2n = 29) (Fig. 2A); C. wolongensis, sp. nov. (2n = 43) (Fig. 2B); C. banshanmenensis, sp. nov. (2n = 29) (Fig. 2C); C. tianbaensis, sp. nov. (2n = 35) (Fig. 2D); C. laojunensis, sp. nov. (2n = 21) (Fig. 2E); C. pudacuoensis, sp. nov. (2n = 21) (Fig. 2F). All species of Chinese Cryptocercus examined in this study were highly similar in external morphology, including female genitalia (see Table 2, e.g. C. luanchuanensis, sp. nov. and C. neixiangensis). All candidate new species are named after their type localities.

Fig. 2. 

Mitotic chromosomes of males of Cryptocercus species. (A) C. chengkouensis, sp. nov. (XLC) (2n = 29); (B) C. wolongensis, sp. nov. (DS) (2n = 43); (C) C. banshanmenensis, sp. nov. (BSM) (2n = 29); (D) C. tianbaensis, sp. nov. (MJG) (2n = 35); (E) C. laojunensis, sp. nov. (JSJLT) (2n = 21); (F) C. pudacuoensis, sp. nov. (PDC) (2n = 21); (G) C. neixiangensis (BTM) (2n = 37); (H) C. luanchuanensis, sp. nov. (LYW) (2n = 39).

Table 2. 

Comparison of female genitalia of 21 Chinese Cryptocercus species

COI and 28S rRNA sequence variation across Asian taxa

We next obtained molecular sequence data from mitochondrial COI and 28S rRNA from multiple Asian taxa. The sequenced lengths of COI and 28S rRNA, excluding primers, were ∼658 bp and 635 bp respectively. All new sequences have been deposited in GenBank with accession numbers KY940572 to KY940657 for COI (accession numbers for 28S in Table 1). Among aligned COI sequences 232 sites were variable, of which 217 were parsimony informative. Among aligned 28S rRNA sequences, 106 sites were variable, of which 103 were parsimony informative.

Intraspecific COI genetic divergence (K2P) ranged from 0.00 to 0.61%, with an average of 0.14%. The greatest intraspecific COI genetic divergence (0.61%) occurred in C. darwini. All intraspecific COI genetic divergence values were less than 0.5%, and many intraspecific COI genetic divergence values were 0.00% (Table S3). Interspecific COI genetic divergence ranged from 2.18% (C. neixiangensis and C. luanchuanensis, sp. nov.) to 20.36% (C. chengkouensis, sp. nov. and C. darwini), with an average of 11.81% (Table S3). The 28S sequences of Cryptocercus had low levels of genetic divergence of 96.27% and lacked sufficient variation to resolve some species (Table S4). 28 species were identified through the use of the barcoding gap, of which nine were new (Fig. 3B).

Fig. 3. 

Maximum-likelihood tree based on COI genes (BPP: Bayesian Posterior Probabilities; MLB: Maximum Likelihood Bootstrap Values). A square (■) at a node indicates that both BPP and MLB >95; a circle (●) indicates that only MLB >95. Numbers after the species name are chromosome numbers, numbers in red are results from this study. Coloured bars indicate different species delimitation by different methods: A, morphology; B, DNA barcoding gap; C, GMYC model; D, ABGD result. We selected 86 samples of 28 Cryptocercus species. The topology shown was very similar to that produced by the BI tree, with some differences (see Fig. S1).

GMYC, ABGD, and other phylogenetic analyses

The likelihoods of the null and GMYC models from COI analysis were 239.49 and 268.42 respectively. The GMYC was an improvement over the null model, and was clustered into 72 (confidence interval: 72–73) entities (likelihood ratio = 57.87) including 28 Cryptocercus and 44 termite species. ABGD based on the same dataset also detected 72 species with a range of prior intraspecific values (P = 0.0046–0.0215) (Fig. 3C).

Phylogenetic analyses performed using ML and BI recovered similar tree topologies in COI analysis with some differences in branching patterns at deeper nodes. Support values from ML analysis were much higher than those from BI analysis (Figs 3, S1). Representatives from candidate species (described primarily on the basis of karyotype number) were found to form monophyletic groups. Analyses based on concatenated COI and 28S sequences revealed similar topologies to those inferred from COI only, with some minor differences.

Asian and American lineages (α and β) formed monophyletic groups, as recovered in BI and ML analyses for COI and combined datasets (Figs 3, S1, S2). For COI analysis (Fig. 3), all species from Qin-Dada Mountains were clustered together with strong support (Clade γ). Species from Manchuria and the Hengduan Mountains were clustered together (Clade δ) with one exception (C. pudacuoensis, sp. nov. (PDC)), although this grouping was not well supported.

Species delimitation

The nine candidate species we identified, primarily on the basis of chromosome numbers, are in addition to the 15 identified previously by Grandcolas (2000), Grandcolas et al. (2005), Wang et al. (2015), and Che et al. (2016). Molecular species delimitation analysis corroborated the existence of each of these species. Our results indicate that DNA-based species-delimitation methods using COI perform well for these subsocial and xylophagous cockroaches, and are likely to be of utility given the lack of defining morphological characters among males, females and nymphs of these organisms. Our study is the first attempt to delimit Cryptocercus species on a large scale, including males, females and nymphs.

As the genus Cryptocercus comprises mainly alpine cockroaches, we were able to obtain only a small number of samples from each locality. This may have contributed to the limited genetic variability within each species. Barcode gap (COI) analysis showed that the maximum intraspecific distance (0.61%) was distinctly lower than the interspecific distance (2.18–20.36%), even for those species (C. neixiangensis and C. luanchuanensis, sp. nov.) with the lowest interspecific divergence (2.18%). Hebert et al. (2004) proposed that the genetic divergence cutoff for species identification should be at least 10 times greater than within species. Our study (average intraspecific distance: 0.14%, average interspecific distance: 11.81%) showed that there is a distinct gap between intraspecific and interspecific distances. DNA-based analyses resolved most Cryptocercus samples from Asia to a putative species. GMYC, ABGD and chromosome numbers enabled their differentiation of the two geographically proximate species C. neixiangensis (BTM) (2n = 37) and C. luanchuanensis, sp. nov. (LYW) (2n = 39) (Fig. 2) (2.18% in COI, 1% in 28S), although they had highly similar morphologies. Some other species pairs were particularly well differentiated despite being geographically proximate to each other (e.g. C. primarius (DCP) and C. pingwuensis (BSG) – 4 km apart).

Uplifting of the Qinghai–Tibetan plateau caused rich diversity of Cryptocercus in the Hengduan Mountains

To date 22 Cryptocercus species have been described worldwide, of which 16 are from China (Bey-Bienko 1938; Grandcolas 2000; Grandcolas et al. 2005; Wang et al. 2015; Che et al. 2016). Another nine new Cryptocercus species, mainly from the Hengduan Mountains (six species), are described in the current study. There are now 25 Cryptocercus species distributed in the mountain forests from 702 m (Gaolingzi, Heilongjiang Province) to 3756 m (Shikaxueshan, Yunnan Province) in China, most of which are in the Hengduan Mountains.

The Hengduan Mountains comprise a large mountainous region in south-western China with an average elevation of more than 4000 m (Spicer et al. 2003). Fifteen Cryptocercus species were discovered here at an altitude of ∼3000 m. The Himalayan and Hengduan Mountains are known to encompass global biodiversity hotspots with high levels of plant and animal biodiversity and endemism (Myers et al. 2000; Wu et al. 2013; Liu et al. 2014). The region is thought to harbour at least 15 000 land plant species, of which more than 29% are endemic, and at least 1141 vertebrate species, of which 15% are endemic (Myers et al. 2000). The rapid and extensive uplift of the Himalaya–Hengduan Mountains is considered to be the major driving force in shaping such high species diversity (Liu et al. 2013; Wu et al. 2013; Wen et al. 2014). The Hengduan Mountains can be considered to be one of the hotspots of Cryptocercus diversity in China.

Since the early Miocene, extensive uplifting of the Qinghai–Tibetan plateau occurred over at least four major periods: 25–17, 15–13, 8–7, and 3.5–1.6 million years ago (Allégre et al. 1984; Spicer et al. 2003; Royden et al. 2008). The Cryptocercus lineages in the Hengduan Mountains began to diversify 20.7–11.7 million years ago, and the most recent diversification between C. arcuatus (SKXS) and C. habaensis (HBXS) occurred ∼5.82–1.58 million years ago (Che et al. 2016), which corresponds well with the major uplifts of the Qinghai–Tibetan plateau. This region is characterised by a series of parallel mountain ranges (Yin and Harrison 2000), but all members of Cryptocercus are wingless and lack the ability to migrate over long distances. The physical isolation caused by these mountain ranges is likely to have contributed to interspecific and intraspecific divergences. The localities of two close geographic species, C. habaensis (Fig. 14E) and C. meridianus (Fig. 14D), are found just 28 km apart, but they are separated by two snow mountains (Haba Snow Mountain, 5396 m; Yulong Snow Mountain, 5596 m) and one gorge (Tiger Leaping Gorge) (Fig. 14A). Tiger Leaping Gorge passes between the two Snow Mountains in a dramatic vertical drop of 3500 m, through which a rapid stream flows and xeric shrubland exists on both banks along the Jinsha River (Qikun Bai, pers. obs.). These natural barriers prevent gene flow among Cryptocercus species. Significant genetic differences (9.50% in COI) can be detected between these two species. Similarly, there is no vegetation on the peaks of ridges at altitudes of 4000–5000 m (Wang et al. 2004), which would prevent migration by Cryptocercus from one peak to another (Fig. 14B, 1C).

Fig. 14. 

Collection sites for Cryptocercus at Mt Habaxueshan and Mt Yulongxueshan: (A) the whole topographical map of Mt Habaxueshan and Mt Yulongxueshan, map data: Google, Copernicus; (B) photograph of Mt Yulongxueshan; (C) photograph of Mt Habaxueshan; (D) habitat of C. meridianus; (E) habitat of C. habaensis. All photographs by Qikun Bai.

Compared with the lineages from the Hengduan Mountains (2124–3756 m, average 2986 m), the lineages from the Qin-Daba Mountains (1470–2600 m, average 1814 m) and Manchuria (702–1199 m, average 935 m) occur at lower elevations. Relatively poor Cryptocercus diversity has been found in these two regions (Qin-Daba Mountains lineages: 7 species; Manchuria lineages: 2 species). Manchuria is a large north-eastern Asian geographic region of ∼130 × 10⁶ km² which mostly consists of vast plains (Elliott 2000). C. relictus has been collected by us at Gaolingzi, Shuangfenglinchang, and Mudanfeng as well as from Russia (Bey-Bienko 1935) and South Korea (Park et al. 2004). This speices therefore appears to be widely distributed in north-eastern Asia. The relatively simple geography of this region compared with the Hengduan Mountains is likely the primary cause of the low diversity in this area, because there are no major barriers to hinder the gene flow among different populations. Changbai Mountain, the highest in the Manchuria region, is an active volcano situated on the boundary between China and North Korea, and was formed within a short period between 2.77 and 0.31 million years ago (Wan 2012; Wei et al. 2007). It is therefore unlikely to have influenced the evolution of Cryptocercus in the way that the Hengduan Mountains have.