Phylogeographic structures of the host insects of Ophiocordyceps sinensis

Zoology 134 (2019) 27–37

Contents lists available at ScienceDirect

Zoology journal homepage: www.elsevier.com/locate/zool

Phylogeographic structures of the host insects of Ophiocordyceps sinensis a,b

a

a,b,c

a,b

a

Yongdong Dai , Changkui Wu , Yuanbing Wang , Yao Wang , Luodong Huang , Xijun Danga, Xuanxue Mod, Pusheng Zenge, Zhuliang Yangf, Darong Yangg, Canming Zhangg, ⁎ Paul Lemettia, Hong Yua,b,

T

a

The International Joint Research Center for Sustainable Utilization of Cordyceps Bioresources in Southeast Asia, Yunnan University, Kunming, Yunnan, 650091, China Yunnan Herbal Laboratory, Institute of Herb Biotic Resources, School of Life Sciences, Yunnan University, Kunming, Yunnan, 650091, China c Yunnan Herbal Biotech Co., Ltd., Kunming, Yunnan, 650203, China d School of Earth Sciences and Resources, China University of Geosciences, Beijing, 100083, China e National Research Center for Geoanalysis, Beijing, 100037, China f Key Laboratory for Plant Diversity and Biogeography of East Asia, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, Yunnan, 650201, China g YiKangBao Biotech Co., Ltd., Shangri-La, Yunnan, 674400, China b

ARTICLE INFO

ABSTRACT

Keywords: Ophiocordyceps sinensis Host large-scale sampling phylogeographic structure molecular clock Qinghai-Tibetan movements

A fungus-insect complex, known as DongChong XiaCao, is formed from the infection of the hepialid larvae by the fungus Ophiocordyceps sinensis, which is endemic to the Qinghai-Tibetan Plateau (QTP). Due to previously limited sample collection size, the data about the diversity and structure of the host insect was insufficient and lacked details. The purpose of this study was aimed to discuss the diversity and phylogeography of the host insects of O. sinensis with a large-scale sampling. The mitochondrial cytochrome oxidase I gene (cox1) was sequenced and analyzed among 710 samples representing 88 geographic locations. 205 haplotypes of cox1 were identified from all the 710 samples and 4 phylogenetic clades with 12 subclades were identified. Instead of the single latitude-based divergence suggested previously, three distribution patterns were deduced to correspond to the phylogeographic structures, including but not limited to the co-existence of a wide and specific local phylogeographic distribution structures. Two separate genetic diversity and differentiation centers, namely the northwestern Yunnan and the southeastern Tibet were identified. Dating analyses from three calibrations supported that the divergence of the 4 clades occurred in the Oligocene-Miocene period (30.54–13.66 million years ago) (Ma), which were connected with the second and third geological movements of the QTP (17–25, 8–13 Ma). Our results provide a more detailed understanding of the divergence and distribution patterns of the host insects of O. sinensis.

1. Introduction The entomogenous fungus Ophiocordyceps sinensis (Berk.) Sung et al. (Ascomycota, Ophiocordycipitaceae), commonly known as the Chinese caterpillar fungus, parasitizes the larvae of the ghost moth (Lepidoptera, Hepialidae) and converts it into sclerotium from which the fungus body grows. The host larva is regarded as the indispensable medium for this fungus-insect complex, which in China is known as the DongChong XiaCao (DCXC). It has been used by the traditional Chinese medicine more than 1,300 years (The Northwest Plateau Institute of Biology and Chinese Academy of Sciences and Institute for Drug Control of Qinghai, 2008; Shrestha et al., 2010; Zhang et al., 2012; Hu et al., 2013). DCXC is endemic to the Qinghai-Tibetan Plateau (QTP) and the alpine regions of the Himalayan and Hengduan Mountains in cold and

⁎

arid environmental conditions at altitudes ranging from 3000 to 5200 meters (Li et al., 2011; Dai et al., 2018). Due to the unique parasitism activity of the moth’s larva in the arid alpine environmental conditions and human over exploitation, the Ophiocordyceps sinensis complex has been designated as “endangered” in China. Recently it has been changed to “vulnerable” species in the China Biodiversity Red List 2018. Due to the limited flight range, short life cycle and geographic isolation, the host insects usually occupy specific mountain ranges and at times very narrow geographic regions. The locality of the insects in specific altitudes on the same mountain ranges has led to an abundant genetic diversity, with at least 74 species of ghost moths that have been suggested as a potential host for the fungus (Yang et al., 1996; Yang, 1997; Wang and Yao, 2011; Qiu et al., 2015; Li et al., 2017). The

Corresponding author at: No. 2 Green Lake North Road, Kunming, Yunnan, 650091, China. E-mail addresses: [email protected], [email protected] (H. Yu).

https://doi.org/10.1016/j.zool.2019.03.003 Received 24 December 2018; Received in revised form 20 March 2019; Accepted 25 March 2019 Available online 02 April 2019 0944-2006/ © 2019 Elsevier GmbH. All rights reserved.

Zoology 134 (2019) 27–37

Y. Dai, et al.

hepialid insects are in larvae instead of adult stage when infected by the fungus O. sinensis, and only the larvae’s residual intact exoskeletons remain after the mycelia sclerotium forming in the caterpillar. As larval morphology couldn’t provide sufficient identification features, the host insect is almost impossible to directly identify or link to the well-defined species. The mitochondrial cytochrome oxidase I gene (cox1) and mitochondrial cytochrome b gene (cytb) have been selected as the DNA marker for insects (Brown and Paterson, 1999; Hebert et al., 2003; Ratnasingham and Hebert, 2007; Pratheepa et al., 2014). They are considered as powerful molecular markers for identification, taxonomy and phylogeny of insect lineage. Thus, it is possible to explore the species identification and diversity of host insects by extracting the DNA from the sclerotia of O. sinensis complex. The cox1 and cytb have been widely used to discuss the species and genetic diversity of the host insects of O. sinensis, and to a certain extent, some diversity characteristics and geographical distribution pattern were identified. However, if the complex ecological geographical environmental conditions are taken into full consideration of the host insects in their wide distribution, the previous small and geographically limited sampling efforts inevitably led to insufficient results, or perhaps even some bias (Cheng et al., 2007; Zhang, 2013; Zhang et al., 2014; Quan et al., 2014). In the current study, a more extensive sampling of O. sinensis complex was carried out across the QTP. The cox1 and cytb sequences were generated and analyzed in 88 geographic locations with a total of 710 samples of the host insects derived from the sclerotia of wild O. sinensis complex. The 3 areas of focus were aimed at: (a) to describe the diversity and the phylogenetic structure of host insect of O. sinensis; (b) to explore the geographical spatial distribution patterns; (c) to estimate the divergence time followed by the discussion on the correlations between the species differentiation and the paleogeologic movements of the QTP.

2.2. Genomic DNA extraction, PCR amplification, sequencing and alignment The insects’ genomic DNA were extracted with the CTAB method (Rogers and Bendich, 1998) from the sclerotium sections, mostly from the head and leg areas of the host cadaver of DCXC. The genomic DNA (> 20 ng/ul) were used as templates to amplify DNA molecular fragments. The 2 loci of cox1 and cytb were amplified using the primer pairs shown in Table 1 (Simon et al., 1994). The polymerase chain reaction (PCR) mixture was performed in PCR Thermocycle equipment (BioRAD 100) based on the manufacturer’s instructions, and the PCR procedure was performed according to Quan et al. (2014) and Zhang et al. (2014). The PCR products were separated by electrophoresis on 1.2% agarose gels (TransGen Biotech, Beijing, China). The bands were cut and removed from the gels for purification using a Gel Band Purification Kit (BioTeke Co., Ltd, Beijing, China), then directly sequenced with an automatic sequence analyzer (ABI 7500) in Biosun Co., Ltd (Shanghai, China). The sequences were assembled and checked with SeqMan v7.1. Each locus of cox1 and cytb was aligned using MAFFT v3.7. 2.3. Genetic diversity, phylogenetic structure and phylogeographic pattern All the aligned 710 cox1 sequences were assigned to unique haplotypes in DnaSP v5.1 (Librado and Rozas, 2009). All indels were treated as single mutational events. Molecular genetic diversity was estimated by calculating haplotype diversity (Hd) and Nei’s nucleotide diversity (π) and other indices by DnaSP v5.1. All haplotypes generated from this study have been deposited in GenBank, and the new haplotypes accession numbers were MH644540–MH644592. The phylogenetic trees of the host insects were constructed by using cox1 and cytb haplotypes. The best model was calculated by using Modeltest v3.7 (Posada and Crandall, 2001). Two methods of Maximum likelihood and Bayesian probabilities were used to construct the phylogenetic structures implemented in MrBayes v3.2 (Ronquist et al., 2012) and MEGA v7.0. Likelihood settings from the best-fit model (TIM + I + G) was selected by AIC with Modeltest v3.7. The prior parameter data set was calculated by the Modeltest v3.7. Three Endoclita species were chosen as the out groups due to its relationships to the host insects in Hepialidae (Yang et al., 2015). The Analysis of Molecular Variance (AMOVA) was analyzed within and between the host insects by Arlequin v3.1 (Excoffier et al., 2007). Variance components within and among populations were calculated with 1000 random permutations. In order to explore the geographical spatial distributive patterns, the geographical mapping of haplotypes was performed in Google Earth v7.1 according to the phylogenetic structures.

2. Materials and methods 2.1. Sampling The 88 geographic locations were sampled across the wide range of the known distribution areas of O. sinensis complex, including Qinghai, Tibet, Sichuan, Gansu, and Yunnan provinces in China and three adjacent countries of Bhutan, India and Nepal. Samples were collected during the harvesting seasons (April-June), ranging from 26.46 °N–38.03 °N to 80.17 °E–103.78 °E, and at altitudes ranging from 2300–5200 meters (Fig. 1). The field survey showed O. sinensis complex had diversified habitat, mainly along the alpine scrub (dominant genus Rhododendron, Salix, Caragana and Hippophae), alpine meadow (dominant genus Potentilla, Kobresia and Pedicularis) and alpine forests (dominant genus Picea and Abies) (Li, 2012). The alpine meadows were dominant in the highland of the QTP. The alpine forests were principal characters at the northwest of Yunnan. In comparison, the alpine scrub was rich in abundance in the north of Gansu-Qinghai regions. After being dug out from the soil or obtaining them from local people, O. sinensis complex samples were placed into sample bags, placed into a portable refrigerator and later stored at 4 °C in the laboratory before analysis. All materials were stored in the Institute of Yunnan herbal biotech resources, Yunnan University in Kunming, China. This study summarized an overall data size of 710 individual samples which were collected from a total of 88 geographic locations (Supplementary data, Table S1). We collected additional 194 individual samples for DNA sequencing from 43 geographic locations, of which 26 locations were identical to previously published data (Quan et al., 2014; Zhang et al., 2014). Our collections added 17 new geographic locations to previous data, which were mainly composed in the northwestern Yunnan, the northern Qinghai, and Bhutan, Nepal. The geographic distributions were identified in Fig. 1.

2.4. The divergence time estimation The divergence time was estimated using a Bayesian method in BEAST v1.8.4 by using the mitochondrial total 13 protein encoding genes (nd2, cox1, cox2, atp8, atp6, cox3, dn3, dn5, dn4, dn4L, dn6, cytb, dn1) (Bouckaert et al., 2014). The 174 different taxa from the super-order Amphiesmenoptera were selected, whose whole mitochondrial genome had been sequenced and annotated in GenBank. Seven host insects of O. sinensis were carefully analyzed. These were Thitarodes renzhiensis and Ahamus yunnanensis (Cao et al., 2012), T. gonggaensis (Shi et al., 2015), T. xiaojinensis (Chen et al., 2015), T. pui (Yi et al., 2016), T. sejilaensis (Zou et al., 2017) and Thitarodes sp. XK2016 (Kang et al., 2017) (Table 2, Supplementary data, Table S2). In order to get a more accurate divergence time of the abovementioned lineages, the molecular clock was calculated within three calibrations. (a) Archaeolepis mane (190 Ma, standard deviation: 10%) was used as the crown calibration point for Lepidoptera (Whalley, 28

Zoology 134 (2019) 27–37

Y. Dai, et al.

Fig. 1. The geographic locations of the 88 populations of the host insects of O. sinensis. The red color represents our collections, the yellow color represents the samples from the previous research by others, the semi-red and semi-yellow represents the data from us and the previous research by others, simultaneously. All samples were collected from the Qinghai-Tibetan Plateau (QTP), Hengduan Mountains (HMs) and Himalayan Mountains at an altitude range from 3000 to 5200 m elevation.

calibration point was used to calibrate the crown group of Nymphalidae (90 Ma, standard deviation: 10%) from Wahlberg et al. (2009). The strict molecular clock was rejected (p < 0.001) and the uncorrelated relaxed clock was used. Two Markov Chain Monte Carlo (MCMC) analyses were calculated for 100 million generations with parameters being sampled every 2000 generations. Convergence of the runs were assessed using Tracer 1.5 (Rambaut and Drummond, 2009), which indicated that most parameter values had effective sample size well above 200. The two separate runs datasets were combined using LogCombiner v2.1.2 (Bouckaert et al., 2014). TreeAnnotator v2.1.2 was used to calculate node ages and the upper and lower bounds of the 95% highest posterior density interval (HPD) with the burnin of initial age of 20% trees (Bouckaert et al., 2014). The chronogram was visualized using Figtree v1.5 (Rambaut, 2012).

Table 1 The primers used for DNA amplification and sequencing of the host insects of O. sinensis. Product Primer name

Primer sequence (5’–3’)

cox1

GGTCAACAAATCATAAAGATATTGG Simon et al. (1994) TAAACTTCAGGGTGACCAAAAAATCA TATGTAYTACCWTGRGGRCAAATATC ATWACACCYCCTAATTTATTAGGRAT

cytb

Hep-cox1F Hep- cox1R Hep-cytbF Hep-cytbR

Reference

1986); (b) the Early Cretaceous unnamed fossil was assigned to Gelechiidae reported by Ross et al. (2010) that was used to calibrate the age of the crown clade within the superfamily Gelechioidea (120 Ma, standard deviation: 10%) (Wahlberg et al., 2013); (c) a secondary Table 2 The 178 taxa information used to estimate divergence time of Lepidoptera. Order Coleoptera Amphies-menoptera

Trichoptera Lepidoptera

Sub- Order

Intra-order

Zeuglop-tera Glossata, Neolepido-ptera

29

Exoporia Heteron-eura

Super-family

Number

Hepialoidea Nepticuloidea Bombycoidea Gelechioidea Geometroidea Hesperioidea Noctuoidea Papilionoidea Pyraloidea Tineoidea Tortricoidea Yponomeutoidea Zygaenoidea Others 7 super-family

4 13 1 9 2 20 9 8 10 13 40 17 5 13 3 4 7

Zoology 134 (2019) 27–37

Y. Dai, et al.

Fig. 2. The geographical distribution patterns of the gradient diversity difference of host insects of O. sinensis. Five gradients have been set (Hd = 0, 0–0.3, 0.3–0.6, 0.6–0.9, 0.9–1.0). Different colors and shapes were used to mark the gradients.

3. Results

phylogenetic status of these 27 species were further clarified (Fig. 3b). Thitarodes armoricanus, T. kangdingroides, T. gonggaensis, T. renzhiensis, T. baimaensis, T. ferrugineus, T. yushuensis, T. xiaojinensis, T. oblifurcus were clustered into Clade I-1. T. cingulatus were grouped to Clade I-2. While Thitarodes sejilaensis and T. pui were contained in Clade II-2 and II-4, separately, while Ahamus menyuanicus, A. luquensis, A. yunnanensis and A. jianchuanensis were clustered into Clade IV. In general, based on the phylogeny of the large-scale samplings of the host inscts of O. sinensis (710 samples), the four clades of host insects were all assigned to Thitarodes and Ahamus (Lepidoptera, Hepialidae) (Zou et al., 2010). Among them, Clade I, II and III belonged to Thitarodes, while Clade IV was identified within the genus Ahamus. Overall, the identification of the four clades uncovered that the host insects of O. sinensis were not monophyletic but belonged to the two genera of Thitarodes and Ahamus in Hepialidae. while the hosts were separated by the non-parasite species in Ahamus.

3.1. The genetic diversity and phylogenetic structure A total of 710 cox1 sequence data was obtained from 710 individual samples from 88 geographic locations. 205 haplotypes were identified (Supplementary data, Table S3). All haplotypes were deposited in the GenBank, and the 53 new haplotypes are available within the accession numbers MH644540–MH644592. Compared with the parasite fungus O. sinensis (Dai et al., 2018), the host insects exhibited higher level of genetic diversity [Hd = 0.9307 vs 0.8035 (fungus); π = 0.04856 vs 0.01137 (fungus)]. 57 out of the 88 locations (64.77%) had haplotype and nucleotide diversity. The high level of haplotype diversity (Hap > 10) coupled with high nucleotide diversity (π > 0.05) were identified in three locations (Gongbo'gyamda, Nyingchi in Tibet) and Deqin in Yunnan). The comparison of all samples data identified only one haplotype within each of the 10 populations, and no genetic diversity was observed within the populations, which includes Minle (Gansu), Menyuan, Huangzhong, Qilian, Maqen, Qumalai, Henan (Qinghai), Riwoqe, Lasaduilong Deqin, Shigatse Dongga (Tibet) (Hd = 0, π = 0) (Fig. 2; Supplementary data, Table S1). The four divergent clades (Clade I-IV) with 12 subclades (Clade I-1 to Clade IV-3) were presented based on the phylogenetic structure with the 205 haplotypes (Table 3; Fig. 3a). The results uncovered the basic diversity and differentiation characteristics of the host insects. Additionally, a total of 168 cytb sequences, which represented 67 geographical locations that selectively generated from the 12 subclades were assigned by cox1 phylogeny (Supplementary data, Table S4). A total of 74 haplotypes were identified that exhibited a high genetic diversity (Hd = 0.947, π = 0.05329). We combined the cox1 and cytb sequences of 168 sample data. 118 new haplotypes were identified. The phylogenetic structure of these haplotype dataset was high in accordance with that of the cox1, presenting 12 branches (Fig. 3b). The subtle difference was that the position of Clade III, that was near to Clade IV, had a close relationship to Clade II within the single cox1 tree. The 27 well-defined Thitarodes and Ahamus species were added while constructing the 2-locus phylogenetic tree. Therefore, the

3.2. The geographical and phylogeographic distribution patterns In order to explore the ecological geographical variation patterns of the host insect of O. sinensis, the genetic differentiation was analyzed within each 9 geographic regions (Table 4; Fig. 4). Except Sect.⑦ (only 1 sample was obtained) each of the 8 regions exhibited genetic diversity. Thus, it was concluded that the southeast of Tibet (Sect.⑤) contained the most diversity or seven subclades where four were endemic (Clade II-1, II-2, II-4, and II-5). This was followed by the northwest of Yunnan (Sect.⑥), which contained four subclades, where two were endemic (Clade IV-1 and IV-3). In contrast, low diversity regions were represented in the geographic Sect.②, ③ and ④. The results revealed inconsistent genetic differentiation among the nine geographic regions but suggested that the modern differentiation center of host insects of O. sinensis occurred in the southeast of Tibet. However, the gene flow was not blocked and showed some genetic differentiation in the highland of the QTP (Sect. ③ and ④). To explore the genetic structure and the spatial genetic pattern of the host insects of O. sinensis, we carried out the AMOVA analyses with phylogenetic clades and geographic regions to definite the main source of the genetic variation. The results suggested that there wasn't any significant differentiation between the 9 geographic regions 30

Zoology 134 (2019) 27–37

Y. Dai, et al.

Table 3 Genetic diversity within 12 subclades of the host insects of O. sinensis by cox1 sequences. Clade

No.of pop.

No. of seq.

Polymorph.

Hap

Hd

π

Fst

Nm

I-1 I-2 II-1 II-2 II-3 II-4 II-5 III-1 III-2 IV-1 IV-2 IV-3 Whole

68 4 3 4 4 3 1 1 4 6 17 3

450 11 10 17 10 24 14 9 14 25 122 4 710

99 52 49 26 15 27 14 21 44 38 43 13 199

117 10 6 8 5 8 7 4 7 13 17 3 205

0.881 0.982 0.911 0.860 0.667 0.801 0.846 0.778 0.846 0.917 0.287 0.833 0.9307

0.01209 0.02752 0.03604 0.01158 0.00875 0.01045 0.00676 0.01903 0.03041 0.01574 0.00401 0.01484 0.04856

0.6360 0.0391 0.4851 0.5371 0.7208 0.2534 – – 0.8287 0.5022 0.9144 0.9999 0.6664

0.1431 6.1439 0.2654 0.2155 0.0968 0.7366 – – 0.0517 0.2478 0.0234 0.0001 0.1252

88

(Fct = 0.3976), whereas the differentiation of the 12 divergent subclades was the main source of genetic variation (Fct = 0.8381) (Table 5). The mapped the geographical distributions of each 12 subclades were shown in (Fig. 5). Subclade I-1 was a ubiquitous group, which was represented across the vast distribution areas of DCXC, including the location of Pithoragarh-India (the west edge), Wenxian-Gansu (the east edge), Huangzhong-Qinghai (the north) and Xuebang snow MountainsYunnan (the south edge). Subclade I-2 distributed in four locations (Luqu, Zhouqu in Gansu, Ngawa, Xiaojin in Sichuan), the border of Gansu and Sichuan. The geographic locations were on the eastern edge of the DCXC distribution areas. All samples of Subclade II-1 to II-5 came from Nyingchi and Lhoka prefecture of the southeastern Tibet. There were great overlaps in geographical distribution among these 5 subclades, indicating dramatic differentiation in these regions, in particular the Subclade II-5 was only found in Nyingchi. Subclade III-1 was found only in the northwestern part of Nepal, and Subclade III-2 was mainly distributed from Bhutan and Naidong (Tibet) while both of the two subclades were narrowly distributed in the southern margin. Subclade IV-1 was mainly distributed in the northwestern Yunnan, exhibiting geographic specificity. The distribution of Subclade IV-2 presented a complex geographical structure. The samples were mainly identified in the surrounding area of the Qinghai Lake, as well as found in Tarlag, Chindu (The middle of Qinghai), Ngawa, Litang (Sichuan), Baima snow Mountains, Haba snow Mountains and Lanping (Yunnan). This specific distribution presented a long-distance dispersal distributive structure. Subclade IV-3 were from three locations of Zayu, Biluo snow Mountains, Xuebang snow Mountains as well as along the south edge of the distributions. In general, based on the phylogeographic structures of the 12 subclades, 3 distribution patterns can be summarized: (a) the wide distribution around the QTP (Subclade I-1), suggests that not all of the Thitarodes species were narrowly distributed as the mountain-river geographical barriers didn’t always lead to genetic isolation or speciation of the host insects; (b) the long-distance dispersal with latitudebased structure (Subclade IV-2); and (c) the specific local geographic locations with latitude-based divergence (Subclade I-2 and IV-1) and longitude-based divergence (Subclade II-1 - II-5, III-1, III-2 and IV-1).

xiaojingensis, T. gonggaensis, T. renzhiensis and Thitarodes sp. XK-2016 were clustered into Clade I. Ahamus yunnanensis may not be the host, but it belonged to Clade IV and had a close relationship to the confirmed hosts in the genus Ahamus. Based on the three calibrations, the divergence time of each node was estimated from Amphiesmenoptera- Lepidoptera (order) – Exoporia and Heteroneura, and the crown of the order Lepidoptera was presented at 191.70 Ma (95% HPD: 177.56–206.56) (Table 6). The divergence times for the seven species of Thitarodes and Ahamus were shown in (Fig. 7). Ahamus yunnanensis was firstly differentiated from these host cluster groups in the middle Oligocene (30.54 Ma, 95% HPD = 24.16–38.40), which represented a possible divergence time for the Clade IV and the other three clades. The calculated differentiation time was approximately 25.26 Ma (95% HPD = 19.36–32.04) between the monophyletic group of T. sejilaensis, T. pui and T. xiaojingensis, T. gonggaensis, T. renzhiensis, Thitarodes sp. XK-2016. It provided a rough divergence time between the Clade I and II. As having occurred around 13.97 Ma (95% HPD = 10.36–18.17) between T. xiaojingensis and T. gonggaensis, T. renzhiensis, Thitarodes sp. XK-2016, it allowed a rough divergence time definition of the crown Clade I-1. In comparison, the divergence time of T. sejilaensis and T. puii occurred around 13.66 Ma (95% HPD = 9.57–18.34), which also could represent the divergence of the crown clade II. In general, the divergence of the four clades was calculated between 30.54–13.66 Ma, which was coinsides with the second and third uplift movements in the QTP (25–17, 13–8 Ma) (Harrison et al., 1992; Chen, 1992; Zhong and Ding, 1996; Li and Fang, 1998; Li et al., 2001). The differentiations of host insects were consistent with the geological events of the QTP during the Oligocene-Miocene periods. 4. Discussion 4.1. A large-scale sampling explained the diversity and the differentiation of the host insects of O. sinensis Due to the disproportionate sampling in each of the nine geographic regions, there was some bias regarding the previous ideas for diversity and differentiation of host insect derived from the wild O. sinensis complex (Zhang et al., 2014; Quan et al., 2014). In this study, a largescale sampling of O. sinensis complex was conducted around the wide range of the distribution areas. The previously under-represented regions, such as the northwest of Yunnan, was one of our targeted focus. A total of 710 cox1 sequences were obtained representing 88 geographic locations of the host insects of O. sinensis, which ensured an abundant sample size and the additional data in order to fill the gaps. A relatively balanced sample size was kept for each of the nine geographic regions to prevent the misaligned genetic statistics providing more reliability for the overall analyses (Table 4). A relatively high diversity was uncovered among the 710 samples of

3.3. Divergence time estimation The divergence time was estimated for the order of Lepidoptera with emphasis on the differentiation of the host insects of O. sinensis based on the mitochondrial 11 protein encoding genes (atp8 and dn6 were excluded in the final analyses because of the big variation) (Fig. 6). The following six species were confirmed to be hosts for O. sinensis. Thitarodes sejilaensis and T. pui were assigned to Clade II, and T. 31

Zoology 134 (2019) 27–37

Y. Dai, et al.

Fig. 3. The phylogenentic tree of host insects of O. sinensis constructed by the haplotypes of cox1 and 2-locus of cox1 and cytb. The phylogenetic tree of 205 haplotypes of cox1 from 710 samples data; (b) The phylogenetic tree of 111 haplotypes and 35 well-defined host species of 2-locus (cox1 and cytb) from 203 samples data. Both phylogenetic trees were constructed with Maximum Likelihood and Bayesian proportion methods.

Table 4 Genetic diversity indices of the host insects of O. sinensis among 9 geographical regions.

Population No.of sequences polymorphic sites Hap Hap group Hd π

All

Sect. ①

Sect.②

Sect.③

Sect. ④

Sect.⑤

Sect.⑥

Sect.⑦

Sect. ⑧

Sect. ⑨

88 710 199 205 12 0.9307 0.04856

11 134 54 14 2 0.391 0.02247

12 59 108 31 3 0.944 0.03886

17 170 80 23 3 0.521 0.00776

11 84 33 23 1 0.629 0.00336

14 152 142 68 7 0.951 0.06088

14 88 108 51 4 0.980 0.03663

1 1 – 1 1 – –

5 14 56 8 2 0.901 0.04633

3 8 64 4 2 0.643 0.03559

32

Zoology 134 (2019) 27–37

Y. Dai, et al.

Fig. 4. The genetic diversity among the 9 geographic regions of the host insects of O. sinensis. Sect.①. the north Qinghai-Gansu regions (surrounding Qinghai Lake); Sect. ②. the west Sichuan regions (the north HMs); Sect.③. the south-middle Qinghai regions (the north QTP); Sect.④. the northeast of Tibet regions (the south QTP); Sect.⑤. southeast Tibet regions (the southwest HMs); Sect.⑥. northwest Yunnan regions (the south HMs); Sect.⑦. the west Himalayas regions; Sect. ⑧. the middle Himalayas regions; Sect. ⑨. the east Himalayas regions. Different colors dots were used to mark the different geographic regions. Table 5 AMOVA analyses among the geographical regions and the phylogenetic clades of the host insects of O. sinensis.

9 geographic sections

4 phylogenetic clades

12 phylogenetic subclades

Source of variation

d.f.

Sum of squares

Variance components

Percentage of variation

F- statistics

Among 9 sections Among populations within sections Within populations Total Among 4 clades Among populations within clades Within populations Total Among 12 subclades Among populations within sunclades Within populations Total

8 79 622 709 3 114 592 709 11 106 592 709

3946.435 3127.788 2978.600 10052.823 6177.934 3039.620 835.269 10052.823 7694.723 1522.831 835.269 10052.823

6.09250 Va 4.44208 Vb 4.78875 Vc 15.32332 16.37275 Va 4.28467 Vb 1.41093 Vc 22.06834 18.73984 Va 2.21021 Vb 1.41093 Vc 22.36089

39.76 28.99 31.25

Fct = 0.3976 Fsc = 0.4812 Fst = 0.6875

74.19 19.42 6.39

Fct = 0.7419 Fsc = 0.9361 Fst = 0.7523

83.81 9.88 6.31

Fct = 0.8381 Fsc = 0.6104 Fst = 0.9369

Fct: genetic differences among groups defined a priori; Fsc: genetic differences among population within groups; Fst: genetic differences among populations.

host insects (Hd = 0.9307; π = 0.04856). Four distinct clades with 12 subclades were identified by phylogenetic reconstruction, uncovering the basic differentiation status of the host insects.

world-class biodiversity hotspot in both of these areas and having complex and unique natural geographical conditions, the insects could gather, exchange and differentiate (Sun, 2002). These two areas were suggested to be the glacial refugia for the alpine hepialid species and effectively avoided extinction during the glacial period of the Quaternary. The genetic diversity was well protected by the local conditions. After the glacial period, the hepialid species re-expanded to higher altitudes and presented new spatial distribution structures.

4.2. The diversity and differentiation centers occurred both in northwest Yunnan and southeast Tibet In the previous studies, the southeast of Tibet, especially Nyingchi distinct was identified as a genetic diversity of the host insects of O. sinensis (Zhang et al., 2014; Quan et al., 2014). However, the data did not clarify the importance of the northwestern region of Yunnan due to limited and unbalanced samplings. For example, only 17 samples of two locations were collected from the northwestern Yunnan in Quan region (2014), while we collected 88 samples representing 14 locations in total. Compared to the previous result, the higher genetic diversity regions were identified in both the northwest of Yunnan (Hd = 0.980, π = 0.03663, subclades = 4) and the southeast of Tibet (Hd = 0.951, π = 0.06088, subclades = 7). It was demonstrated that both regions were the diversity and differentiation centers for the host insects. As

4.3. At least three geographical distribution patterns explained the phylogeography of the host insects of O. sinensis In this study, we analyzed the distributions patterns of 710 samples from 12 divergent subclades, and at least three geographic distribution patterns were discovered to correspond to the following phylogeographic structures: (a) the wide distributions around the QTP; (b) the long-distance dispersal with latitude-based distributions; and (c) the specific local distributions. It has been previously verified that there was geographical variation exhibiting latitude-based population divergence of the host insects of O. 33

Zoology 134 (2019) 27–37

Y. Dai, et al.

Fig. 5. Phylogeographic structures of the 12 subclades of host insects of O. sinensis. (a) The geographical distribution overall structure of the 12 subclades; (b) The geographical distribution structure of Subclades I-1, I-2, II-1, IV-1; (c) The geographical distribution structure of Subclades II-2, II-3, IV-2, IV-3; (d) The geographical distribution structure of Subclades II-4, II-5, III-1, III-2. Different colors dots were used to mark the different subclades.

sinensis (Quan et al., 2014). The geographic distribution of Subclade I-2 and IV-1 were in accordance with this divergence pattern. However, our results also supported that there had at least three geographic distribution patterns, instead of a single latitude-based structure as was previously suggested. The wild distribution structure type was displayed from the dominant haplotypes (Hap.1 and Hap.6) and the phylogenetic subclade I-1, which strongly revealed the complex mountain-river geographical barriers in Hengduan Mountains didn’t always lead to genetic isolation and speciation of the host insects. In additional, the long-distance dispersal structure (Subclade IV-2) was also discovered. It is possible that

there might be an accessible dispersal channel around the edge of regional areas of the host insects. Finally, the specific local geographic locations revealed that the geographic and environmental conditions caused the speciation and differentiation of the host insects. 4.4. The differentiation of the host insects of O. sinensis strongly linked to geological movements of the QTP The uplift of the Himalayas was a complex process characterized by poly-stages, non-uniform velocity and nonuniformity, which has been divided into four stages: 45–38 Ma, 25–17 Ma, 13–8 Ma, and 3 Ma up to Fig. 6. Time-calibrated phylogeny of Lepidoptera lineages obtained from Bayesian analyses. The divergence time was estimated using a Bayesian method in BEAST v1.8.4 by using the mitochondrial total 13 protein encoding genes (nd2, cox1, cox2, atp8, atp6, cox3, dn3, dn5, dn4, dn4L, dn6, cytb, dn1) (Bouckaert et al., 2014). The 174 different taxa from the super-order Amphiesmenoptera were aligned and analyzed.

34

Zoology 134 (2019) 27–37

Y. Dai, et al.

Table 6 Divergence time of the main Lepidoptera lineages based on 3 calibrations. Node

Lineage

Time (Ma)

95% HPD

Fossil calibration

Wahlberg et al. (2013)

1 2* 3

Amphiesmenoptera Lepidoptera (Exoporia-Heteroneura) Exoporia crown (The crown Endoclita signifier- Napialus hunanensis- Thitarodes –Ahamus) Heteroneura, Ditrysia crown Gelechioidea crown Papilionoidea crown Pyraloidea crown Bombycoidea crown Noctuoidea crown Hesperioidea crown Geometroidea crown Tortricoidea crown Nymphalidae 7 Thitarodes –Ahamus (the crown group of host insects of O. sinensis) 6 Thitarodes except Ahamus yunnanensis 4 Thitarodes xiaojingensis - Thitarodes gonggaensis/renzhiensis/sp. XK-2016

322.54 191.70 61.00

287.59–358.73 177.58–206.55 44.12–77.04

190 (sd = 10%)

190 (180–200)

153.61 106.20 105.82 103.13 93.01 88.63 81.97 74.80 64.18 87.52 30.54 25.22 13.97

140.44–166.24 96.53–115.58 98.83–113.16 94.13–110.91 84.45–101.49 79.75–97.36 72.44–91.53 65.08–84.50 55.00–74.59 80.34–94.70 24.16–38.40 19.36–32.04 10.36–18.17

120 (sd = 10%)

106 (95–117) 104 (95–114) 93 (80–105) 84 (74–93) 82 (73–92)

4 5* 6 7 8 9 10 11 12 13* 14 15 16

90 (sd = 10%)

83 (72–93) 68 (52–86) 90 (80–100)

The highlight with "bold" represent the host insects taxa and the higher orthologous lineages of this groups. * Represent three time calibrations.

present time. The data also showed that large-scale intra-continental convergence took place approximately 30 Ma around the Indian and Asian continents (Zhong and Ding, 1996; Yuan, 2003). Due to the lack of fossil records on Hepialidae, the standard insect mtDNA clock (2.3% pairwise sequence divergence per million years) was commonly used to estimate the divergence time (Brower, 1994). However, different treatments confuse estimates of mtDNA substitution rates in opposing directions and obscure lineage-specific differences in rates (Papadopoulou et al., 2010). In this study, the divergence time of Lepidoptera lineages with the focus on the Hepialidae clade were estimated based on the 11 mitochondrial genes. The divergence time of the eight super-families in Heteroneura were all consistent with the result of Wahlberg et al. (2013). As such, we deduced a more comprehensive understanding of the evolutionary of the host insects of O. sinensis based on multiple calibrations. The divergence between Ahamus (Clade IV) and Thitarodes (Clade I-III) occurred approximately during 30.54 Ma. The differentiation came about 25.26 Ma between Clade I and II and the differentiation of crown Clade I and II occurred nearly at the same time (I: 13.97 Ma; II: 13.66 Ma).

The divergence was consistent with the geological strike-slip and the second and third uplift movements in the QTP (Harrison et al., 1992; Chen, 1992; Zhong and Ding, 1996; Li and Fang, 1998; Li et al., 2001). We believed that the QTP geological movements were an important force for the differentiation of the host insects of O. sinensis. However, the differentiation of the host insects of O. sinensis was connected with geological event in the QTP, which made this group a good indicator for the geochronology. We deduced a new way of thinking and methodology to estimate the time of geological changes, whereby it is possible to cross-examine the geochronology based on paleogeology and molecular biology at a more precise scale. 4.5. The host insects of O. sinensis were not monophyletic but belonged to two genera (Thitarodes and Ahamus) in Hepialidae based on large sampling Based on the criterion of the overlaps at the distribution ranges of the altitudes between the insects and the fungus O. sinensis, at least 74 hepialid species were proposed to the host of O. sinensis that belong to the Thitarodes and Ahamus genus and some were members of Pharmacis, Magnificus, Bipectilus, Endoclita, Gazoryctra, Parahepialus (Wang and

Fig. 7. Correlations between the differentiation of the host insects of Ophiocordyceps sinensis and the geological uplift of the QTP. Three phylogenetic divergence of the host insects were consistent with the geological strike-slip and the second and third uplift movements in the QTP. 35

Zoology 134 (2019) 27–37

Y. Dai, et al.

Yao, 2011; Li et al., 2017). However, due to the characteristics of our research methodology, it was speculated rather than confirmed that the host insects can be defined as “potential” vs. “definite”. Although 88 samples were collected in northwestern Yunnan, the distribution regions of A. yunnanensis and A. jianchuanensis, we found no experimental data to confirm that A. yunnanensis and A. jianchuanensis acted as the host of O. sinensis (Fig. 3). This experimental evidence confirmed that the overlapping criterion methodology was not practical to recognize the host insects of O. sinensis. In our study, a large-scale sampling of 710 data samples were conducted around the wide distribution of DCXC complex. It was shown that all the host insects did belong to Thitarodes and Ahamus, and no others insect genus was found to host of O. sinensis. It is suggested, therefore that the present assumed species diversity has overestimated the host insects of O. sinensis. However, there isn't enough experimental data nor practical theory to support at this time any definitive conclusions.

palaeogeologic history of pan-Himalayas. We thank Yu L. Li (Qinghai University) for helpful scheduling suggestions for Qinghai sampling in 2017. We thank Yong J. Zhang (Shanxi University), and Qing M. Quan (Tongji University) for cox1 and cytb sequences with locality information. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.zool.2019.03.003. References Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.H., Xie, D., Drummon, A.J., 2014. Beast 2: a software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 10, e1003537. Brower, A.V., 1994. Rapid morphological radiation and convergence among races of the butterfly Heliconius erato inferred from patterns of mitochondrial DNA evolution. Proc. Natl. Acad. Sci. USA 91, 6491–6495. Brown, B., Paterson, A.R., 1999. Mitochondrial COI and II provide useful markers for Wiseana Lepidoptera: Hepialidae. species identification. Bull. Entomol. Res. 89, 287–293. Cao, Y.Q., Ma, C., Chen, J.Y., Yang, D.R., 2012. The complete mitochondrial genomes of two ghost moths, Thitarodes renzhiensis and Thitarodes yunnanensis: the ancestral gene arrangement in Lepidoptera. BMC Genomics 13, 1–13. Chen, F.B., 1992. Hengduan event, an important tectonic event of the late Cenozic in eastern Asia. Mountain Res. 10, 195–202. Chen, S., Shi, P., Zhang, D., Lu, Z., Li, L., 2015. Complete mitochondrial genome of Hepialus xiaojinensis Lepidoptera: Hepialidae. Dna Seq. 28, 305–306. Cheng, Z., Geng, Y., Liang, H.H., Yang, X.L., Li, S., Zhu, Y.G., Guo, G.P., Zhou, T.S., Chen, J.K., 2007. The phylogenetic relationship between cordyceps sinensis and its host, bat moth, was studied using mitochondrial Cytb gene sequence. Progress in Nat. Sci. 17, 1045–1052. Dai, Y.D., Wu, C.K., Yuan, F., Wang, Y.B., Huang, L.D., Chen, Z.H., Zeng, W.B., Wang, Y., Yang, Z.L., Mo, X.X., Lemetti, P., Yu, H., 2018. Genetic structure and biogeography of an alpine endemic caterpillar fungus Ophiocordyceps sinensis. BMC Evol. Biol (Unpublished results). Excoffier, L., Laval, G., Schneider, S., 2007. Arlequin (version 3.0), An integrated software package for population genetics data analysis. Evol. Bioinform. Online 1, 47. Harrison, T.M., Copeland, P., Kidd, W.S., Yin, A., 1992. Raising Tibet. Science 255, 1663–1670. Hebert, P.D., Ratnasingham, S., Dewaard, J.R., 2003. Barcoding animal life: cytochrome c oxidase subunit 1 divergences among closely related species. P. Roy. Soc. B-Biol. Sci 270, 96–99. Hu, X., Zhang, Y.J., Xiao, G.H., Zheng, P., Xia, Y.L., Zhang, X.Y., St. Leger, R.J., Liu, X.Z., Wang, C.S., 2013. Genome survey uncovers the secrets of sex and lifestyle in caterpillar fungus. Chi. Sci. Bul 58, 2846–2854. Kang, X., Hu, Y., Hu, J., Hu, L., Wang, F., Liu, D., 2017. The mitochondrial genome of the Lepidopteran host cadaverThitarodes sp. of Ophiocordyceps sinensis, and related phylogenetic analysis. Gene 598, 32–42. Li, H., 2012. The Resources of Ophiocordyceps sinensis in Tibet, China. Yunnan Science and Technology Press. Li, J.J., Fang, X.M., 1998. Research on the uplift of Qinghai - Tibet plateau and its environment change. Chin. Sci. Bull. 43, 1569–1574. Li, J.J., Fang, X.M., Pan, B.T., Zhao, Z.J., Song, Y.G., 2001. Late cenozoic intensive uplift of Qinghai-Xizang plateau and its impacts on environments in surrounding area. Quaternary Sci. 21, 381–391. Li, Y., Wang, X.L., Jiao, L., Jiang, Y., Li, H., Jiang, S.P., Lhssumtseiring, N., Fu, S.Z., Dong, C.H., Zhan, Y., Yao, Y.J., 2011. A survey of the geographic distribution of Ophiocordyceps sinensis. J. Microbiol. 49, 913–919. Li, W.J., Zhang, Z.Y., Li, Q.P., Lyu, Y.H., Jiang, C.J., Wei, Z.H., 2017. Research progress of host insects of Cordyceps and its raising technique. World Chin. Med. 12, 3142–3150. Librado, P., Rozas, J., 2009. DnaSP v5, a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. Papadopoulou, A., Anastasiou, I., Vogler, A.P., 2010. Revisiting the insect mitochondrial molecular clock: the mid-Aegean trench calibration. Mol. Biol. Evol. 27, 1659–1672. Posada, D., Crandall, K.A., 2001. Selecting the best-fit model of nucleotide substitution. Syst. Biol. 50, 580–601. Pratheepa, M., Jalali, S.K., Arokiaraj, R.S., Venkatesan, T., Nagesh, M., Panda, M., Pattar, S., 2014. Insect barcode information system. Bioinformation 10, 98–100. Qiu, Y., Chen, Y.L., Peng, C., Wan, D.G., Shen, C.H., Yi, B., Hou, F.X., Guo, J.L., 2015. The advance of the host insects of Ophiocordyceps sinensis in China. Lishizhen Medicine and Materia Medica Research 26, 720–722. Quan, Q.M., Chen, L.L., Wang, X., Li, S., Yang, X.L., Zhu, Y.G., Wang, M., Cheng, Z., 2014. Genetic diversity and distribution patterns of host insects of caterpillar fungus Ophiocordyceps sinensis in the Qinghai-Tibet plateau. Plos One 9, e92293. Rambaut, A., 2012. FigTree. Tree Figure Drawing Tool, v. 1.4.0. Available from:. http:// tree.bio.ed.ac.uk/software/figtree/. Rambaut, A., Drummond, A., 2009. Tracer: MCMC Trace Analysis Tool, Version 1.5.0. Available from:. http://tree.bio.ed.ac.uk/software/tracer/UniversityofOxford. Ratnasingham, S., Hebert, P.D.N., 2007. Bold: the barcode of life data system www.barcodinglife.org. Mol. Ecol. Notes 7, 355–364.

5. Conclusions A total of 710 samples data were collected which represent 88 geographic locations of the host insects derived from O. sinensis complex. A total of 205 haplotypes were identified and clustered into 4 clades with 12 subclades, which fell into the two genera of Thitarodes and Ahamus in Hepialidae (Lepidoptera) instead of a single monophyletic clade. Three spatial geographical distribution patterns contributed to the genetic differentiation of the host insects within their coexistence of wide, long-distance dispersal and the specific geographical distribution. The host insects experienced a non-single mode of diffusion. Both the northwestern of Yunnan and southeastern regions of Tibet are considered as the centers of diversity and differentiation. Both have been suggested to coincide with the glacial refugia in the Quaternary. Dating analyses based on the three calibrations that support the divergence of the four clades of the host insects that occurred in the Oligocene-Miocene period (30.54–13.66 Ma). This was strongly related to the geological strike-slip and the uplift movements of the QTP, suggesting that the geographical movements had a large influence on the differentiation of the host insects of O. sinensis. At least three geographical distribution patterns were identified within the phylogeographic structures. However, the exact mechanism of the geographical distributive patterns is still subject to debate and may require additional studies. We maintain a non-conclusive view about the specific mechanism about the phylogeographic distribution of the host insects. The exact pattern of distribution does require additional understanding of geology, stratigraphic structures, geochemical elements and the ecological factors. Declaration of interest None. Funding This work was supported by the National Natural Science Foundation of China (31870017), the Yunnan Science and Technology Foundation of the International Joint Research Center for Sustainable Utilization of Cordyceps Bioresources in Southeast Asia (KC1810172), the National Natural Science Foundation of China (K1021110), the Science and Technology Development Fund of Guidance from the Central Government to Locals (KC1610530), and the China Postdoctoral Science Foundation (2017 M613017). Acknowledgements We thank De Y. Hong (Beijing institute of Botany, Chinese Academy of Sciences) for his kindly suggestion of the diversity and 36

Zoology 134 (2019) 27–37

Y. Dai, et al. Rogers, S.O., Bendich, A.J., 1988. Extraction of DNA from Plant Tissues. Plant Molecular Biology Manual A6. Kluwer Academic Publishers, Dordrecht, pp. 1–10. Ronquist, F., Teslenko, M., Mark, P.V.D., Ayres, D.L., Darling, A., Höhna, S., 2012. Mrbayes 3.2, efficient bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 61, 539–542. Ross, A.J., Mellish, C., York, P., Crighton, B., 2010. Burmese Amber. In: Penny, D. (Ed.), Biodiversity of Fossils in Amber from the Major World Deposits. Siri Scientific Press, Manchester, UK, pp. 208–235. Shi, P., Lu, Z., He, Y., Chen, S., Qin, S., Qing, Y., Yan, J., 2015. Complete mitochondrial genome of Hepialus gonggaensis Lepidoptera: Hepialidae., the host insect of Ophiocordyceps sinensis. Dna Seq. 27, 4205–4206. Shrestha, B., Zhang, W.M., Zhang, Y.J., Liu, X.Z., 2010. What is the Chinese caterpillar fungus Ophiocordyceps sinensis (Ophiocordycipitaceae)? Mycology 1, 228–236. Simon, C., Frati, F., Beckenbach, A., Crespi, B., Liu, H., Flook, P., 1994. Evolution, weighting, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Ann. Entomol. Soc. Am. 87, 651–701. Sun, H., 2002. Tethys retreat and Himalayas-Hengduanshan mountains uplift and their significance on the origin and development of the Sino-Himalayan elements and alpine flora. Acta Botanica Yunnanica 24, 273–288. The Northwest Plateau Institute of Biology, Chinese Academy of Sciences & Institute for Drug Control of Qinghai, 2008. Zhong Guo Chong Cao: History, Resources and Scientific Research. Shaanxi Science and Technology Press. Wahlberg, N., Leneveu, J., Kodandaramaiah, U., Pena, C., Nylin, S., Freitas, A.V.L., Brower, A.V.Z., 2009. Nymphalid butterflies diversify following near demise at the Cretaceous/Tertiary boundary. P. Roy. Soc. B-Biol. Sci. 276, 4295–4302. Wahlberg, N., Wheat, C.W., Peña, C., 2013. Timing and patterns in the taxonomic diversification of Lepidoptera Butterflies and Moths. Plos One 8, e80875. Wang, X.L., Yao, Y.J., 2011. Host insect species of Ophiocordyceps sinensis: a review. ZooKeys 127, 43–59. Whalley, P., 1986. A review of the current fossil evidence of Lepidoptera in the Mesozoic.

Biol. J. Linn. Soc. Lond. 28, 253–271. Yang, D.R., 1997. The relationships between species diversity in Hepialidae and environmental evolution in the Qinghai-Tibet Plateau. Acta Biol.Plateau Sinica 13, 107–119. Yang, D.R., Li, C.D., Shu, C., Yang, Y.X., 1996. Studies on the Chinese species of the genus Hepialus and their geographical distribution. Acta Entomol. Sinica 39, 413–422. Yang, X., Yang, X., Xue, D., Han, H., 2015. The complete mitochondrial genome of Endoclita signifer Lepidoptera, Hepialidae. Mitochondrial Dna a DNA Mapp. Seq. Anal. 27, 4620–4621. Yi, J.Q., Que, S.Q., Xin, T.R., Xia, B., Zou, Z.W., 2016. Complete mitochondrial genome of Thitarodes puiLepidoptera: Hepialidae. Mitochondrial Dna a DNA Mapp. Seq. Anal. 27, 109–110. Yuan, D.Y., 2003. Tectonic Deformation Features and Space-time Evolution in Northeastern Margin of the Qinghai-Tibetan Plateau since the Late Cenozoic Time. Ph. D. dissertation. Institute of Geology, China Earthquake Administration. Zhang, S., 2013. The population genetics of Ophiocordyceps sinensis and its evolutionary relationship with host insects. Ph. D. dissertation. Shanxi University. Zhang, Y.J., Li, E.W., Wang, C.S., Li, Y.L., Liu, X.Z., 2012. Ophiocordyceps sinensis, the flagship fungus of China, terminology, life strategy and ecology. Mycology 3, 2–10. Zhang, Y.J., Zhang, S., Li, Y.L., Ma, S.L., Wang, C.S., Xiang, M.C., Liu, X.Z., 2014. Phylogeography and evolution of a fungal-insect association on the Tibetan plateau. Mol. Ecol. 23, 5337–5355. Zhong, D.L., Ding, L., 1996. The uplift progress and mechanisms of the Qinghai-Tibetan plateau. Science in China Earth Sci. 26, 289–295. Zou, Z.W., Liu, X., Zhang, G.R., 2010. Revision of taxonomic system of the genus Hepialus Lepidoptera, Hepialidae. currently adopted in China. J. Hunan Univers. Sci. Technology 25, 114–120. Zou, Z.W., Min, Q., Cheng, S.Y., Xin, T.R., Xia, B., 2017. The complete mitochondrial genome of Thitarodes sejilaensis Lepidoptera: Hepialidae., a host insect of Ophiocordyceps sinensis and its implication in taxonomic revision of Hepialus adopted in China. Gene 601, 44–55.

37