All adult samples described were collected between June and July 2016, at Yangliuping (29°41'2.54"N, 101°53'32.24"E), inside Yanzigou glacier, Mt. Gongga, Sichuan (Fig. 2, Fig. 8). Since Thitarodes species do not congregate at light traps and mating flights have not been observed for this particular species, collection was undertaken by thoroughly searching through habitat vegetation with flashlights. Adults can be found hanging at the edge of low vegetation, especially near where pupae molts can be seen. Collection for adults was also attempted at other glacial valleys along eastern Mt. Gongga from June to July between 2016 and 2018 but were unsuccessful. Pupae and larvae were collected at multiple glacial valleys along eastern Mt. Gongga (Fig. 8, see Suppl. material 2: Table S2 for names of seven valleys and their coordinates) as well as a habitat east of the Kangding-Moxi fault (Yajiageng, 29°53'53.12"N, 102° 0'37.87"E) from May to June, 2016 to 2018, with assistance from the local communities (see Acknowledgments). In each valley, at habitats where caterpillar fungus is harvested by local people, multiple 50 × 50 × 30 cm grassland plots were searched to discover ground-boring Thitarodes larvae and pupae. Grass coverage were placed back after sampling. Samples were preserved in 90% ethanol. Photographs of specimens were taken with a Nikon Coolpix 4500 digital camera. Dissections were performed after softening adult genitalia in heated 10% NaOH solution for 10–30 mins and transferring the genitalia on dissection slides with glycerol. Wing venation slides were obtained by softening the wings in 30% dish detergent for 30–60 mins, and gently brushing the scales off the wings. Slides were examined and photographed under a Carl Zeiss Stemi 2000-CS stereoscope system and Carl Zeiss SteREO Discovery V12 stereoscopic microscope system. Although previous authors have used the Comstock-Needham venation nomenclature to describe new species (e.g., Maczey et al. 2010; Zou et al. 2011), Wootton’s venation nomenclature (1979) was used in our description to reflect differences in wing homology.

Genomic DNA of 134 samples of Thitarodes were extracted from leg (adult) or thoracic (larvae, pupae) tissue with Qiagen DNeasy Blood and Tissue Kits (Qiagen Inc.). The COI region of 48 samples was amplified and sequenced with LCO1490 (Folmer et al. 1994) and Nancy (Bogdanowicz et al. 1993) primers, following the protocol for Lepidoptera COI sequencing outlined in Boyle et al. (2014) (with the exception that the annealing temperature was set at 55 °C, repeating for 35 cycles). Samples were sequenced at Thermo Fisher Scientific Inc, Shanghai. These sequences were aligned in Geneious Prime (using the MAFFT algorithm) together with nine other COI sequences of described Thitarodes species in China and available on GenBank (five of which are segments of mitogenome sequences). Outgroups were selected according to Grehan (2011, 2012). We used either COI sequences or the COI segment of the mitogenome sequence of the East Asian genus Napialus Chu & Wang, 1985, the Neotropical genus Phassus Walker, 1856, and the Australian genus Oxycanus Walker, 1856 (Hepialidae) (see Suppl. material 1: Table S1 for all sequences used). A Maximum Likelihood tree was inferred using the IQ-Tree algorithm (Nguyen et al. 2015) with a free rate variation model (Kalyaanamoorthy et al. 2017) and ultrafast bootstrap (UFBoot, Hoang et al. 2018) integrated 1000 times and SH-aLRT branch test (Guindon et al. 2010). All 134 extracted DNA samples were digested with EcoRI enzyme and library preparation followed the wet lab protocol of Etter et al. (2012). Samples were sent for 150bp pair-end sequencing on an Illumina HiSeq4000 platform located at 1gene Inc, Hangzhou. Illumina reads were demultiplexed and processed by STACKS (Rochette and Catchen 2017). For each sample, a locus was identified with a minimum depth coverage of three reads, and a 2bp distance was allowed between each locus. Loci having a distance less than 1bp across individuals were merged into consensus loci. The resulting matrix of all SNPs from all samples was converted into a presence/absence matrix, where missing reads (indicative of divergent enzymatic cut sites) were noted by 0 and otherwise 1, as in de Medeiros and Farrell 2018.The resulting matrix rows were arranged by hierarchical clustering (Ward 1963) to facilitate visualization.

A total of 541 available sequences of Thitarodes was used to build a summary tree. These sequences came from three sources: CO1 and mitogenome sequences of known Thitarodes species used in generating the Maximum Likelihood CO1 tree in this study (Suppl. material 1: Table S1), 26 Cytb sequences from known Thitarodes species as described in Zou et al. (2017), 226 published sequences of different segments extracted from caterpillar fungus sclerotia (Zhang et al. 2014; Quan et al. (2014a, b). The latter sequences were not associated with clear species names. CO1, Cytb and wg sequences from a single sample in Zhang et al. (2014) were concatenated, and CO1, Cytb and CO2 sequences from a single sample in Quan et al. (2014a, b) were concatenated. All sequences, after concatenation, were aligned in Geneious Prime (MAFFT algorithm). For mitogenome sequences, only the CO1, CO2 and Cytb segments were used in alignment. The Maximum Likelihood tree was inferred the using the same method CO1 sequence tree was inferred (see previous section). In the resulting phylogeny, each tip within 0.025 distance from each other were clustered and labelled as a single cluster to assess the degree of uneven sampling across the phylogeny.

COI sequences of 48 Thitarodes samples collected from six glacial valleys were successfully sequenced (Figs 8, 9A, Suppl. material 1: Table S1). The best model for the parsimony tree assessed by the BIC value determined by ModelFinder (Kalyaanamoorthy et al. 2017) is “K3Pu+F+I+G4”. All Thitarodes sequences form a highly supported monophyletic group, sister to non-Thitarodes outgroups. These samples cluster into three highly supported clades in the phylogeny when other available CO1 sequences of Thitarodes is taken into account (Fig. 9A): (1) All samples from the Yajiageng valley (GenBank accession number MK226959) form a monophyletic clade that exhibits less than 3% divergence from a deposited sequence of T. gonggaensis (Fu & Huang, 1991) in GenBank by Shi et al. (2016). (2) All samples collected at a lower altitude habitat at Gangbogeng valley (GenBank accession number MK226960) formed a highly supported clade with no known sister group that is less than 10% divergent. (3) The rest of the samples, including all samples collected at the habitat of T. shambalaensis and other samples collected in six different glacial valleys, form a highly supported clade (GenBank accession number MK226958) with more than 5% divergence from any known sequence of Thitarodes species.

Figure 8. 

The distribution of Thitarodes species sampled in this study and historical records of Thitarodes species in the Kangding – Mt. Gongga region.

A total of 128 out of 134 samples from seven glacial valleys was successfully sequenced according to the RAD-seq protocol (Davey and Blaxter 2010); this generated on average 4,945,132 reads (SD = 2,841,511), and 582,986 SNPs (SD = 357,473) per sample. A total of 3.9 million SNPs was counted (Fig. 9B, Suppl. material 2: Table S2). Hierarchical clustering of reads coverage showed three groups of samples sharing SNPs: (1) All six samples collected from Yajiageng valley share 0.3 of the 3.9 million SNPs. These samples’ COI sequences are similar to that of T. gonggaensis (GenBank accession number MK226959). (2) All 16 samples from the lower altitude habitat at Gangbogeng valley, collected in both 2016 and 2017, share 0.3 of the 4 million SNPs not present in any other samples. The COI sequences of these samples (GenBank accession number MK226960), form a unique group with no phylogenetically close sister taxon. (3) All 106 samples collected in valleys across the eastern side of Mt. Gongga share the majority of the 4 million SNPS that are not shared by (1) and (2). The COI sequence of this group corresponds to that of T. shambalaensis described in this study (GenBank accession number MK226958).

Figure 9. 

A Maximum Likelihood tree of CO1 sequences of species sampled in this study (yellow, red and blue dots) and other known species of the genus Thitarodes (grey dots). Outgroups are in black dots. The corresponding position of the genus Ahamus designated by Zou et al. (2010) is labeled B RAD-Seq SNP coverage of samples from three species (and degraded samples).

A total of 538 sequences (excluding three sequences from non-Thitarodes outgroups) was used to construct a summary phylogenetic tree (Suppl. material 3: Table S3, 142 CO1 sequences,163 Cytb sequences, 33 CO2 sequences, 103 wingless (wg) sequences, 6 mitogenomes, 91 multiple locus haplotypes (MLH) based on CO1). These sequences represent genetic material from at least 175 individuals belonging to the genus Thitarodes (excluding MLH sequences). Of the 175 individual samples only 42 have a species name associated, and only eight of which, including two from this study, can be verified in a publication. The best model for parsimony tree by BIC value determined by ModelFinder (Kalyaanamoorthy et al. 2017) is “TIM2+F+I+G4”. All sequences from known Thitarodes moth individuals and caterpillar fungus sclerotium extractions form one monophyletic group, clearly separated from the other Hepialidae genera Napialus Chu & Wang, 1985, Oxycanus Walker, 1856 and Phassus Walker, 1856. MLH sequences from Quan et al. (2014a) potentially overlap with individual CO1 sequences provided in Quan et al. (2014b). It is unclear whether Quan et al. (2014b) provided individual-level or population-level consensus sequences, but their list of localities is partly identical to Quan et al. (2014a). These MLH sequences (91) were dropped from the analysis. Of the resulting 184 individual samples, 154 have their closest sister group within 0.025 genetic distance from another individual. Cytb sequences of two samples, both labeled T. armoricanus, are placed into two distinct clades. The structure of the summary tree is similar to the phylogenetic tree generated from only CO1 sequences, except with the placement of samples collected from Gangbogeng valley.

Figure 10. 

Summary tree for all known species of Thitarodes. Samples in this study are labeled in color bars (blue, red and yellow). Each tip represents a cluster of all samples that are within 2.5% genetic distance from each other, the length of the grey bar represents the number of individual samples falling into the cluster. Numbers after the dots enumerate the number of samples that fall into a particular sequence type. All named individual samples are labeled. The genus Ahamus as designated by Zou et al. (2010) is labeled.

Our analysis shows that T. shambalaensis is not only morphologically, but also phylogenetically distinct from other known species in the genus Thitarodes. Even according to the generic description by Zou et al. (2010), it is still nested amongst other known Thitarodes species. Maczey et al. (2010) suggested that based on morphology, it is possible that many currently described species (such as T. gonggaensis (Fu & Huang, 1991), T. jialangensis (Yang, 1994), T. namensis (Chu & Wang, 2004), T. namlinensis (Chu & Wang, 2004), and H. hainanensis Chu & Wang, 2004) are subspecific variations of a single widespread species ranging from Sichuan to Nepal, and this hypothesized species would include the later described T. jiachaensis and T. shambalaensis. Our phylogenic analysis shows that these latter two taxa are distinct. Although the genitalic structure and wing venation of T. shambalaensis is most similar to T. jiachaensis, the two taxa are phylogenetically distant.

Some other Thitarodes species also possess triangular pseudotegumina with fan-shaped, indented sclerotization at the pseudotegumenal arms, such as T. jialangensis (Yang, 1994), T. pratensis (Yang, 1994), T. callinivalis (Liang, 1995), T. litangensis (Liang, 1995), T. kangdingensis (Chu & Wang, 1985), and T. markamensis (Liang, Li & Shen, 1992). Although morphologically distinct from T. shambalaensis (e.g., T. jialangensis Yang, 1994 has a dark forewing without spots; T. pratensis has an elongated 23^rd antenna segment and orange eyes), it is possible that these species, along with T. shambalaensis, are subspecific variations of a single widespread species ranging across the Hengduan Mountains. This hypothesis could not be further tested without a thorough revision of the genus.

Our phylogenies are largely consistent with the subdivision of the genus into Thitarodes and Ahamus by Zou et al. (2010): the latter genus, containing species such as T. jianchuanensis and T. yunnanensis (Yang, Li & Shen, 1992), form a monophyletic group in both our CO1 phylogeny and our summary phylogeny. One major inconsistency between our results and Zou et al. ’s (2010) treatment is that both T. jiachaensis and T. namensis in our phylogeny would fall into the clade of Ahamus: this is inconsistent because the defining morphological trait for genus Ahamus (Latin for “no hook”) is the lack of a hook-like structure at the ventral side of the valvae, which T. jiachaensis and T. namensis clearly possess. We note that since “no hook” is a plesiomorphic trait, it is problematic to use it to justify a monophyletic group within Thitarodes group (i.e., Ahamus). It is also possible, since CO1 sequences of T. jiachaensis and T. namensis used in this study are unpublished, that the assignment of the samples to these sequences are in error. We encourage other workers in the field to provide evidence regarding whether T. jiachaensis and T. namensis phylogenetically fall within the genus Ahamus (Zou et al. 2010) and whether there are additional, closely related groups to T. shambalaensis described in this study.

Phylogenetic placement and SNP coverage visualization from our study both suggest that some larvae and pupae samples collected from a lower elevation habitat in Gangbogeng valley across two consecutive years belong to an undescribed taxon (T. shambalaensis has been collected along this same valley, but in a habitat at a higher altitude). Our attempts to collect adults of this unknown taxon from 2016 to 2018 have not been successful, nor is it found in any other valley that we sampled. The limited range of this taxon compared to the parapatric and relatively widespread T. shambalaensis and its evolutionary history is intriguing. If the group is indeed a member of the genus Ahamus, as defined by Zou et al. (2010), we predict that the adults should not possess a hook-like structure on the ventral sides of its valvae.

CO1 sequences of our samples collected at Yajiageng valley are closely related (within 2% divergence) to known sequences labeled as T. gonggaensis in the study of Shi et al. (2016). The type locality of the species is “Kangding” (Fu and Huang 1991); Shi et al. (2016) collected the sample at Gonggashanxiang (29°32'24.0"N, 101°35'24.0"E), 65 km south west of Kangding. The locality of Shi et al. (2016) is also 57 km southwest of the samples of T. gonggaensis collected in this study, thus explaining the divergence with samples collected in our studies.

According to Wang and Yuo (2011), four other species have “Kangding” as the type locality or have been recorded in “Kangding”: T. kangdingroides (Chu & Wang, 1985), T. kangdingensis (Chu & Wang, 1985), T. oblifurcus (Chu & Wang, 1985), and T. sichuanus (Chu & Wang, 1985). It is unlikely that these samples were collected at Kangding city (30°01'21.1"N, 101°57'27.6"E). Although the city is one of the most prominent trading centers for caterpillar fungus, it is located at a relatively low elevation (2,500 m) and is not a known habitat of caterpillar fungus. These species, like samples of T. gonggaensis in this study and in Shi et al. (2016), were most likely collected in the mountains around Kangding city. Specimens of these four other species are difficult to access and no genomic data has been published at this time.

The type species of the genus, T. armoricanus, is also recorded as collected (and reared) by Tao et al. (2015), at Mt. Zheduo, 15 km northwest of Kangding, although no specimens have been dissected to verify the identification. This suggests the existence of at least eight species of Thitarodes (counting T. shambalaensis and the unknown species from Gangbogeng valley sequenced in this study) in the vicinity of Kangding-Mt. Gongga region, six of which are endemic. To our knowledge, much of the caterpillar fungus habitat around the Kangding-Mt. Gongga area, especially toward the western and north-western slope of Mt. Gongga, has not been thoroughly surveyed. We would expect more species to be discovered in this region following a systematic survey. It is hypothesized that vicariance generated by the mountain topology of the Hengduan mountains, of which Kangding is part, contributed to the rate of speciation of this biodiversity hotspot (Xing and Ree 2017). Regarding the genus Thitarodes, species diversity and endemism is not unique to the Kangding-Mt. Gongga region: the snow mountains of Yulong, Meili, Laojun in Yunnan, also part of the Hengduan mountains, is habitat to highly endemic species described to be living in close proximity (Yang 1994; Liang 1995).

A common concern in any study of the population genetics of non-model organisms is whether the analyzed samples come from different populations of the same species, or whether multiple species are involved. In our study, we still need to verify whether all Thitarodes samples collected across the eastern slope of Mt. Gongga are the same species. Without molecular evidence, this is difficult to do, since adult Thitarodes samples are difficult to identify, and most samples are obtained in larval or pupal forms, which are not sufficiently taxonomically informative for species delineation.

CO1 sequences show that our samples cluster into three clades, with inter-clade divergence falling below the 3.78% commonly observed in populations within a species, and intra-group divergence falling within the 11.06 % range of species within the same genus boundary (see review of CO1 sequence divergence by Kartavtsev 2011). We consider this strong evidence supporting for the existence of three distinct species, however, this level of interspecific (intra-population) divergence is not sufficient to distinguish patterns of population evolution history (e.g., evolutionary relationships between population of T. shambalaensis in Yanzigou valley and T. shambalaensis in Gangbogeng valley in association with the history of glaciation).

RAD‐Seq has traditionally been considered useful in population level studies, while species-level divergence would reduce the amount of shared SNP coverage across species (Davey and Blaxter 2010; Cariou et al. 2013; Eaton et al. 2017). Our results suggest that RAD‐Seq provides enough SNPs within species to further analyze the population history of each species but is not informative in analysis at the interspecific level, as samples with SNP data form three distinct clusters, but these are not shared among all samples. The clustering of those taxa with SNP data supports the CO1 species delineation indicating three main lineages.

We conclude that in this study, both the phylogeny based on CO1 sequences and the visualization of genome-wide SNP coverage provide evidence for the presence of three species in the eastern slopes of Mt. Gongga: T. shambalaensis occupies most valleys in the eastern slopes of Mt. Gongga; T. gonggaensis is found in Yajiageng valley; another undescribed species is isolated at a low elevation habitat at Gangbogeng valley.

Our summary tree shows that the genus Thitarodes is monophyletic. All known moth sequences extracted from caterpillar fungus sclerotia, despite the difficulty of assigning them to discrete species, nevertheless cluster within the genus Thitarodes.

Among the 538 sequences from 184 individuals of the genus Thitarodes, only 42 could be associated with a species name, and only eight of these are verified in publications. Of the 54 species summarized by Chu et al. (2004), only 39 have matching genetic information, most of which are unpublished Cytb sequences. 140 of the individuals are caterpillar fungus sclerotia samples, 66 of which are not within 2.5% phylogenetic distance of a known species. This suggests that many additional taxa remain to be identified and described.

In cases where multiple sources of sequences are all labeled as the same species, species identification has shown to be quite consistent: the Cytb sequence of T. yunnanensis matches to the mitogenome of T. yunnanensis, the Cytb sequence of T. jianchuanensis clusters together with CO1 sequence of T. jianchuanensis (as both are aligned with other mitogenomes), the CO1 sequence of T. gonggaensis matches its mitogenome, and the CO1 sequence of T. renzhiensis matches its mitogenome. Slight variations exist between the segment sequences and their corresponding flanking regions on the mitogenome, but we could find no phylogenetic inconsistencies. The only exception is that the two entries of the Cytb sequence of the type species of the genus, T. armoricanus, are distinctly different. While we have reason to believe that the type specimen of this sample was collected around Kangding (Chu 1965), Wang and Yao (2011) noted that this species ranges from Gansu, Tibet to Xinjiang. As we discussed in the previous section on the endemism of Thitarodes in the Hengduan mountains, we are of the opinion that the range estimate T. armoricanus of Wang and Yao (2011) includes undescribed cryptic species that have distinctive Cytb sequences, as revealed in our phylogenetic analysis.

Our summary tree also reveals inadequate sampling in many clades. When clustered by a 2.5% genetic distance (what we consider to be a reasonable threshold for inter-species variation, see Kartavtsev 2011), only 30 genetic clusters were recovered. Out of the 184 individual samples, 154 fall into 15 out of the 30 clusters; the other 15 clusters are known by either a single “orphan taxon” (T. shambalaensis, T. anomopterus, T. xiaojinensis (Tu, Ma & Zhang, 2009), T. sp. from Gangbogen in this study) or just a few sequences from a caterpillar fungus sclerotium sample. Additional research will help to elucidate the nature of these lone taxa. It is likely that “close sisters” of these taxa might be found if all known species of Thitarodes could be sequenced.

The bulk of the sampling effort so far has focused on only half of the known diversity of Thitarodes. This uneven phylogenetic sampling is not simply the result of relying on sequences of caterpillar fungus sclerotia, which have no morphological species description. The uneven sampling problem extends to the phylogeny of named species as well: the reference sequences of 39 out of the 42 named species fall into one of those 30 2.5%-distance clusters. One such cluster (including T. renzhiensis, T. oblifurcus, etc.) includes 17 known species! Most of these 17 species are known by one Cytb sequence only, which might be the result of sampling a particularly invariant region of the mitogenome. However, this issue has not been brought to attention in previous attempts to summarize genetic data for the genus (Zou et al. 2010, 2017). We noticed that in the comprehensive sampling across known range of caterpillar fungus, Zhang et al. (2014) only recovered seven well-supported clades and two “orphan taxa”. It is unclear by what genetic distance Zhang et al. (2014) have considered to be a clade, but their estimate might be a closer estimate of the number of phylogenetically well-established species in Thitarodes.

To summarize, our current understanding of the genus Thitarodes derives from three major sources:

1 The analysis of samples of caterpillar fungus sclerotia, and sometimes just Thitarodes larvae or pupae, has provided crucial information about the habitat and ecology of the species in the genus, but often fails to provide information about the morphology of the Thitarodes adult. We encourage workers in this field to include more detailed descriptions of the habitats where samples of caterpillar fungus are collected (including vegetation, climate and soil type), and make an effort to understand adult Thitarodes biology in these known habitats.

2 Many described species have types deposited in institutions around China; earlier descriptions of these samples often do not include photos of wing venation and genitalia. Revisions of these taxa with updated photos and molecular data for comparative analysis are critical.

3 New descriptions of species in the group (see Jiang et al. 2016), and NGS-based analysis of several species of the group (see Guo et al. 2016; Li et al. 2016), should provide reference sequences along with genitalia dissections, so that any known species can be linked to its position in the molecular phylogeny.