|Home | About | Journals | Submit | Contact Us | Français|
Hirsutella rhossiliensis is a parasite of juvenile nematodes, effective against a diversity of plant-parasitic nematodes. Its global distribution on various nematode hosts and its genetic variation for several geographic regions have been reported, while the global population genetic structure and factors underlying patterns of genetic variation of H. rhossiliensis are unclear. In this study, 87 H. rhossiliensis strains from five nematode species (Globodera sp., Criconemella xenoplax, Rotylenchus robustus, Heterodera schachtii, and Heterodera glycines) in Europe, the United States, and China were investigated by multilocus sequence analyses. A total of 280 variable sites (frequency, 0.6%) at eight loci and six clustering in high accordance with geographic populations or host nematode-associated populations were identified. Although H. rhossiliensis is currently recognized as an asexual fungus, recombination events were frequently detected. In addition, significant genetic isolation by geography and nematode hosts was revealed. Overall, our analyses showed that recombination, geographic isolation, and nematode host adaptation have played significant roles in the evolutionary history of H. rhossiliensis.
IMPORTANCE H. rhossiliensis has great potential for use as a biocontrol agent to control nematodes in a sustainable manner as an endoparasitic fungus. Therefore, this study has important implications for the use of H. rhossiliensis as a biocontrol agent and provides interesting insights into the biology of this species.
It is well known that fungi possess diverse strategies to obtain nutrients, as saprobes, pathogens (of plant and animal hosts), and symbionts (with other microbes, plants, or animals) (1). One of the important living strategies is parasitoidism, where parasitoids sterilize and/or kill their hosts before the hosts reach reproductive age. Plant-parasitic nematodes are worldwide agricultural pests and are responsible for global agricultural losses amounting to an estimated $157 billion each year (2). Endoparasitic fungi are predominantly soil-dwelling organisms that have the ability to sterilize plant-parasitic nematodes (3). They produce characteristic spores that adhere to, penetrate, and consume nematode bodies (4, 5).
Hirsutella rhossiliensis (Ophiocordycipitaceae, Hypocreales, Ascomycota) (6) lives in the soil, and it has been mainly isolated from five nematodes including Criconemella xenoplax (7), Heterodera schachtii (8), Heterodera glycines (9), Rotylenchus robustus (10), and Globodera (11). Though a dominant parasitoid of H. glycines J2s (second-stage juveniles) in the United States, mainly from Minnesota, H. rhossiliensis occurs less frequently in northeast China (12). In the agricultural fields, the plants, nematodes, and endoparasitic fungi constitute a tight tripartite system with the plants providing food to nematodes and the nematodes serving as nitrogen sources to endoparasitic fungi. This existing natural interaction suggests that endoparasitic fungi could have great potential for use as biocontrol agents to control nematodes in a sustainable manner (13,–15).
To date, H. rhossiliensis has been detected in many localities in Europe, the United States, and China (11, 12). Genetic variation among different isolates of H. rhossiliensis has been detected using randomly amplified polymorphic DNA (RAPD) markers (16) and analyses of two housekeeping genes, internal transcribed spacer (ITS) of ribosomal DNA (rDNA) and the mitogen-activated protein kinase (MAPK) gene (13). Because of its global distribution and wide host range, H. rhossiliensis provides an ideal model to study the effects of geography and host nematodes on fungal population genetic divergence and factors driving genetic differentiation of fungi.
Investigating the evolutionary patterns and processes of nematode endoparasitic fungi over their geographical and ecological distribution provides valuable information on the radiation between nematode endoparasitic fungi and hosts (17). While phenotypic plasticity and genomic heterogeneity can lead to broad geographic and ecological distribution, low levels of gene flow often result in microbial speciation (18, 19). Population genetic analyses based on genotyping and genealogical analyses can help to infer past evolutionary events and evaluate the potential processes that generated the patterns of genetic variation (20, 21). An analysis of the population genetics of Hirsutella minnesotensis indicated that this dominant species in China had a clonal population structure and likely experienced a founder effect after colonization of China but without a significant bottleneck (22). However, the population genetic structure, possible geographic origins, and host nematode adaptation of H. rhossiliensis are largely unknown.
In recent years, single-nucleotide polymorphism (SNP) analyses have been widely used to detect intraspecific variation in microorganisms, including endoparasitic fungi (23, 24). When analyzing the genetic diversities of human fungal pathogens, for instance, Histoplasma capsulatum (25), Candida albicans (26), Aspergillus fumigatus (23), and Cryptococcus neoformans (27), it has been demonstrated that SNP analyses were highly discriminative in these species. Based on multilocus sequence analyses (MLSAs), genetic diversity in plant-pathogenic fungi was also reported for Magnaporthe grisea and Ustilaginoidea virens (28, 29).
Field experiments demonstrated that H. rhossiliensis was effective for control of H. glycines in greenhouses (30). Current biological control applications using H. rhossiliensis, however, have shown relatively limited success in agricultural fields. We think that an effective application of H. rhossiliensis or other fungi for biological control requires a thorough understanding of their population structure and their modes of reproduction. In this study, we used a multilocus sequence analysis (MLSA) scheme based on SNP markers to analyze population genetics of H. rhossiliensis by sampling from three geographic populations and five host nematode-associated populations. We aimed to understand (i) the level and partitioning of genetic variation within H. rhossiliensis, (ii) genetic structure and modes of reproduction in H. rhossiliensis, and (iii) whether geographic isolation and host nematode adaptation lead to rapid differentiation among populations of H. rhossiliensis. This study not only provides essential insights into the population differentiation and host adaptation of H. rhossiliensis but also helps scientists to better understand the ecological mechanism underlying the natural control of H. rhossiliensis in its role as a biocontrol agent against nematodes.
We employed all available H. rhossiliensis isolates (87 in total) around the world for this study. They were classified into three geographic populations (Europe, the United States, and China) and five host nematode-associated populations (Globodera, H. schachtii, C. xenoplax, R. robusta, and H. glycines). Except for 8 isolates from the CBS-KNAW Fungal Biodiversity Centre and 17 from the USDA-ARS Collection of Entomopathogenic Fungal Cultures (ARSEF), the other 62 isolates were all recovered by us from parasitized H. glycines in Heilongjiang Province, China, or northern Minnesota, USA (Table 1 and Fig. 1; see also Table S1 in the supplemental material).
To isolate H. rhossiliensis from parasitized H. glycines, we collected rhizosphere soil samples from soybean fields in the main soybean-producing regions across northeast China from 2007 to 2008 and across Minnesota (USA) from 2010 to 2011. The sucrose flotation and centrifugation method (31) was used to extract the second-stage juveniles (J2s) that naturally occur in the soil, and these J2s were then transferred to wells of a 6-well tissue culture plate. The extracted J2s were examined with an inverted microscope (Olympus, Japan) at ×40 magnification, washed twice with sterilized water, and then placed on Difco potato dextrose agar (PDA) plates at 25°C for 2 weeks. All of the filamentous fungi isolated from these J2s were cultured on PDA plates covered with a piece of cellophane paper, and 0.05 to 0.1 g of mycelia was harvested and used to extract genomic DNA with the cetyltrimethylammonium bromide (CTAB) method (32). Fungal cultures were identified based on their morphology (30) and nonribosomal DNA (nrDNA) ITS sequences (13).
We initially chose three housekeeping genes as molecular markers: the nrDNA ITS region, the beta-tubulin gene (tub), and the translation elongation factor 1α gene (tef). However, the number of phylogenetically informative sites provided by the three fragments was insufficient. Therefore, novel DNA markers were determined by searching the genome data (15). Specifically, 27 fragments of genes putatively involved in the infection processes of H. rhossiliensis into nematodes and 15 fragments of interest as a result of a genomic comparison between H. minnesotensis and H. rhossiliensis were chosen as candidate markers. PCR primers were designed for each of the 42 fragments using the online software Primer3 (http://frodo.wi.mit.edu/primer3/). From 12 randomly chosen isolates of different origins (CBS104.94, CBS105.94, ARSEF2006, ARSEF3746, ARSEF2788, ARSEF2789, ARSEF3755, ARSEF3756, CN-13-1, CN-30-1, USA-87-2, and USA-92-1) (see Table S1 in the supplemental material), we amplified and sequenced each of the 42 DNA fragments. Fragments with a high frequency of nucleotide variation were finally chosen as markers used in this study. A neutrality test for the selected markers was performed by DnaSP 5.10 (33).
Nine molecular markers, including ITS, tub, tef, and six newly identified markers, were amplified from each of the 87 H. rhossiliensis isolates. Each PCR mixture (50 μl) was composed of 5 μl of 10× EasyTaq buffer, 5 μl of deoxynucleoside triphosphate (dNTP) mixture (2.5 mM), 2 μl of each primer (10 μM), 1 μl of genomic DNA, 0.5 μl of EasyTaq DNA polymerase (TransGen Biotech, Beijing, China), and 36.5 μl of double-distilled water (ddH2O). The PCRs were performed in a TGradient thermocycler (Biometra, Germany) under the following conditions: denaturation at 94°C for 3 min; 35 cycles of denaturation at 94°C for 40 s, annealing at variable temperatures for 50 s, and elongation at 72°C for 1 min; and a final 10-min elongation step at 72°C. After confirmation of the PCR products by agarose gel electrophoresis, the amplicons were cleaned using the 3S Spin PCR product purification kit (Biocolor Bioscience & Technology Company, China). The purified PCR products were then sequenced using an ABI 3730 XL DNA sequencer (Applied Biosystems, USA) with Sanger dideoxy terminator sequencing. In the case of poor sequencing quality resulting from the direct sequencing of PCR products, the representative PCR fragments were cloned, and individually cloned fragments were then sequenced again.
Sequences of the fragments were aligned with the Muscle algorithm of MEGA 6.0 (34) and trimmed to remove ambiguously aligned regions. The nucleotide diversity (π), haplotype diversity, number of polymorphic sites, and allele numbers were estimated for both individual fragment and combined sequences using the program DnaSP 5.10 (33). By employing DnaSP 5.10, alleles at each locus were assigned and combined into an allelic profile designated a haplotype (HA) for each strain.
We performed a partition homogeneity test (PHT) to detect potential phylogenetic conflicts between different loci using PAUP version 4.0b10 (Sinauer Associates, Sunderland, MA, USA; D. Swofford, 2003). PHT used informative characters and simple stepwise-addition heuristic searches with 1,000 replicates. A P value lower than 0.001 indicated statistically significant differences (35). After identification of loci with no significant conflict, the phylogenic relationship among 87 strains was constructed using two methods: Bayesian inference (BI) and maximum likelihood (ML). To conduct phylogenetic analyses, we used the program PartitionFinder 1.1.1 (36) to evaluate the best partitioning scheme and determine suitable substitution models under the Bayesian information criterion (see Table S5 in the supplemental material). For partitioned data sets, model parameters across different partitions were unlinked, and the overall rate of evolution among partitions was allowed to vary. The Bayesian analysis, which was performed in MrBayes 3.2.2 (37), was initiated with random starting trees, run for 100 million generations with four incrementally heated chains (Metropolis-coupled Markov chain Monte Carlo ), and sampled at intervals of 10,000 generations. To avoid entrapment in the local optima, two independent Bayesian analyses were run, and the log-likelihood scores were compared for convergence (38). An ML analysis was performed in RAxML version 8.0.0 (39) and conducted under the GTRGAMMA model, and the ML support was assessed by 1,000 bootstrap replicates. Two H. minnesotensis strains, CN3608 and CBS113353, were used as the outgroup.
A phylogenetic network was constructed using the median-joining method implemented in the program SplitsTree 4.12 (40) to further assess the relationships among worldwide H. rhossiliensis haplotypes. This model-free method uses a parsimony approach based on pairwise differences to connect each sequence to its closest neighbor and allows for the creation of internal nodes (median vectors), which are interpreted as unsampled or extinct ancestral genotypes to link the existing genotypes in the most parsimonious way (41).
To determine whether the number of loci was sufficient to represent the genotypic diversity of the populations, we plotted the relationship between the number of loci and the genotype diversity using MultiLocus1.3b (42). A population genetic analysis was also performed within and between the three geographic populations and within and between the five host nematode-associated populations. For the within-population genetic analysis, we analyzed the genetic variation patterns with DnaSP 5.10, including the number of haplotypes, nucleotide diversity, and haplotype diversity. For the cross-population genetic analysis, pairwise differentiation coefficients (Fst) (43) were calculated and tested for significance against 1,000 bootstrap replicates using Arlequin 3.1 (44), and the average gene flow (Nm) (45) was calculated using DnaSP 5.10.
Sexual reproduction generates random associations between alleles at different loci because genes from different individuals are mixed in a given population by sexual reproduction. Therefore, when recombination occurs among populations, random associations between alleles may be observed at different loci, even if sexual reproduction has not been observed in the natural populations (46). To identify evidence for recombination in H. rhossiliensis, two methods were used: the proportion of compatible pairs of loci (PrCP) (47) and the Φw test (or pairwise homoplasy index [PHI]) (48). Tests for recombination based on the principle of compatibility have been proven to be the most powerful (49, 50). Two loci are compatible (PrCP = 1) if it is possible to account for all the observed genotypes by mutations without inferring homoplasy (reversals, parallelisms, or convergences) or recombination; otherwise, the loci are incompatible (PrCP < 1). The Φw test means that phylogenetic incompatibility is an indicator of recombination at the population level; the lack of phylogenetic incompatibility, in contrast, implies asexual reproduction. The tests for PrCP were calculated using MultiLocus1.3b (42) with 1,000 randomizations; the Φw test was calculated by SplitsTree 4.12.
An analysis of molecular variance (AMOVA) was performed to determine the relative contribution of genetic variation among and within populations using GenAlEx 6.5 (51). The AMOVA was conducted for both the three geographic populations and the five host nematode-associated populations.
Mantel tests were performed in GenAlEx 6.5 for three groups of correlations: between the fungal genetic distance and longitudinal-latitudinal coordinates of geographic distance, between the fungal genetic distance and longitudinal distance, and between the fungal genetic distance and latitudinal distance. For the genetic distances, we used the pairwise population PhiPT values, which represent genetic distances among all pairs of populations as a tri-matrix, and these distances were then compared with the geographic distances between populations.
All of the nucleotide sequences obtained in this study have been submitted to the GenBank database, and the accession numbers are KP885720 to KP886014 (see also Table S2 in the supplemental material).
To determine appropriate markers, we successfully amplified the 42 candidate fragments from each of the 12 test isolates. Six fragments (HIR_00714, HIR_06686, HIR_08121, HIR_08679, HIR_01569, and HIR_09423) showed a high frequency of nucleotide variation and existed as single copies as indicated by BLAST searches against the H. rhossiliensis genome. The six fragments, together with three commonly used markers (nrDNA ITS, tub, and tef), were finally chosen to amplify from each of the 87 isolates (Table 2). Direct sequencing of PCR products from any of the nine fragments indicated no heterogeneity in sequencing chromatograms. Three to 11 alleles were identified in our samples depending on the locus. Moreover, the result of a neutrality test showed that all nine markers were neutral (Table 2).
A significant conflict was detected by PHT between HIR_00714 and the other loci (see Table S4 in the supplemental material). Therefore, HIR_00714 was excluded from the concatenated data set in the following analyses. From the concatenation of the remaining eight loci, a total of 280 variable sites within the combined 4,658-bp sequence and 44 multilocus haplotypes were identified from the 87 isolates (Table 2).
Bayesian (BI) and maximum likelihood (ML) phylogenetic trees were constructed based on the 8-locus concatenated data set. H. rhossiliensis strains tend to cluster by geography and nematode hosts, and each of the two phylogenetic approaches identified six clusters: European strains from Globodera sp. (NL-GL), American strains from C. xenoplax (USA-CX), American strains from R. robusta (USA-RR), American strains from H. schachtii (USA-HS), American strains from H. glycines (USA-HG), and Chinese (northeastern) strains from H. glycines (CN-HG). However, there were subtle topological differences between the two approaches. For example, two American strains from H. schachtii (CBS567.92 and ARSEF3761) did not group into the USA-HS subclade, and the positions of two Chinese strains from H. glycines (CN-30-2-JX and CN-3-10-JL) and one American strain from C. xenoplax (ARSEF3746) had inconsistent placements between the BI and ML trees (Fig. 2; see also Fig. S2 in the supplemental material).
The network diagram showed a linear structure among the 87 strains (Fig. 3), with Chinese strains from H. glycines, American strains from H. glycines, European strains from Globodera sp., American strains from H. schachtii, American strains from R. robusta, and American strains from C. xenoplax each forming an individual cluster.
A high value for the genetic differentiation coefficient Fst (P < 0.01) revealed that significant genetic differentiation occurred among all of the geographic or host nematode-associated populations except between NL-P and USA-P, the value for which was not significant (P > 0.01) (Tables 3 and and4).4). Consistently, the estimated gene flow was low among the three geographic populations (Nm = 0.14) and among the five host nematode-associated populations (Nm = 0.28).
The PrCP index provided evidence for recombination except for two populations, American strains from C. xenoplax (USA-CX) and American strains from R. robusta (USA-RR). When the 87 strains were analyzed together, recombination signals were found in 3/5 locus pairs (Table 3; see also Table S6 in the supplemental material).
Similarly, the Φw test found statistically significant evidence for recombination (P < 0.001). However, there was a slight difference in the Φw results for recombination among the three geographic populations and five host nematode-associated populations compared with the PrCP results. For example, the Φw test did not find statistically significant evidence for recombination in the USA-P, USA-CX, USA-HS, or USA-RR population (Table 3).
The AMOVA results showed that 44% of the genetic variation could be attributed to variation among the three geographic populations and 56% could be attributed to variation among the individual isolates within populations (Table 5). For the five host nematode-associated populations, the AMOVA results showed that 61% and 39% of the genetic variation could be attributed to variation among the populations and among individual isolates within populations, respectively. Both of these sources of variation were significant (P < 0.01) (Table 6).
The Mantel tests revealed significant genetic isolation by distance irrespective of the horizontal geographical distances (based on longitudinal-latitudinal coordinates, along a longitudinal gradient, or along a latitudinal gradient) or population divisions (three geographic populations) employed (Table 7).
We collected 87 H. rhossiliensis strains from five host nematodes in Europe, the United States, and China. Although a large number of strains were obtained from H. glycines, at least three strains were included from each of the other nematode hosts (Table 1; see also Table S1 in the supplemental material). We would explain that our samples from four nematodes (Globodera sp., H. schachtii, C. xenoplax, and R. robusta) were isolated in the 1980s and 1990s, while samples from H. glycines were isolated during 2007 to 2011. Because some nematodes did not occur in China or we did not find this fungus from them (11, 52,–54), we failed to isolate H. rhossiliensis from other host nematodes except from H. glycines. It was also difficult to obtain strains from the above four host nematodes (except H. glycines) in the United States. The cultures from culture collection centers (CBS and ARSEF) were preserved in lyophilized forms, which largely decreased strain mutations over the past 20 to 30 years (see Table S1 in the supplemental material). Although there might be some genetic variations of fungus in nature, we believe that the variation during 20 to 30 years could not affect the main results obtained from this study.
In this study, recombination was found in most of the H. rhossiliensis populations and when the whole samples were considered. It can be inferred and supported by two lines of evidence, the PrCP values and the Φw test (Tables 3 and and4).4). The PrCP test measures whether two loci are compatible, and the Φw test examines phylogenetic incompatibility. So, it is reasonable that there is conflict between the two results (49). Indeed, the sexual cycles of most fungi are difficult to observe in nature (55), and inferences of the potential for a sexual cycle to occur have largely relied on analyses of gene and genotype frequencies in natural populations. The sexual stage of H. rhossiliensis has not been discovered in nature or in the laboratory, although we detected mating type genes in its genome (unpublished data). In plant fungal pathogens, a high level of recombination provides an advantage in rapidly generating many new combinations of virulence genes to counterbalance corresponding resistance genes in the host (56, 57). In arbuscular mycorrhizal fungi, recombination or recombination-like events in addition to clonality have greatly contributed to genetic diversity (58). This study has demonstrated evidence of recombination among H. rhossiliensis strains, and recombination appears to be one of the factors shaping the evolution of H. rhossiliensis. Similar results were also observed in Candida albicans, which showed both clonality and recombination, even though a complete sexuality stage is not known to exist in this fungus (59). However, sexual recombination may not be the sole factor because we cannot exclude the possibility that parasexual recombination has contributed to the observed linkage equilibrium and phylogenetic incompatibility (60).
Our population genetics analyses of H. rhossiliensis were based on two types of population divisions: three geographic populations (NL-P, USA-P, and CN-P) and five host nematode-associated populations (NL-GL, USA-CX, USA-RR, USA-HS, and USA/CN-HG) (Table 1). The significance of geography and the host nematodes in shaping the structure of H. rhossiliensis populations on a large scale was supported by the strong divergent signals in the phylogenetic and haplotype network analyses (Fig. 2 and and3).3). In this study, a significant positive correlation was observed between H. rhossiliensis genetic distance and latitudinal and longitudinal geographic distance, thus indicating that geography is an important barrier to gene flow and has promoted the genetic divergence within fungal species (Table 7). Moreover, differentiation among the three geographic populations was also observed from high Fst values and low Nm values. We suggest that the North Pacific Ocean may be a boundary separating the Chinese population and the American population of H. rhossiliensis because of limited opportunities to disperse via wind or soil. A similar geography-based separation has been reported in Colletotrichum truncatum from chili peppers in China that were genetically differentiated into southern and northern populations (61). Previously, we found significant differences among local populations of H. minnesotensis from China and the United States, where isolation by geography plays a key role in genetic differentiation (22). Based on these results, we conclude that geographic distance contributes to driving population genetic differentiation of H. rhossiliensis.
Among the five host nematode-associated populations studied here, we identified positive correlations between H. rhossiliensis genetic distance at different latitudinal and longitudinal distances. We observed 61% of the molecular genetic variance among those nematode-associated populations. This might result from host-related selection driving population genetic differentiation at neutral and selected loci (62). Another explanation was that host shifts might occur frequently between neighboring nematodes, which increased the mating probability of H. rhossiliensis across the nearby hosts and also increased the genetic variation within and among populations.
Our analyses revealed that both nematode hosts and geography contributed to genetic differentiation among H. rhossiliensis populations (Tables 5 to to7).7). We, however, could not perform a two-way AMOVA to further look at their relative effects because population assignments of our samples did not meet the requirement for performing such an analysis. In this study, we confirmed that difference in host nematode was a key factor in regulating genetic divergence of nematode parasitic fungus from the phylogenetic approaches and the haplotype network (Fig. 2 and and3).3). Similarly, the strong association of genetic groups with different geographic locations suggested the important role of geographic isolation. Previous studies on microorganisms such as Bacillus simplex concluded that genetic differentiation was strongly attributed to geographic locations (63). Numerous studies found that geographic region and host species contribute equally to the population genetic differentiation, for instance, Grosmannia clavigera (64), Claviceps purpurea (65), Venturia inaequalis (66), Arthrobotrys oligospora (67), Colletotrichum gloeosporioides (68), and Metarhizium anisopliae (69).
In conclusion, our study identified the population genetic structure of the nematode endoparasitic fungus H. rhossiliensis and found a high level of genetic variation among its populations worldwide. Using multiple sources of evidence, we demonstrated that recombination among strains, isolation by geography, and isolation based on the host nematode were factors driving the genetic differentiation of H. rhossiliensis. Therefore, the findings of this study on the population structure of H. rhossiliensis and the associated multiple driving factors affecting its genetic diversity provide an essential understanding of the mechanisms of evolution and differentiation of the fungus on a global scale. These findings are expected to provide important insights into nematode control by this fungus.
This study was supported by the National Natural Science Foundation of China (grant no. 31430071) and the National Basic Research Program of China (973 program, grant no. 2013CB127506).
We thank Richard A. Humber, who kindly provided 15 strains from the USDA-ARS Collection of Entomopathogenic Fungal Cultures (ARSEF) database at the USDA-ARS Biological Integrated Pest Management Research Unit, Cornell University, New York, USA.
Supplemental material for this article may be found at http://dx.doi.org/10.1128/AEM.01708-16.