Accepted Article

Received Date : 23-Oct-2013 Revised Date : 18-Sep-2014 Accepted Date : 22-Sep-2014 Article type : Original Article

Phylogeography and Evolution of a Fungal-insect Association on the Tibetan Plateau

Yongjie Zhang1,2,3, Shu Zhang2, Yuling Li4, Shaoli Ma4, Chengshu Wang5, Meichun Xiang1, Xin Liu4, Zhiqiang An6, Jianping Xu3,7*, and Xingzhong Liu1*

1

State Key Laboratory of Mycology, Institute of Microbiology, Chinese Academy of

Sciences, Beijing 100101, China;

2

School of Life Sciences, Shanxi University, Taiyuan

030006, China; 3 Department of Biology, McMaster University, Hamilton, Ontario, L8S 4K1, Canada; 4 Institute of Grassland, Qinghai Academy of Animal & Veterinary Sciences, Xining 810016, China;

5

Institute of Plant Physiology and Ecology, Shanghai Institutes for

Biological Sciences, Chinese Academy of Sciences, Shanghai 200032, China;

6

Institute of

Molecular Medicine, University of Texas Health Science Center at Houston, Houston, TX 77030, USA;

7

Laboratory for Conservation and Utilization of Bio-Resources, Yunnan

University, Kunming, 650091, P. R. China.

Keywords: Ophiocordyceps sinensis, ghost moth, parasitoidism, phylogeography, evolution, the Tibetan Plateau

Corresponding authors: Jianping Xu, Department of Biology, McMaster University, Hamilton, Ontario, L8S 4K1, Canada; Phone +1 905 525 9140-27934; email [email protected]; Xingzhong Liu, State Key Laboratory of Mycology, Institute of This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1111/mec.12940 This article is protected by copyright. All rights reserved.

Accepted Article

Microbiology, Chinese Academy of Sciences, Beijing 100101, China; Phone +86 10 6480 7505; email [email protected]

Running title: O. sinensis-ghost moth evolution

Abstract Parasitoidism refers to a major form of inter-species interactions where parasitoids sterilize and/or kill their hosts typically before hosts reach reproductive age. However, relatively little is known about the evolutionary dynamics of parasitoidism. Here we investigate the spatial patterns of genetic variation of Chinese cordyceps, including both the parasitoidal fungus Ophiocordyceps sinensis and its host insects. We sampled broadly from alpine regions on the Tibetan Plateau and obtained sequences on seven fungal and three insect DNA fragments from each of the 125 samples. Seven and five divergent lineages/cryptic species were identified within the fungus and host insects respectively. Our analyses suggested that O. sinensis and host insects originated at similar geographic regions in southern Tibet/Yunnan, followed by range expansion to their current distributions. Cophylogenetic analyses revealed a complex evolutionary relationship between O. sinensis and its host insects. Significant congruence was found between host and parasite phylogenies and the time estimates of divergence were similar, raising the possibility of the occurrence of cospeciation events, but the incongruences suggested that host shifts were also prevalent. Interestingly, one fungal genotype was broadly distributed, consistent with recent gene flow. In contrast, the high-frequency insect genotypes showed limited geographic distributions. The dominant genotypes from both the fungus and the insect hosts may represent ideal materials from which to develop artificial cultivation of this important Chinese traditional medicine. Our results demonstrate that both historical and contemporary events have played important

This article is protected by copyright. All rights reserved.

Accepted Article

roles in the phylogeography and evolution of the O. sinensis-ghost moth parasitoidism on the Tibetan Plateau.

Introduction Interactions between organisms sculpt a myriad of phenomena during the evolution of life (Thompson 2009). At present, evolutionary dynamics between interacting partners have been documented for multiple categories of species interactions, including mutualism (Six & Paine 1999; Aanen et al. 2002; Schultz & Brady 2008), parasitism (Lively 1989; Hafner et al. 1994; Lively & Dybdahl 2000; Desdevises et al. 2002), and predation (Brodie et al. 2002). Although these studies do not show evidence for extensive cospeciation, cophylogenetic signals are often detected (de Vienne et al. 2013). As a major form of biological interaction, parasitoidism differs from typical parasitism in that parasitoids usually sterilize or kill their hosts before the hosts reach their reproductive age (Godfray 1994). Although parasitoids occur in about as many classes of organisms as do parasites, unlike obligate parasites, parasitoids can live independently at certain stage(s) of their life cycle and do not entirely dependent on their hosts for survival. Common examples of parasitoidism are found between invertebrates, such as between the jewel wasp Ampulex compressa and cockroaches and between the braconid wasp Cotesia congregatus and tomato or tobacco hornworms. In addition, some parasitoids have shown capable of regulating the population densities of their hosts and have been used as bio-control agents (Dahlsten et al. 1998). Parasitoids may also manipulate host’s behavior, and thus a close evolutionary relationship between parasitoids and hosts is expected (Libersat et al. 2009; Libersat & Gal 2013). There is currently little knowledge, however, on the population genetics and evolution for interacting partners involved in parasitoidism.

This article is protected by copyright. All rights reserved.

Accepted Article

Fungi and insects are two of the largest groups of organisms on our planet. Fungusinsect interactions range from agonistic to mutualistic (Vega & Blackwell 2005), but the evolutionary dynamics between fungi and insects are mostly studied for mutualistic systems, such as fungus-growing termites and their mutualistic fungal symbionts (Aanen et al. 2002), attine ants and the fungi that they cultivate in their nests (Currie et al. 2003), and Dendroctonus bark beetles and the mycangial fungi that they both consume and disperse (Six & Paine 1999). Among the described fungi, more than 700 species in 90 genera have been recognized as insect pathogens (Roberts & Humber 1981). Some of these fungi cause aberrant behavior of the hosts that can lead to host sterilization and/or host deaths (Libersat et al. 2009; Pontoppidan et al. 2009). One prominent example of the fungus-insect parasitoidism is the Chinese caterpillar fungus Ophiocordyceps sinensis (syn. Cordyceps sinensis) and its insect hosts, the ghost moths that belong to the family Hepialidae (Order Lepidoptera, Class Insecta).

The Chinese caterpillar fungus–ghost moth association is among the best-known entomopathogenic fungus–insect associations (Vilcinskas 2011; Zhang et al. 2012b). Both the fungus and its host insects are endemic to alpine regions along the eastern and southern parts of the Tibetan Plateau (Winkler 2008). The product resulting from parasitoidism of O. sinensis on ghost moths, commonly known as the “Chinese cordyceps” (Zhang et al. 2012b) and “Himalayan Viagra” (Stone 2008), is a renowned traditional Chinese medicine that has shown a diversity of pharmaceutical functions (Zhu et al. 1998a, b; Buenz et al. 2005; Chen et al. 2013; Lo et al. 2013; Yan et al. 2014). Due to huge market demand, our inability to artificially cultivate the association, and its high economic value (currently worth more than its weight in gold), natural populations of Chinese cordyceps are heavily exploited, leading to decreased production and resulted in O. sinensis on the Endangered Species List sanctioned

This article is protected by copyright. All rights reserved.

Accepted Article

by the State Council of the People’s Republic of China (Stone 2008; Zhang et al. 2012b). A comprehensive understanding of the evolutionary relationships between O. sinensis and its host insects will thus provide fundamental knowledge for not only parasitoidism interactions but also for the artificial cultivation of this prominent Chinese medicine.

Though much of the life cycles of the fungi and the host insect remain unknown in nature, investigations have identified that O. sinensis is highly dependent on ghost moths for survival and sexual reproduction. Indeed, the fungal reproductive cycle shows a high degree of synchrony with the life cycle of its host insects (Fig. 1). For example, the fungal fruiting season and the release of ascospores from a stroma (i.e., a dense structural tissue that emerges from a sclerotium and produces fruiting bodies) coincides with the timing when host insects shed cuticles and become vulnerable to fungal infections (Yang et al. 1989). Host insects infected by O. sinensis move from soil tunnels to sub-soil surface with their heads upward, which is advantageous for fungal stroma development from insect heads and for fungal ascospore dispersal (Zhang et al. 2012b). O. sinensis also turns infected caterpillar into a sclerotium (i.e., the mummified insect body) that could help the fungus survive environmental extremes on the Tibetan Plateau during the cold months from November to the following April. However, the fungus can grow on artificial media and a recent investigation detected O. sinensis in plant roots and rhizophere based on PCR assays (Zhong et al, 2014), suggesting additional factors other than the host insect may influence the survival and reproduction of O. sinensis in nature.

Recent research identified significant genetic variation among strains of O. sinensis (Chen et al. 1999; Kinjo & Zang 2001; Liang et al. 2005, 2008; Zhang et al. 2009). These studies have used markers such as RAPD, ISSR, or sequences at the ITS gene and their

This article is protected by copyright. All rights reserved.

Accepted Article

strains were from relatively limited geographic regions. Wang & Yao (2011) summarized the proposed morphological species of host insects of the fungus reported in the literature up to 2011 while Quan et al (2014a) used mainly the COI sequences to investigate the genetic diversity of host insects. While informative at certain levels, the above studies based on RAPD, ISSR, or single copy gene markers and relatively limited samples are insufficient to distinguish strains and cryptic species. Indeed, understanding the evolutionary dynamics between these two partners has been hampered by the difficulty in collecting samples, the uncertainties in host insect identification based merely on larval characteristics, and a lack of polymorphic molecular markers for analysis. Furthermore, although some researchers have compared the phylogenies of O. sinensis and host insects (Chen et al. 1998; Quan et al. 2014b), cophylogenic and spatial genetic patterns between O. sinensis and its host insects have not been critically evaluated.

In this study, we sampled from almost all known distribution areas of the Chinese cordyceps and employed multiple fragment analyses to investigate the evolutionary histories of the fungus and host insects and their evolutionary relationships. In view of their special distribution on alpine regions of the Tibetan Plateau, we also investigated their spatial genetic patterns. We hypothesize that O. sinensis and its host insects have a strong cophylogenetic signal, and that both have pronounced spatial genetic structures. Cophylogenetic patterns were evaluated using both distance- and event-based approaches. Spatial genetic patterns were evaluated by Mantel test, partial Mantel test and spatial autocorrelation analyses. Our results provide valuable information for future artificial cultivation of O. sinensis on host insect larvae, as well as for the conservation and sustainable utilization of this medicinally important resource.

This article is protected by copyright. All rights reserved.

Accepted Article

Materials and methods Sample source Chinese cordyceps samples were collected during harvesting seasons (May and June) on the Tibetan Plateau between 2005 and 2012. After being dug out from the soil, Chinese cordyceps samples were placed into zip plastic bags, transported to the laboratory in a portable refrigerator and then stored at 4°C before use. Pure O. sinensis cultures were isolated from fresh Chinese cordyceps samples within two weeks after sample collection, following the protocols described in Zhang (2012). Our current sampling spanned the known distribution range of Chinese cordyceps, including southwestern China (Xizang Autonomous Region, Qinghai Province, Sichuan Province, Yunnan Province, and Gansu Province), Nepal, and India. A total of 125 Chinese cordyceps samples were included in this study (Supplementary Table S1). They originated from 114 different localities in 81 different counties.

DNA extraction, PCR amplification and sequencing Genomic DNA was extracted with the CTAB method (Zhang et al. 2010) from either the Chinese cordyceps (stroma or sclerotium sections) or O. sinensis pure cultures. The genomic DNA samples were used as templates to amplify seven nuclear fragments of O. sinensis. These gene fragments included the nrDNA ITS region, a mating type gene (MAT12-1), a serine protease-encoding gene (csp1), and four random DNA fragments (OSRC14, OSRC17, OSRC27, and OSRC32) (Zhang et al. 2008, 2009, 2011, 2012a). The four random DNA fragments have either unknown function (OSRC14), or encode hypothetical proteins (OSRC17, OSRC27) or a WD-repeat [also called the beta-transducin repeat, a short structural motif of approximately 40 amino acids, often terminating in a tryptophan-aspartic acid (W-D) dipeptide] containing protein (OSRC32) (Zhang et al. 2012a). Blast analyses against the

This article is protected by copyright. All rights reserved.

Accepted Article

genome of O. sinensis (Hu et al. 2013) indicated that MAT1-2-1, csp1 and the four random fragments are all single copy genes. Although nrDNA exists as variable numbers of tandem repeats, direct sequencing of PCR products of the nrDNA ITS region indicated no heterogeneity (Zhang et al. 2009). Three host insect genes were amplified directly from the total DNA extracted from each of the 125 Chinese cordyceps (the sclerotium section). The three genes include two from the mitochondrial genome (COI and cytb) and one from the nucleus (the wingless gene, wg). Primer sequences as well as annealing temperatures are listed in Supplementary Table S2. PCR products were sequenced using an ABI 3730 XL DNA sequencer with BigDye 3.1 Terminators (Applied Biosystems, Carlsbad, CA, USA). In case of poor sequencing quality resulting from direct sequencing of PCR products or singleton alleles, representative PCR fragments were cloned and individual cloned fragment was then sequenced again. Our manual checking of chromatograms and the cloning and resequencing to confirm singleton alleles minimize the potential PCR and sequencing errors.

Phylogenetic analyses DNA sequences for each locus were aligned using Muscle 3.8 (Edgar 2004) and edited visually in MEGA 6.0 (Tamura et al. 2013). All regions were included in the initial alignment. However, the first intron in csp1 was excluded from subsequent analysis because of ambiguous alignment. Kimura two-parameter (K2P) genetic distances were calculated using MEGA 6.0 (Tamura et al. 2013) with gaps treated as “complete deletion”. Samples were assigned to alleles (i.e., unique DNA sequences) or MLHs (multiple locus haplotypes, i.e., unique combined sequences) with the program DnaSP 5.10 (Librado & Rozas 2009) with information from all indels included.

This article is protected by copyright. All rights reserved.

Accepted Article

We performed partition homogeneity test (PHT) to detect potential conflicts between gene trees of different loci using PAUP version 4.0b10 (Sinauer Associates, Sunderland, MA, USA). PHT used only informative characters and simple stepwise-addition heuristic searches with 1,000 replicates (TBR; maxtrees = 500). A P-value lower than 0.001 indicates statistically significant differences in tree topologies (Cunningham 1997a, 1997b).

Three tree-building methods were used, namely maximum parsimony (MP), Bayesian inference (BI), and maximum likelihood (ML). Only one representative sequence of each MLH was used, and redundant sequences were removed from the combined data sets using TOPALI 2.5 (Milne et al. 2009). Parsimony phylogenetic analyses were conducted using PAUP version 4.0b10 (Sinauer Associates, Sunderland, MA, USA). Gaps were treated as a fifth state character in heuristic searches that were conducted following 1,000 replicates of random stepwise addition and tree bisection–reconnection for branch-swapping. To allow each locus to contribute equally, a weighted MP analysis was performed with characters weighted inversely proportional to the total number of phylogenetically informative sites per locus. When multiple MP trees were produced, the tree chosen for display in a figure was the one determined to be most likely using substitution models identified by PartitionFinder 1.1.1 (Lanfear et al. 2012). Branch support was estimated by performing 1,000 bootstrap replicates with a heuristic search consisting of five random-addition replicates for each bootstrap replicate.

To conduct Bayesian phylogenetic analyses with MrBayes 3.2.2 (Ronquist et al. 2012), we used the program PartitionFinder 1.1.1 (Lanfear et al. 2012) to evaluate the best partitioning scheme and to determine substitution models under the Bayesian Information Criteria (Supplementary Table S3). Two simultaneous MCMC runs with four chains each (3

This article is protected by copyright. All rights reserved.

Accepted Article

hot and 1 cold) were performed for 10 million generations and sampled every 1,000 generations discarding a burn-in of 20%. For partitioned data sets, model parameters across different partitions were unlinked, and the overall rate of evolution among partitions was allowed to vary. Tracer 1.6 (Rambaut & Drummond 2009) was used to ensure the Markov chains had reached stationary states by plotting log-likelihood score and all other model parameter values against the generation number. In addition, we examined the effective sample size (ESS) values (i.e., ESS > 200), the average standard deviation of split frequencies (ASDSF) values (i.e., ASDSF 50% (Supplementary Fig. S3). The ancestral reconstruction of the node representing the seventh clade (Clade I) is ambiguous and suggests five possible ancestral ranges (i.e., G, west Sichuan; I; H, north Tibet; C, north Qinghai; and L, Himalaya).

For host insects, the ancestral range of root node was suggested as area I when using BI tree, but areas K, L, A were also suggested when using MP or ML trees (Supplementary Fig. S4). For the individual host insect clades, Clade D was estimated to have originated in

This article is protected by copyright. All rights reserved.

Accepted Article

area I; Clades A, B, C, and E were estimated most likely at areas G or J (North Yunnan), B (south Gansu), A, and J, respectively.

Spatial genetic pattern The Mantel tests indicated significant IBD within the fungus according to the horizontal geographical distance based on longitudinal–latitudinal coordinates (Pearson correlation coefficient r = 0.182, P = 0.001; Supplementary Table S6). Further analysis found this correlation existed mainly along a latitudinal gradient (r = 0.230, P = 0.001), but not along a longitudinal gradient (r = 0.062, P = 0.122). For the host insects, significant correlations existed along both the latitudinal (r = 0.219, P = 0.001) and longitudinal (r = 0.122, P = 0.001) gradients (Supplementary Table S6). We found no significant effect of altitudinal gradients on genetic distances for either the fungus or the insects (r = 0.048–0.119, P = 0.176–0.320). The results of the partial Mantel tests showed that the relationship revealed by Mantel test were stable even when controlled for genetic differentiation of one partner in the O. sinensis-host insect association (Table 3).

To further evaluate the effects of geographic distance on the genetic signatures obtained in the present study, we evaluated and found significant positive spatial autocorrelation among sampling sites for both O. sinensis and host insects based on horizontal distances no matter what distance matrices (longitudinal–latitudinal distance, latitudinal distance alone, or longitudinal distance alone) were used (Fig. 5; Supplementary Fig. S5). This correlation was strongest at the shortest distances followed by a decrease with increasing distance, suggesting that IBD is strongest among neighboring populations and that the positive relationship between genetic and geographic distance deteriorates with increasing distance among sites. For O. sinensis, spatial autocorrelation coefficient was significant until

This article is protected by copyright. All rights reserved.

Accepted Article

1,200 km (P = 0.001–0.017) based on longitudinal–latitudinal distance, 900 km (P = 0.001– 0.002) based on latitudinal distance, and 1,100 km (P = 0.001–0.029) based on longitudinal distance although there was weak P-value support at 0–400 km (P = 0.057) and 0–900 km (P = 0.054) of longitudinal distance (Fig. 5). For host insects, significant positive spatial autocorrelation was observed for distance intervals up to 1,500 km (P = 0.001–0.023) based on longitudinal–latitudinal distance, 1,000 km (P = 0.001–0.010) based on latitudinal distance, and 1,500 km (P = 0.001–0.032) based on longitudinal distance (Supplementary Fig. S5). Along the altitudinal gradient, neither global significant spatial autocorrelation nor a pattern of decreasing autocorrelation was observed for the fungus and host insects. Interestingly, a significant positive autocorrelation was detected for host insects at certain distance classes, namely at 0–200 m (P = 0.042), 0–500 m (P = 0.014), or 0–650 m (P = 0.022).

Cophylogenetic analysis between O. sinensis and host insects Based on the tanglegram between O. sinensis and host insect MLH phylogenies (Fig. 6), multiple O. sinensis MLHs were found associated with a single host insect MLH, and multiple host insect MLHs were associated with a single O. sinensis MLH. In total, 11 of the 53 fungal MLHs and 8 of the 82 insect MLHs were each associated with more than one host insect and fungal MLH, respectively. The most dominant fungal MLH, MLH 1 represented by GS09-111, was associated with 26 insect MLHs. Similarly, the three most dominant insect MLHs, MLHs 3, 6 and 46 represented by GS09-131, GS09-229 and XZ07-133 respectively, were associated with two, seven, and three fungal MLHs, respectively.

While certain fungal MLHs seemed to be generalists that infected multiple host MLHs and certain host insect MLHs were susceptible to infections by multiple fungal MLHs,

This article is protected by copyright. All rights reserved.

Accepted Article

there was also strong evidence for some specificity between O. sinensis and host insect at both the clade level and individual MLH level. For example, 42 out of the 53 (79.2%) fungal MLHs were found to be associated with only one host MLH (Fig. 6). In addition, fungal MLHs in Clades V and VI were found to be associated only with host MLHs in Clade C. Fungal MLHs in Clades I, II, III, IV, and VII were predominantly associated with host MLHs in Clades A, A, D, D, and B respectively. From the viewpoint of host insects, Clade B was predominantly associated with fungal Clade VII; Clade D with fungal Clades III and IV; Clades A with fungal Clades I and II, and Clade C with fungal Clades V & VI.

The distance-based approach ParaFit indicated that the global congruence between O. sinensis and host insect trees was highly significant (P = 0.001). Among the 97 links between the 53 O. sinensis MLHs and the 82 host insect MLHs, 58 (60%) were statistically significant (P < 0.05; Supplementary Table S7). Separate analyses of each clade in the fungal or host insect phylogenetic trees revealed significant global congruence for every fungal clade and associated insect MLHs, and for four of the five insect clades (except Clade E) and associated fungal MLHs (Supplementary Table S7; Fig. 6). All fungus-insect links were statistically significant (P < 0.05) in fungal Clades IV–VII and corresponding insect MLHs, and in the insect Clade B and corresponding fungal MLHs. For the other fungal and insect clades, 27– 80% of corresponding fungus-insect links were consistent with global congruence (Supplementary Table S7).

Scenarios reconstructed by the event-based approach Jane 4 varied according to the cost scheme employed with 14–28 cospeciations (or codiversification, referring to the simultaneous fungal and host insect divergence) retrieved (Table 4). The number of inferred “failure to diverge” event, which is 44, was stable across different cost schemes, and the

This article is protected by copyright. All rights reserved.

Accepted Article

sorting event (i.e., a fungal MLH that lives on a host species remains on only one of the resulting species after a host divergence) was always the most frequent even when a higher cost was assigned to sorting event. Many host-switching events (always considered subsequent to duplication in Jane) were also deduced while duplication events (i.e., independent fungus divergence) were almost always the least frequent evolutionary event. All Jane generated scenarios indicated significant congruence of phylogenies since zero percent of random sample solutions returned a lower cost than the original problem solution (Supplementary Fig. 6a). Because other evolutionary events in total were more common than cospeciation events, cospeciation might be more costly and we thus tested whether a higher cost assigned to cospeciation would influence our results. However, our test showed that a two-fold cost for cospeciation still returned 14 cospeciation events. Together, these results provided a statistically significant support (P < 0.001) for non-random patterns of the two phylogenies.

Using the adaptive cost evaluation approach of CoRe-PA, 102 Pareto-optimal reconstructions (quality value: 0.01–0.99) were obtained from 10,000 random cost models. We investigated the first 10 preferred reconstructions and found they were surprisingly similar in their cost models (Table 5). These 10 reconstructions had quality values of 0.01– 0.039 and contained 15–23 cospeciation, 20–28 duplication, 5–10 host switching and 74–106 sorting events (optimal costs = cospeciation: 0.19–0.27, duplication: 0.17–0.23, host switch: 0.49–0.57, and sorting: 0.04–0.06). The statistical significance test with 1,000 randomized host-parasite associations resulted in 98.4% of the solutions having less than 20 cospeciations which is the average number of cospeciation events from those 10 reconstructions (Supplementary Fig. 6b). The results further indicate an overall significant amount of congruence between the phylogenetic trees of O. sinensis and its host insects.

This article is protected by copyright. All rights reserved.

Accepted Article

The significant correlation between O. sinensis and host insects was also corroborated by a partial Mantel test where geography was set as a constant (Table 3). This correlation remained stable when either horizontal or vertical geographic distance was held constant.

Local adaptation and density dependence Our above genetic analyses suggested that one fungal MLH was capable of infecting multiple host insect MLHs. To experimentally confirm this observation and test the hypothesis that fungal strains are more infective to hosts from the same geographic region, we collected O. sinensis isolates and corresponding host insects from two distantly separated regions and performed a laboratory cross-inoculation experiment. Each of the two fungal isolates had the same genotype as GS09-111 within Clade I (i.e., the most dominant fungal MLH identified in this study). One of the two host insects had a genotype very similar to QH09-11 and QH09-151 (with 1–2 nucleotide difference) within insect Clade A, and the other to GS09-229 and GS09-311 (with 1–2 nucleotide difference) within insect Clade C. Our results indicated that both isolates were capable of infecting host insects from both regions. Furthermore, as expected, the infection rates of sympatric fungus–insect combinations (mean = 23.1%; n = 2) were about 10 times higher than that of allopatric combinations (mean = 2.3%; n = 2) (Fig. 7a), indicating a significantly greater infectivity of fungal parasitoids for local host insects. The close ecological dependence of the fungi on host insects at a local scale was also supported by the positive correlation between the numbers of naturally occurring O. sinensis fruiting bodies and the density of host insect larvae in soil (r = 0.99; P = 3.49E-11; two-tailed t-test) (Fig. 7b).

This article is protected by copyright. All rights reserved.

Accepted Article

Discussion The Tibetan Plateau is the highest and largest plateau in the world and arguably the most prominent topological feature on Earth. Various regional fauna and flora on the Tibetan Plateau have been reported to have originated before the major uplift of the Tibetan Plateau in the last 10 million years, and diversified during the uplift that created a diversity of ecological niches and gene flow barriers (see reviews in Yang et al. 2009 and Favre et al. 2014). However, little is known about how the uplift of the Tibetan Plateau affected the evolution of microorganisms and species interactions involving microorganisms. Here we used a prominent Tibetan Plateau endemic, the Chinese caterpillar fungus O. sinensis and its ghost moth hosts, as a model system to investigate the potential effects of the geological changes on microbial evolution and species interaction on the Tibetan Plateau.

Our analyses identified significant divergence within both the fungus O. sinensis and its ghost moth hosts, with seven and five statistically well-supported lineages respectively. For example, the ITS sequence divergence within O. sinensis is about four times of that within another caterpillar fungus Cordyceps militaris (Wang et al. 2008). Our results suggest the taxonomy of both the fungus and the host insect needs significant revision. At present, all the Chinese cordyceps fungal isolates are classified into one morphological species (i.e., Ophiocordyceps sinensis) while about 50 morphological species have been proposed for the host insect. According to Wang & Yao (2011), the host insects mainly belong to the genus Thitarodes in the family Hepialidae although several are found in other genera in Hepialidae. While COI has been widely accepted as a barcode of animals, few COI sequences are available for host insects of O. sinensis especially those with high-quality type specimens. Thus, in the absence of corresponding morphological and molecular sequence information for well-curated and representative specimens, we prefer to describe the observed

This article is protected by copyright. All rights reserved.

Accepted Article

genotypes or genotype groups as multilocus haplotypes or clades. However, we believe our study provides robust framework from which to analyze morphology and speciation for both groups of organisms.

Our study revealed a significant global congruence between the fungal and host insect phylogenies. All four methods, including the distance-based ParaFit, event-based Jane 4 and CoRe-PA, and the partial Mantel test showed an overall consistent pattern. Multiple cospeciation events were identified by both Jane 4 and CoRe-PA. In addition, time estimate for the origin of the fungus based on the ITS dataset was very close to that of host insects based on the COI data set. The origination time estimate of the fungus based on the 5-locus dataset, however, differed from that based on the ITS dataset, but was similar to the time estimate of host insects reported by Quan et al (2014a). Using a slower substitution rate, Quan et al (2014a) revealed that host insects first diverged around 4.1 Mya. In the absence of fossil evidence and/or more robust molecular clocks for these two organisms, the above time estimates must be considered tentative. However, our current analyses can’t reject the hypothesis that the fungus and host insects might have originated at about the same time. In addition, the biogeographic inference suggested that the fungus and host insects have likely originated in similar areas, around south Tibet/Yunnan, followed by dispersion and diversification to other areas.

Aside from similar origins and cospeciation events between the two partners, the event-based analyses also identified a large number of other events like duplication and host switch. Given the nature of parasitoidism (sterilization and/or killing the hosts), frequent host shifts were expected. Because the fungal parasitoid can reproduce both asexually and sexually and generate a large number of progeny that can spread quickly, in the long-term,

This article is protected by copyright. All rights reserved.

Accepted Article

this fungal-ghost moth association might not be evolutionary stable. However, other factors could contribute to the relative stability of this parasitoidism in the short and intermediate future, as it has been for the last several million years. For example, the changing topography of the Tibetan Plateau and the differences in dispersal barriers between the fungus and host insects have contributed significantly to the evolutionary dynamics between these two interacting partners. Additional analyses of the spatial and temporal dynamics at the local level could reveal not only the relative frequencies of the different evolutionary events within individual ecological niches but also the potential long-term evolutionary trends between the fungus and host insects.

While the observed genetic isolation by geographic distance is consistent with the significant effects of historical geological events during the uplift of the Tibetan Plateau (see reviews in Yang et al. 2009 and Favre et al. 2014) for both the fungus and host insects, the genotype data also suggest contemporary gene flows among regions. Specifically, the sharing of multilocus genotypes among fungal isolates from diverge geographic locations was consistent with frequent, contemporary gene flow. The gene flow was likely due to longdistance spore dispersals by wind, and indeed, wind-assisted dispersals were often suggested as the main mechanism for dispersal in many fungi (Aylor 1986; Allen et al. 1989). However, human activities such as collecting and trading the Chinese cordyceps could also have facilitated the dispersals. In contrast to the abundant evidence for long-distance dispersals in the fungus, long-distance dispersals for host insects seemed very limited as most shared genotypes were between strains from adjacent geographic regions. Similar results were reported by Yang et al. (1996) who found that most host insect species of O. sinensis had a very narrow distribution on the Tibetan Plateau and host insect species might vary among different mountain ranges and even from different sides and habitats of the same mountain.

This article is protected by copyright. All rights reserved.

Accepted Article

Our results are also similar to what have been found in other insects on the Tibetan Plateau (Yan et al. 2006; Schmidt et al. 2012).

Aside from the geological and climatic factors, other factors have likely played important roles in the biogeography and evolution of the O. sinensis–ghost moth interactions. For example, the larvae of host insects live underground and feed on the roots of plants in alpine meadows on the Tibetan Plateau. Plants in at least 19 angiosperm families are suspected food sources for the ghost moth larvae (Zhu et al. 2004). In addition, a very recent study by nested PCR indicated that O. sinensis is present in the plant roots during its anamorphic life cycle (Zhong et al. 2014). However, it remains unknown if different fungal and host insect lineages prefer the roots of different plants. If such preferences exist, the distributions of plants can have a significant impact on the fungal and/or host insect communities. Close ecological and evolutionary relationships among interacting partners have been reported for tripartite and quadripartite interactions (Weiblen & Bush 2002; Currie et al. 2003; Lopez-Vaamonde et al. 2005; Wilson et al. 2012).

Understanding the evolutionary and ecological relationships among the plant communities and the fungal and host insect genotypes will greatly facilitate our efforts to raise the appropriate host insects in captivity and domesticate the fungal–insect complex. Furthermore, our knowledge about the genotype distributions of both O. sinensis and the host insects as well as their patterns of phylogenetic relationships lay a solid foundation towards artificially cultivating the O. sinensis–host insect complex. For example, the fungal MLH 1 accounts for 37.6% of the analyzed individuals, and this MLH was found in almost all distribution areas from south to north and from west to east on the Tibetan Plateau. Crossinoculation experiments and the fungus-insect phylogenetic tanglegram both showed that the

This article is protected by copyright. All rights reserved.

Accepted Article

fungal MLH 1 could infect different MLHs of host insects. These data suggest the fungal MLH 1 should be an excellent fungal candidate for artificial cultivation of Chinese cordyceps. For host insects, MLHs 3 (represented by GS09-131), 6 (represented by GS09-229), and 46 (represented by XZ07-133) were numerically the most dominant genotypes. Even though these three insect MLHs were geographically relatively limited, each of the three MLHs could be infected by more than one fungal MLH. Of special interest are the observed associations between insect MLHs 3 and 46 and the fungal MLH 1. Even though we did not find an individual containing the insect MLH 6 and the fungal MLH 1, their association is highly possible in nature because one of the host insects used in the cross-infection experiment differed from MLH 6 by only one nucleotide at the COI gene fragment. Taken together, the results suggest that all three insect MLHs could be used for artificial cultivation practices with fungal strains with the MLH 1 genotype. In addition, since the three insect MLHs were from different geographic regions, artificial cultivations in different regions could adopt different fungal-ghost moth genotype combinations. For example, in north Tibet, the preferred genotype combination might be fungal MLH 1 and insect MLH 46. However, in south Gansu, middle and south Qinghai, the preferred combination might be fungal MLH 1 and insect MLH 3.

In conclusion, our study suggests that the evolution and diversification of O. sinensis and its host insects have been impacted by each other and by abiotic factors such as the uplift of the Tibetan Plateau and climatic events. Our analyses identified the large-scale spatial patterns of genetic variation and inferred its tentative temporal diversification patterns. The phylogeographic and evolutionary patterns identified here provide a broad framework from which to investigate the effects of other abiotic (e.g., river systems) and biotic (e.g., plant communities) factors on the evolution of O. sinensis and its host insects. Our study also

This article is protected by copyright. All rights reserved.

Accepted Article

suggested candidate genotypes from which to develop artificial cultivations of the Chinese cordyceps. A method of successful cultivation of the O. sinensis–ghost moth association should alleviate the harvesting pressure on this natural resource and provide a reliable supply of this important Chinese traditional medicine to a wider consumer community.

Acknowledgements This work was supported by the National Plans of Sciences and Technology (grants number 2007BAI32B04), the National Outstanding Youth Foundation (30625001), the cooperative project

between

Guangdong

Province

and

the

Chinese

Academy

of

Sciences

(2009B091300015), the National Science Foundation of China (81102759, 31140013), the Yunnan Provincial Department of Science and Technology (2010CI106), and the Specialized Research Fund for the Doctoral Program of Higher Education (20101401120007). We thank all the people who helped us collect Chinese cordyceps samples: R. Chettry from Nepal, B. Dhar from India, and the following people from China: W. Luo, S. Ma, P. Cai, S. Wei, M. Wang, Z. Luo, X. Tian, H. Wen, L. Guo, L. Xu, F. Li, B. Sun, L. He, X. Jiang, F. Bai, Y. Hao, W. Gong, and Z. Hao. We also thank N. Wahlberg (Finland) for his kind presentation of genomic DNA of Hepialus humuli and Z. Ren, X.-Y. Zhang, X.-L. Zhang, J. Sun, and Y. Li for their help and suggestions in data analysis. The authors thank Profs. J.W. Taylor (University of California, Berkeley), B. Jaffee (University of California, Davis), J.S. Quinn, B.J. Evans, and B. Golding (McMaster University) for their critical comments.

References Aanen DK, Eggleton P, Rouland-Lefèvre C et al. (2002) The evolution of fungus-growing termites and their mutualistic fungal symbionts. Proceedings of the National Academy of Sciences of the United States of America, 99, 14887–14892.

This article is protected by copyright. All rights reserved.

Accepted Article

Allen MF, Hipps LE, Wooldridge GL (1989) Wind dispersal and subsequent establishment of VA mycorrhizal fungi across a successional arid landscape. Landscape Ecology, 2, 165-171. Aylor DE (1986) A framework for examining inter-regional aerial transport of fungal spores. Agricultural and Forest Meteorology, 38, 263–288. Berbee ML, Taylor JW (2010) Dating the molecular clock in fungi – how close are we? Fungal Biology Reviews, 24, 1–16. Bonnet E, Van de Peer Y (2002) ZT: a software tool for simple and partial Mantel tests. Journal of Statistical Software, 7, 1–12. Bouckaert R, Heled J, Kuhnert D et al. (2014). BEAST 2: A software platform for bayesian evolutionary analysis. PLoS Computational Biology, 10, e1003537. Brodie EDJ, Ridenhour BJ, Brodie EDI (2002) The evolutionary response of predators to dangerous prey: hotspots and coldspots in the geographic mosaic of coevolution between garter snakes and newts. Evolution, 56, 2067–2082. Brower AV (1994) Rapid morphological radiation and convergence among races of the butterfly Heliconius erato inferred from patterns of mitochondrial DNA evolution. Proceedings of the National Academy of Sciences of the United States of America, 91, 6491–6495. Buenz E, Bauer B, Osmundson T, Motley T (2005) The traditional Chinese medicine Cordyceps

sinensis

and

its effects

on

apoptotic

homeostasis. Journal of

Ethnopharmacology, 96, 19–29. Charleston

MA

(2011)

TreeMap

3b.

Available

from

https://sites.google.com/site/cophylogeny. Chen PX, Wang S, Nie S, Marcone M (2013) Properties of Cordyceps sinensis: A review. Journal of Functional Foods, 5, 550–569.

This article is protected by copyright. All rights reserved.

Accepted Article

Chen YJ, Yang DR, Yang YX, Zhang YP (1998) Molecular evolution of Cordyceps sinensis from seven regions on the Tibetan Plateau. In: Proceedings of the Third Youth Annual Conference of the China Association for Science and Technology (Materials Science and Engineering) (ed. China Association for Science and Technology), pp 87–89. China Science and Technology Press, Beijing. Chen YJ, Zhang YP, Yang YX, Yang DR (1999) Genetic diversity and taxonomic implication of Cordyceps sinensis as revealed by RAPD markers. Biochemical Genetics, 37, 201–213. Conow C, Fielder D, Ovadia Y, Libeskind-Hadas R (2010) Jane: a new tool for the cophylogeny reconstruction problem. Algorithms for Molecular Biology, 5, 16. Cunningham CW (1997a) Can three incongruence tests predict when data should be combined? Molecular Biology and Evolution, 14, 733–740. Cunningham CW (1997b) Is congruence between data partitions a reliable predicator of phylogenetic accuracy? Empirically testing an iterative procedure for choosing among phylogenetic methods. Systematic Biology, 46, 464–478. Currie CR, Wong B, Stuart AE et al. (2003) Ancient tripartite coevolution in the attine antmicrobe symbiosis. Science, 299, 386–388. Dahlsten D, Rowney D, Copper W et al. (1998) Parasitoid wasp controls blue gum psyllid. California Agriculture, 52, 31–34. Desdevises Y, Morand S, Jousson O, Legendre P (2002) Coevolution between Lamellodiscus (Monogenea: Diplectanidae) and Sparidae (Teleostei): the study of a complex hostparasite system. Evolution, 56, 2459–2471. Edgar RC (2004) MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acid Research, 32, 1792–1797.

This article is protected by copyright. All rights reserved.

Accepted Article

Favre A, Packert M, Pauls SU et al. (2014) The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. Biological Reviews, doi: 10.1111/brv.12107 [Epub ahead of print]. Godfray HCJ (1994) Parasitoids, Behavioral and Evolutionary Ecology. Princeton University Press, Princeton. Hafner MS, Sudman PD, Villablanca FX et al. (1994) Disparate rates of molecular evolution in cospeciating hosts and parasites. Science, 265, 1087–1090. Hu X, ZhangYJ, Xiao GH et al. (2013) Genome survey uncovers the secrets of sex and lifestyle in caterpillar fungus. Chinese Science Bulletin, 58, 2846–2854. Kinjo N, Zang M (2001) Morphological and phylogenetic studies on Cordyceps sinensis distributed in southwestern China. Mycoscience, 42, 567–574. Lanfear R, Calcott B, Ho SY, Guindon S (2012) Partitionfinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Molecular Biology and Evolution, 29, 1695–1701. Legendre P, Desdevises Y, Bazin E (2002) A statistical test for host-parasite coevolution. Systematic Biology, 51, 217–234. Liang HH, Cheng Z, Yang XL et al. (2008) Genetic diversity and structure of Cordyceps sinensis populations from extensive geographical regions in China as revealed by intersimple sequence repeat markers. The Journal of Microbiology, 46, 549–556. Liang HH, Cheng Z, Yang XL et al. (2005) Genetic variation and affinity of Cordyceps sinensis in Qinghai Province based on analysis of morphologic characters and intersimple sequence repeat markers. Chinese Traditional and Herbal Drugs, 36, 1859– 1864. Libersat F, Delago A, Gal R (2009) Manipulation of host behavior by parasitic insects and insect parasites. Annual Review of Entomology, 54, 189–207.

This article is protected by copyright. All rights reserved.

Accepted Article

Libersat F, Gal R (2013) What can parasitoid wasps teach us about decision-making in insects? The Journal of Experimental Biology, 216, 47–55. Librado P, Rozas J (2009) DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics, 25, 1451–1452. Lively CM (1989) Adaptation by a parasitic trematode to local populations of its host. Evolution, 43, 1663–1671. Lively CM, Dybdahl MF (2000) Parasite adaptation to locally common host genotypes. Nature, 405, 679–681. Lo HC, Hsieh C, Lin FY, Hsu TH (2013) A systematic review of the mysterious caterpillar fungus Ophiocordyceps sinensis in DongChongXiaCao (冬蟲夏草 Dōng Chóng Xià Cǎo) and related bioactive ingredients. Journal of Traditional and Complementary Medicine, 3, 16–32. Lopez-Vaamonde C, Godfray HCJ, West SA, Hansson C, Cook JM (2005) The evolution of host use and unusual reproductive strategies in Achrysocharoides parasitoid wasps. Journal of Evolutionary Biology, 18, 1029–1041. Meier-Kolthoff JP, Auch AF, Huson DH, Göker M (2007) CopyCat: cophylogenetic analysis tool. Bioinformatics, 23, 898–900. Merkle D, Middendorf M, Wieseke N (2010) A parameter-adaptive dynamic programming approach for inferring cophylogenies. BMC Bioinformatics, 11, S60. Milne I, Lindner D, Bayer M et al. (2009) TOPALi v2: a rich graphical interface for evolutionary analyses of multiple alignments on HPC clusters and multi-core desktops. Bioinformatics, 25, 126–127. Peakall ROD, Smouse PE (2012) GENALEX 6.5: genetic analysis in Excel. Population genetic software for teaching and research -- an update. Bioinformatics, 28, 2537–2539.

This article is protected by copyright. All rights reserved.

Accepted Article

Pontoppidan MB, HimamanW, Hywel-Jones N, Boomsma JJ, Hughes DP (2009) Graveyards on the move: the spatio-temporal distribution of dead Ophiocordyceps-infected ants. PLoS ONE, 4, e4835. Quan QM, Chen LL, Wang X et al. (2014a) Genetic diversity and distribution patterns of host insects of caterpillar fungus Ophiocordyceps sinensis in the Qinghai-Tibet Plateau. PLoS ONE, 9, e92293. Quan QM, Wang QX, Zhou XL et al. (2014b) Comparative phylogenetic relationships and genetic structure of the caterpillar fungus Ophiocordyceps sinensis and its host insects inferred from multiple gene sequences. Journal of Microbiology, 52, 99–105. Rambaut

A,

Drummond

AJ

(2009)

Tracer

version

1.5.

Available

from

http://beast.bio.ed.ac.uk/Tracer. Roberts DW, Humber RA (1981) Entomogenous fungi. In: Biology of conidial fungi (eds Cole G & Kendrick B), pp. 201–236. Academic Press, New York. Ronquist F, Teslenko M, van der Mark P et al. (2012) MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology, 61, 539–542. Schmidt J, Opgenoorth L, Holl S, Bastrop R (2012) Into the Himalayan exile: the phylogeography of the ground beetle Ethira clade supports the Tibetan origin of forestdwelling Himalayan species groups. PLoS ONE, 7, e45482. Schultz TR, Brady SG (2008) Major evolutionary transitions in ant agriculture. Proceedings of the National Academy of Sciences of the United States of America, 105, 5435–5440. Six DL, Paine TD (1999) Phylogenetic comparison of ascomycete mycangial fungi and Dendroctonus bark beetles (Coleoptera: Scolytidae). Annals of the Entomological Society of America, 92, 159–166. Stone R (2008) Last stand for the body snatcher of the Himalayas? Science, 322, 1182.

This article is protected by copyright. All rights reserved.

Accepted Article

Sukumaran J, Holder MT. (2010) DendroPy: a Python library for phylogenetic computing. Bioinformatics, 26, 1569–1571. Sung GH, Poinar Jr GO, Spatafora JW (2008) The oldest fossil evidence of animal parasitism by fungi supports a Cretaceous diversification of fungal–arthropod symbioses. Molecular Phylogenetics and Evolution, 49, 495–502. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S (2013) MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0. Molecular Biology and Evolution, 30, 2725–2729. Taylor JW, Jacobson DJ, Kroken S, et al. (2000) Phylogenetic species recognition and species concepts in fungi. Fungal Genetics and Biology, 31, 21–32. Thompson JN (2009) The coevolving web of life. The American Naturalist, 173, 125–140. Vega FE, Blackwell M (2005) Insect-Fungal Associations: Ecology and Evolution. Oxford University Press, Cary, USA. de Vienne DM, Refregier G, Lopez-Villavicencio M et al. (2013) Cospeciation vs host-shift speciation: methods for testing, evidence from natural associations and relation to coevolution. New Phytologist, 198, 347–385. Vilcinskas A (2011) Insect Biotechnology. Springer, Dordrecht. Wang XL, Yao YJ (2011) Host insect species of Ophiocordyceps sinensis: a review. ZooKeys, 127, 43–59. Weiblen GD, Bush GL (2002) Speciation in fig pollinators and parasites. Molecular Ecology, 11, 1573–1578. Wilson JS, Forister ML, Dyer LA et al. (2012) Host conservatism, host shifts and diversification across three trophic levels in two Neotropical forests. Journal of Evolutionary Biology, 25, 532–546.

This article is protected by copyright. All rights reserved.

Accepted Article

Winkler D (2008) Yartsa Gunbu (Cordyceps sinensis) and the fungal commodification of Tibet’s rural economy. Economic Botany, 62, 291–305. Yan JK, Wang WQ, Wu JY (2014) Recent advances in Cordyceps sinensis polysaccharides: Mycelial fermentation, isolation, structure, and bioactivities: A review. Journal of Functional Foods, 6, 33–47. Yan L, Wang G, Liu CZ (2006) Number of instars and stadium duration of Gynaephora menyuanensis (Lepidoptera: Lymantriidae) from Qinghai-Tibetan Plateau in China. Annals of the Entomological Society of America, 99, 1012–1018. Yang DR, Li CD, Shu C, Yang YX (1996) Studies on Chinese species of the genus Hepialus and their geographical distribution. Acta Entomologica Sinica, 39, 413–422. Yang SJ, Dong HL, Lei FM (2009) Phylogeography of regional fauna on the Tibetan Plateau: A review. Progress in Natural Science, 19, 789–799. Yang YX, Yang DR, Shen FR, Dong DZ (1989) Studies on hepialid larvae for being infected by Chinese “insect herb” fungus (Cordyceps sinensis). Zoological Research, 10, 227– 231. Yu Y, Harris A, He XJ (2013) RASP (Reconstruct Ancestral State in Phylogenies) 2.1 beta. Available from http://mnh.scu.edu.cn/soft/blog/RASP. Zhang S, Zhang YJ, Liu XZ et al. (2011) Cloning and analysis of the MAT1-2-1 gene from the traditional Chinese medicinal fungus Ophiocordyceps sinensis. Fungal Biology, 115, 708–714. Zhang YJ (2012) Biology of the Chinese Caterpillar Fungus Ophiocordyceps sinensis. Science Press, Beijing. Zhang YJ, Bai FR, Zhang S, Liu XZ (2012a) Determining novel molecular markers in the Chinese caterpillar fungus Ophiocordyceps sinensis by screening a shotgun genomic library. Applied Microbiology and Biotechnology, 95, 1243–1251.

This article is protected by copyright. All rights reserved.

Accepted Article

Zhang YJ, Li EW, Wang CS, Li YL, Liu XZ (2012b) Ophiocordyceps sinensis, the flagship fungus of China: terminology, life strategy and ecology. Mycology, 3, 2–10. Zhang YJ, Liu XZ, Wang M (2008) Cloning, expression, and characterization of two novel cuticle-degrading serine proteases from the entomopathogenic fungus Cordyceps sinensis. Research in Microbiology, 159, 462–469. Zhang YJ, Xu LL, Zhang S et al. (2009) Genetic diversity of Ophiocordyceps sinensis, a medicinal fungus endemic to the Tibetan Plateau: Implications for its evolution and conservation. BMC Evolutionary Biology, 9, 290. Zhang YJ, Zhang S, Liu XZ, Wen HA, Wang M (2010) A simple method of genomic DNA extraction suitable for analysis of bulk fungal strains. Letters in Applied Microbiology, 51, 114–118. Zhong X, Peng QY, Li SS et al. (2014) Detection of Ophiocordyceps sinensis in the roots of plants in alpine meadows by nested-touchdown polymerase chain reaction. Fungal Biology, 118, 359–363. Zhu HF, Wang LY, Han HX (2004) Fauna Sinica: Insecta, Lepidoptera Hepialidae Epiplemidae Volume 38. Science Press, Beijing. Zhu JS, Halpern GM, Jones K (1998a) The scientific rediscovery of an ancient Chinese herbal medicine: Cordyceps sinensis Part I. Journal of Alternative and Complementary Medicine, 4, 289–303. Zhu JS, Halpern GM, Jones K (1998b) The scientific rediscovery of a precious ancient Chinese herbal regimen: Cordyceps sinensis Part II. Journal of Alternative and Complementary Medicine, 4, 429–457. Zwickl DJ (2006) Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets under the maximum likelihood criterion. Ph.D. dissertation, The University of Texas at Austin.

This article is protected by copyright. All rights reserved.

Accepted Article

Author contributions Y.J.Z., Y.L.L., J.P.X., and X.Z.L. designed the research; Y.J.Z., S.Z., and S.L.M performed the research; Y.J.Z., J.P.X., M.C.X., and X.Z.L analyzed the data; Y.J.Z., J.P.X., Z.Q.A., C.S.W., and X.Z.L wrote the paper.

Data accessibility DNA sequences have been submitted to GenBank under accession numbers JQ325053– JQ326199 and KM197519–KM197558 (See also Supplementary Table S1). The aligned sequences for all individuals and the generated tree files are available via Dryad Digital Repository (doi:10.5061/dryad.13272). The raw data from the infection trials and the density survey are also available in Dryad (doi:10.5061/dryad.13272).

Figure Legends Figure 1 Life cycles of Ophiocordyceps sinensis and its host insect. The fungus completes its life cycle within one year, while the host insect needs more than two years to complete its life cycle. The fungus overwinters in the form of sclerotium (i.e., mummified insect body) and the host insects in the form of larvae at various instars between November and April (shading areas) each year. When infected by the fungus, the life cycle of the host insect larvae is halted. This occurs when the 4th to 5th instar larvae are shedding cuticles and become susceptible to fungal infection (Yang et al. 1989). Ascospores of the fungus germinate by microcycle conidiation (i.e., the direct formation of conidia without involving mycelial growth) or by producing hyphae. At present, the relative importance of the different morphological forms of the fungus during infection is still not clear.

This article is protected by copyright. All rights reserved.

Accepted Article

Figure 2 O. sinensis MLH identifier of each sampling site. Sampling sites above the county level are indicated on the map by colored triangles. MLH identifiers of O. sinensis are shown for each county (see also Supplementary Table S1). Note that the fungal MLH 1 is the most widely distributed across the sampling areas.

Figure 3 Phylogenetic analyses of O. sinensis based on the concatenated 6-locus dataset. Topologies reconstructed by Bayesian inferences (a), maximum parsimony (b), and maximum likelihood (c), are shown. Nodes supported by ≥ 70% parsimony/likelihood bootstrap values or by ≥0.95 Bayesian posterior probability values are indicated. Same numbers of clades were identified by each phylogenetic approach, but there is topological difference for the clades among the three approaches.

Figure 4 Phylogenetic analyses of host insects based on the concatenated 2-locus data set. Topologies reconstructed by Bayesian inferences (a), maximum parsimony (b), and maximum likelihood (c) are shown. Nodes supported by ≥ 70% parsimony/likelihood bootstrap values or by ≥0.95 Bayesian posterior probability values are indicated. Same numbers of clades were identified by each phylogenetic approach, but there is topological difference for the clades among the three approaches. In the ML tree (4c), the MLHs of clade A did not cluster into one monophyletic clade but shared polytomies with the parental nodes of clades B and E. Thus, we used a rightward-oriented triangle to differentiate clade A from clades B and E. The rightward-oriented triangle of clade C was based on a similar consideration.

This article is protected by copyright. All rights reserved.

Accepted Article

Figure 5 Correlograms indicating the average genetic autocorrelation coefficient (r) as a function of increasing distance class size for O. sinensis. Spatial autocorrelation analyses were performed based on various geographical distances, namely longitudinal–latitudinal distance (a), latitudinal distance (b), longitudinal distance (c), and altitudinal distance (d). Error bars indicate the 95% confidence interval around the observed r values, and red dash marks indicate the 95% confidence interval surrounding the null hypothesis of no spatial autocorrelation (r = 0).

Figure 6 The tanglegram between O. sinensis and host insect MLH phylogenies. Fungal (right) and host insect (left) phylogenies from Bayesian inference were used to generate the tanglegram using TreeMap 3.0. The most dominant fungal MLH 1 was represented by GS09111. The three most dominant insect MLHs 3, 6, and 46 were represented by GS09-131, GS09-229, and XZ07-133, respectively.

Figure 7 Evidence for local adaptation and density dependence. Evidence for local adaptation (a) was verified by cross-inoculation experiments using fungal isolates and insect larvae from two distinct sites. Vertical bars represent one standard error of the mean of infection frequency. Columns with different letters at the top represent statistically significant differences in their frequencies of infection (P < 0.01; one-way ANOVA). Evidence for density dependence (b) was verified by comparing the numbers of naturally occurring O. sinensis fruiting bodies with the density of underground host insect larvae from a 14-year field survey. The numbers of naturally occurring O. sinensis fruiting bodies are represented by the average number of Chinese cordyceps samples collected by experienced collectors within each day across the whole harvesting season each year. The density of host insect

This article is protected by copyright. All rights reserved.

Accepted Article

larvae is represented by the numbers of host insect larvae found within 0.25 m2 of soil, and the values are means of five replicates.

Table 1. Nucleotide variations of DNA fragments among 125 individuals of Chinese cordyceps

Fragment a

Aligned length (actual length)

Variable site b

Average genetic distance, d (SE)

Allele or MLH (Frequenc y) c

Singleton allele or MLH

28 (52)

15

22 (75)

14

18 (62)

6

12 (90)

7

15 (71)

6

13 (90)

4

21 (41)

12

O. sinensis nrDNA ITS MAT1-2-1 csp1 OSRC14 OSRC17 OSRC27 OSRC32 Combined 7locus data set Combined 6locus data set excluding OSRC32 Host insect

543 (531538) 878 (873878) 581 (579581) 590 (557590) 481 (479480) 612 (571598) 641 (611638) 4326 (42254258) 3685 (35993643)

COI

608 (608)

cytb

385 (379385)

wg

411 (411)

Combined 3locus data set

1404 (13981404)

56 (36, 8, 12) 54 (38, 11, 5)

71 (30, 0, 41) 158 (103, 22, 33)

0.0102 (0.0022) 0.0047 (0.0009) 0.0067 (0.0014) 0.0057 (0.0015) 0.0057 (0.0015) 0.0048 (0.0011) 0.0477 (0.0055)

466 (278, 60, 128)

0.0119 (0.0009)

64 (25)

49

308 (175, 38, 95)

0.0061 (0.0006)

53 (47)

38

166 (124, 42, 0) 131 (99, 26, 6) 69 (43, 26, 0)

0.0498 (0.0051) 0.0743 (0.0071) 0.0224 (0.0038)

77 (15)

66

60 (18)

48

28 (46)

18

366 (266, 94, 6)

0.0477 (0.0034)

89 (11)

79

37 (33, 2, 2) 60 (19, 8, 33) 30 (19, 9, 2)

This article is protected by copyright. All rights reserved.

Accepted Article

Combined 21019 235 (167, 0.0385 locus data set 82 (14) 71 (1019) 68, 0) (0.0036) excluding cytb a Sequences of the first intron of csp1 are excluded due to ambiguous alignment. No any region was excluded from other fragments. The fungal combined data set without OSRC32 and the insect combined data set without cytb were used in all following analyses. b

Variable sites consist of substitution sites (phylogenetically informative or

uninformative) and indel (insertion/deletion) sites. The total number of variable sites is given outside parentheses; that of phylogenetically informative sites, phylogenetically informative sites and indel sites are given individually within parentheses. c Numbers

outside parentheses represents the number of alleles or multi-locus

haplotypes (MLH), and those within parentheses are the numbers of individuals represented by the most dominant allele or MLH. An allele refers to any variant of DNA sequence observed at a given locus with indels considered. A MLH refers to a group of alleles at multiple loci.

Table 2. Results of the partition homogeneity tests nrDNA ITS MAT1-2-1 csp1 OSRC14

OSRC17

MAT1-2-1

0.340

csp1

0.339

0.216

OSRC14

0.505

0.573

0.109

OSRC17

0.324

0.509

0.306

0.035

OSRC27

0.328

0.018

0.342

0.008

0.024

OSRC32

0.001

0.003

0.001

0.001

0.001

COI

cytb

cytb

0.001

wg

0.022

0.001

This article is protected by copyright. All rights reserved.

OSRC27

0.001

Accepted Article

Table 3 Partial Mantel test among O. sinensis genetic distance, host insect genetic distance and geographic distances Correlation value Matrices Covariate P-value (r) O. sinensis & host insect Longitude & latitude 0.356 0.001 Latitude alone 0.349 0.001 Longitude alone 0.378 0.001 Altitude 0.557 0.011 O. sinensis & longitudinalHost insects 0.109 0.008 latitudinal distance O. sinensis & latitudinal distance Host insects 0.162 0.001 O. sinensis & longitudinal Host insects 0.017 0.355 distance O. sinensis & altitudinal distance Host insects -0.022 0.484 Host insects & longitudinalO. sinensis 0.164 0.001 latitudinal distance Host insects & latitudinal O. sinensis 0.146 0.001 distance Host insects & longitudinal O. sinensis 0.106 0.001 distance Host insects & altitudinal O. sinensis 0.111 0.194 distance Partial Mantel tests considered for the correspondence between the first two matrices (the first column in the table) while setting constant the third (the second column). All tests were based on 100,000 permutations. Significant correlations are indicated by a P-value in bold.

Table 4 Results of cophylogenetic analyses with the event-based approach Jane 4 employing different cost schemes Event Event Total cost 1 costs C D D&S L F Default cost scheme 0, 1, 2, 1, 1 17 7 28 77 44 184 Cospeciation adjusted 1, 1, 2, 1, 1 14 12 26 78 44 200 -1, 1, 2, 1, 1 18 7 27 77 44 164 Duplication adjusted 0, 2, 2, 1, 1 19 3 30 79 44 189 0, 0, 2, 1, 1 17 14 21 88 44 174 Duplication & Host switch adjusted 0, 1, 3, 1, 1 18 14 20 93 44 211 0, 1, 1, 1, 1 15 1 36 73 44 154 Sorting adjusted 0, 1, 2, 2, 1 15 5 32 74 44 261 0, 1, 2, 0, 1 28 24 0 239 44 68 Failure to diverge adjusted 0, 1, 2, 1, 2 17 7 28 77 44 228 0, 1, 2, 1, 0 17 11 24 82 44 141 1 , Order of “event costs” (from left to right) is: C (cospeciation); D (duplication); D&S (duplication + host switch); L (loss/sorting); and F (failure to diverge)

This article is protected by copyright. All rights reserved.

Accepted Article

Table 5 The first 10 Pareto-optimal reconstructions from 10,000 random cost models with the event-based approach CoRe-PA Event Optimal cost Quality value Total cost C D S L C D S L 0.010 18.46 16 26 10 79 0.27 0.18 0.49 0.06 0.014 18.58 22 21 9 84 0.20 0.23 0.51 0.06 0.017 18.40 22 22 8 89 0.21 0.23 0.51 0.05 0.020 18.92 20 22 10 78 0.22 0.22 0.50 0.06 0.021 17.94 15 28 9 87 0.27 0.17 0.50 0.05 0.022 18.75 17 25 10 78 0.25 0.19 0.49 0.06 0.032 18.96 21 20 11 74 0.19 0.23 0.51 0.06 0.035 17.55 23 23 6 100 0.19 0.23 0.53 0.05 0.036 16.41 22 25 5 106 0.19 0.19 0.57 0.04 0.039 17.61 20 26 6 100 0.24 0.19 0.52 0.05 C, cospeciation; D, duplication; S, host switch; L, loss/sorting

This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.

Accepted Article This article is protected by copyright. All rights reserved.