Out in the Cold: Identification of Genomic Regions Associated With Cold Tolerance in the Biocontrol Fungus Clonostachys rosea Through Genome-Wide Association Mapping.

There is an increasing importance for using biocontrol agents in combating plant diseases sustainably and in the long term. As large scale genomic sequencing becomes economically viable, the impact of single nucleotide polymorphisms (SNPs) on biocontrol-associated phenotypes can be easily studied across entire genomes of fungal populations. Here, we improved a previously reported genome assembly of the biocontrol fungus Clonostachys rosea strain IK726 using the PacBio sequencing platform, which resulted in a total genome size of 70.7 Mbp and 21,246 predicted genes. We further performed whole-genome re-sequencing of 52 additional C. rosea strains isolated globally using Illumina sequencing technology, in order to perform genome-wide association studies in conditions relevant for biocontrol activity. One such condition is the ability to grow at lower temperatures commonly encountered in cryic or frigid soils in temperate regions, as these will be prevalent for protecting growing crops in temperate climates. Growth rates at 10°C on potato dextrose agar of the 53 sequenced strains of C. rosea were measured and ranged between 0.066 and 0.413 mm/day. Performing a genome wide association study, a total of 1,478 SNP markers were significantly associated with the trait and located in 227 scaffolds, within or close to (< 1000 bp distance) 265 different genes. The predicted gene products included several chaperone proteins, membrane transporters, lipases, and proteins involved in chitin metabolism with possible roles in cold tolerance. The data reported in this study provides a foundation for future investigations into the genetic basis for cold tolerance in fungi, with important implications for biocontrol.

Introduction

As the human global population grows and plant pathogens and pests increase and spread, sustainable agricultural production, and disease control increases in importance to ensure future food security (Syed Ab Rahman et al., 2018). Biological control agents (BCAs) demonstrate an environmentally sustainable and potent alternative to current industrial pesticide practices in agriculture, due to the risk chemicals may pose to human and animal health (Syed Ab Rahman et al., 2018). The mycoparasitic fungus Clonostachys rosea (Link: Fr.) Schroers, Samuels, Seifert & W. Gams, comb. nov. (Schroers et al., 1999) is a widely studied BCA, effective against several economically important fungal plant pathogens such as Botrytis cinerea, Bipolaris sorokiniana, Fusarium spp., and Sclerotinia sclerotiorum (Knudsen et al., 1995; Sutton et al., 1997; Teperi et al., 1998; Xue et al., 2009; Rodríguez et al., 2011). The sexual form (teleomorph) is described as Bionectria ochroleuca (Schw.) Schroers & Samuels, and certain strains readily produce sexual fruiting bodies (perithecia) in single culture indicating its homothallic mode of reproduction (Schroers et al., 1999). The biocontrol effect of C. rosea can be derived from its ability to directly parasitize plant pathogenic fungi, but also from exploitation competition for space and nutrients and from interference competition via antibiosis (Li et al., 2002; Rodríguez et al., 2011; Jensen et al., 2017).

Active growth in colder temperatures (below 15°C), is an important trait for a BCA to be effective in the field in temperate climates such as in the Nordic countries. This enables rapid colonization of seeds, wounds and other substrates for efficient exploitation competition, but is also important for maintaining high metabolic activity. This trait was originally one of the reasons for the selection of certain C. rosea strains as BCAs, e.g., for control of F. culmorum on wheat and barley at 15°C (Knudsen et al., 1995; Teperi et al., 1998). On the molecular level, mechanisms of cold tolerance in fungi partly overlap with other abiotic stress tolerances. Estimations suggest that approximately 85% of the biosphere is constantly subjected to temperatures of 5°C or lower (Li et al., 2012). Cold tolerance can be associated with changes in cell wall composition (carbohydrates, proteins and pigments), changes in lipid and osmolyte composition of cell membranes and the RNA translation and protein folding machinery (RNA helicases, trehalose, chaperone proteins), alcohol metabolism, transport proteins and mitochondria (Ellison et al., 2011; Blanc et al., 2012; Treseder and Lennon, 2015). Temperature such as cold has also been identified as a major environmental factor for organisms to adapt to, as it is one of the major environmental forces limiting life (D'Amico et al., 2002; Blanc et al., 2012; Li et al., 2012). Cold adaptation has been studied across free-living and plant-associated fungi to assess the temperature effects on adaptation, for example as a driving force for environmental colonization in plant-associated fungi and population divergence in fungi such as Neurospora crassa (Robinson, 2001; Ellison et al., 2011; Li et al., 2012). For fungi, cold tolerance may further be associated with osmotic stress tolerance, and spore survivability in low water, low nutrient conditions (Ruisi et al., 2007).

As genome sequencing has become cheaper and more user friendly, performing large-scale comparative genomics is now a standard practice for the identification of single nucleotide polymorphisms (SNPs), and determining their association with phenotypes of interest using genome-wide association studies (GWAS). The GWAS approach allows for statistical determination of SNPs of interests over a population of individuals with a shared evolutionary history. Previously, GWASs have been performed on fungal pathogens such as Heterobasidion annosum, Zymoseptoria tritici, Parastagonospora nodorum, and Fusarium graminearum (Dalman et al., 2013; Gao et al., 2016; Talas et al., 2016; Hartmann et al., 2017). However, as far as we are aware, no GWAS on fungal BCAs have been performed to date. The C. rosea strain IK726 genome was originally sequenced using short reads of the Illumina platform (Karlsson et al., 2015). However, the modern sequencing technology of the PacBio platform is allowing for more comprehensive and accurate genomics, by sequencing reads that are thousands to hundreds of thousands of base pairs (bp) in length instead of a few hundred bp (Roberts et al., 2013). This provides a very potent framework when combined with short read sequencing for assembling genomes (Schadt et al., 2010; Rhoads and Au, 2015; Chakraborty et al., 2016).

In this study, we evaluate the potential for performing GWAS in a fungal BCA. We present a re-sequenced genome of C. rosea IK726, combining the previously published genome scaffolds with long contigs gained from PacBio sequencing into a second, improved version of the genome, expanding the size from 58.3 to 70.7 Mbp and adding 6,497 predicted genes to a total of 21,246. In addition, we used Illumina technology to re-sequence the genomes of 52 additional C. rosea strains and performed a GWAS for growth rate in cold temperature (10°C). The analysis identified 1,478 SNP markers associated with the trait and certain genes located in the vicinity were predicted to encode several chaperone proteins, membrane transporters, lipases, and proteins involved in chitin metabolism with possible roles in cold tolerance. This genomic study updated the C. rosea IK726 genome, established a population genomic dataset for C. rosea and revealed important genomic regions for future studies of cold tolerance in this very promising BCA.

Materials and methods

Fungal growth and culture conditions

Clonostachys rosea strains (Supplementary Table 1) were revived from glycerol stocks stored at −80°C, and maintained on potato dextrose agar (PDA; Oxoid, Cambridge, UK) at 25°C in darkness. Two growth rate assays were performed in triplicates for each strain by inoculating ½ strength PDA petri dishes with an agar plug with actively growing mycelium, followed by incubation at 10°C and 25°C, respectively, in darkness. Growth measurements of strains at 10°C were made by marking the edges of the fungal mycelium in four places in a square shape on the petri dish at 4 days and 28 days post inoculation. The mean change in millimeters between the 4-day and 28-day marks was measured and divided by 24 to provide data as mm growth per day. Growth measurements at 25°C were done in the same manner, but measured from the edge of the inoculation agar plug to the mycelial front at 5 days post inoculation. Growth rate data at 10°C were analyzed by analysis of variance (ANOVA) using the Welch's test implemented in Minitab ver. 18 (Minitab Inc., State Collage, PA), while pairwise comparisons were made using the Games-Howell method at the 95% significance level. Growth rate data at 25°C were analyzed by ANOVA using a general linear model approach implemented in Minitab ver. 18 (Minitab Inc., State Collage, PA). Pairwise comparisons were made using the Fisher method at the 95% significance level. Linear regression analysis was performed by Pearson correlation in Minitab ver. 18 (Minitab Inc., State Collage, PA). Principle component analysis, using the R programming language, was performed in order to detect whether there was any correlation between geographical location of isolate origin and growth at 10°C.

For DNA extraction, strains were grown in 200 ml liquid Czapek-Dox medium (Sigma-Aldrich, Steinheim, Germany), Vogel's minimal medium (Vogel 1956) or malt extract (1.75 %) with peptone (0.25 %) medium at room temperature, shaking at 120 rpm. Cultures were harvested after 3–13 days, depending on growth rate, by snap freezing in liquid nitrogen, and then freeze-dried.

Genomic DNA extraction and sequencing

High-quality genomic DNA was extracted using CTAB/chloroform extraction for Illumina sequencing or Qiagen-tip 100 (Qiagen, Hilden, Germany) for PacBio sequencing (Menkis et al., 2015). Starting material for both protocols were freeze-dried mycelium ground with sand in a mortar and pestle. The CTAB/chloroform was scaled up to 50 ml falcon tubes with 10 ml CTAB-buffer and two to four rounds of chloroform extraction. Qiagen-tip 100 DNA extractions followed the Qiagen genomic DNA handbook with modifications of the yeast protocol. Instead of making spheroplasts, extraction started with adding G2 buffer to the ground mycelium and thereafter following the protocol. Regardless of protocol, at the precipitation step DNA was harvested with a glass hook and the remaining DNA was pelleted by centrifugation. For PacBio sequencing the harvested DNA was used, while for Illumina sequencing the harvested DNA was mostly used but the precipitated DNA was used in a few cases. DNA quality and concentrations were measured with Nanodrop (Thermo Fisher Scientific, Waltham, MA), Qubit dsDNA BR assay kit (Thermo Fisher Scientific, Waltham, MA) and agarose gel electrophoresis.

Base coverage of the C. rosea IK726 genome was generated using PacBio RSII Technology with an insert length of 20 kbp using standard library preparation kits. Base coverage of additional C. rosea genomes were generated using Illumina HiSeqX paired end sequencing with an insert length of 350 bp and read length of 150 bp using standard library preparation kits. The PacBio long reads were assembled into polished contigs using HGAP and the SMRT workflow (Chin et al., 2013).

Genome assembly and annotation, SNP calling and population structure analysis

The previously reported genome assembly of C. rosea IK726 (Karlsson et al., 2015) was combined with the polished contig output of the PacBio sequencing using quickmerge (Chakraborty et al., 2016). The PacBio contigs and the original assembly were compared using MUMmer ver. 3.23 (Delcher et al., 2003). The merged genome (version 2) was annotated using a pipeline described previously (Karlsson et al., 2015). For the additional 52 C. rosea strains, Illumina reads were aligned against the C. rosea IK726 ver. 2 genome using Bowtie 2 ver. 2.2.4 (Langmead and Salzberg, 2012). The alignments were further analyzed with SAMtools ver. 0.1.18 for filtering out PCR duplicates (Li et al., 2009). Variant calling was performed on filtered alignments using a combination of Freebayes ver. 1.0 and the bamaddrg script for adding read groups to binary alignment map (BAM) files (Garrison and Marth, 2012; Garrison, 2018). Population structure analysis was performed on the 53 isolates using Structure ver. 2.3.4 (Hubisz et al., 2009) using 8,000 randomly sampled SNPs across the population. Linkage disequilibrium (LD) was determined by using vcftools ver. 0.1.15, and to calculate r2 using a 10 kbp window (Danecek et al., 2011).

Genome-wide association study

The GWAS was performed using PLINK ver. 1.90 (Purcell et al., 2007), using the parameters—maf 0.1—hwe 1e-5 for SNP filtering. The SNPs were annotated using the ANNOVAR software (Wang et al., 2010). The resulting beta and beta standard deviation from the PLINK analysis were used as input to the R package ashr, for empirical Bayesian multiple hypothesis testing, for estimating local false sign rate (lfsr) for more appropriate significance estimation (Stephens, 2016). We used an lfsr of 0.05 as cut-off for significance. In addition, we used false discovery rate (FDR) correction for multiple testing as implemented in PLINK. Identification of conserved protein modules and features were made using the conserved domain database (CDD, Marchler-Bauer et al., 2017) and the transporter classification database (TCDB, Saier et al., 2016).

Chitin synthase phylogenetic analysis

The amino acid sequences of the C. rosea chitin synthetases CRV2G00008512 and CRV2G00019295 were used to search GenBank at the National Center for Biotechnology Information (NCBI) for similar sequences using the BlastP algorithm. Amino acid sequence alignment was performed using MUSCLE (Edgar, 2004; Tamura et al., 2013), in the MEGA6 software. Phylogenetic analysis was performed using Maximum Likelihood methods implemented in MEGA6. Pairwise gap deletion was performed and the substitution model used was WAG with gamma distributed rates. Statistical support was performed using 1,000 bootstrapped resamplings of the tree.

Nucleotide accession numbers

Sequencing analysis files for all isolates in this study except C. rosea IK726 (deposited separately described by Karlsson et al., 2015) were deposited under the study accession number PRJEB26874 at the European Nucleotide Archive/European Molecular Biology Laboratory–European Bioinformatics Institute repository. PacBio sequencing files for the re-sequencing of C. rosea strain IK726, and the updated genome annotation were also deposited under ENA project PRJEB26874.

Results

Growth rates at 10°C and 25°C differs between C. rosea strains

Growth rates of 53 C. rosea strains at 10°C were normally distributed, ranged from 0.066 to 0.413 mm/day (Table 1), and there were significant differences between strains (P < 0.001). The lowest growth rates were recorded for strains CBS 708.97, SHW-1-1, 1881, and IK726 originating from USA, China, Slovenia and Denmark, respectively. The highest growth rates were recorded for two strains from Slovenia (1832 and 1882) and strain CBS 188.33 from the Netherlands. No connection between growth rate and geographic origin was detected (Supplementary Figure 1). Growth rates at 25°C were normally distributed, ranged from 0.940 to 3.796 mm/day (Table 1), and there were significant differences between strains (P < 0.001). There was a weak correlation between growth rates at 10°C and 25°C (r2 = 0.178, P = 0.025), indicating that the mechanisms for cold tolerance were distinct from the basis for high growth rate at ambient temperature. Therefore, the growth rate at 10°C for each strain was normalized by the growth rate at 25°C (Table 1) for the GWAS analysis, to enable identification of SNPs and genomic regions specifically contributing to higher growth rate at colder temperature.

Table 1

Growth rates of Clonostachys rosea strains on PDA at 10°C and 25°C, and adjusted growth at 10°C.

Isolate	Cold growth mm/day	PDA growth mm/day	Cold rate/PDA rate	Origin
GG-1-2	0.201	0.940	0.214	Asia
2177	0.326	1.556	0.210	Europe
CBS 705.97	0.399	2.202	0.181	North America
JLB-7-1	0.347	1.923	0.181	Asia
IK216	0.326	1.811	0.180	Europe
2178	0.358	2.020	0.177	Europe
1832	0.410	2.416	0.170	Europe
2176	0.295	1.755	0.168	Europe
1316	0.326	2.022	0.161	Europe
IK371	0.309	2.035	0.152	Europe
B13	0.396	2.683	0.148	Oceania
1701	0.260	1.881	0.138	Europe
CBS 188.33	0.410	3.000	0.137	Europe
CBS 376.55	0.399	2.952	0.135	North America
CBS 193.94	0.372	2.850	0.130	South America
GL-1-1	0.313	2.474	0.126	Asia
1882	0.413	3.290	0.126	Europe
1829	0.194	1.550	0.125	Europe
1833	0.233	1.857	0.125	Europe
CBS 178.28	0.333	2.800	0.119	Europe
CBS 287.78	0.226	1.921	0.117	North America
CMI192798	0.344	2.930	0.117	Europe
GJS89-34	0.295	2.525	0.117	South America
CBS 154.27	0.326	2.828	0.115	North America
CBS 907.72E	0.351	3.061	0.115	Asia
CBS 907.72D	0.281	2.458	0.114	Asia
1884	0.208	1.856	0.112	Europe
CBS 100502	0.233	2.141	0.109	Europe
CBS 148.72	0.278	2.588	0.107	Europe
1827	0.174	1.627	0.107	Europe
SDT-5-1	0.208	1.985	0.105	Asia
CBS 421.87	0.285	2.825	0.101	Europe
CBS 222.93	0.219	2.210	0.099	South America
1830	0.215	2.175	0.099	Europe
SYP-4-2	0.191	1.971	0.097	Asia
1883	0.278	2.892	0.096	Europe
CBS 704.97	0.160	1.724	0.093	North America
NHH-61-2	0.188	2.074	0.090	Asia
STG-21-1	0.142	1.586	0.090	Asia
CBS 569.69	0.191	2.200	0.087	Europe
JXLS-1-1	0.267	3.240	0.083	Asia
CBS 649.8	0.194	2.383	0.082	Africa
SHW-1-1	0.122	1.561	0.078	Asia
IBT7519	0.194	2.531	0.077	Europe
CBS 100000	0.201	2.776	0.073	Oceania
NHH-48-2	0.174	2.426	0.072	Asia
CBS 907.72G	0.222	3.315	0.067	Asia
1885	0.174	2.700	0.064	Europe
CBS 706.97	0.184	2.969	0.062	North America
CBS 277.5	0.188	3.045	0.062	North America
1881	0.128	2.446	0.053	Europe
IK726	0.132	3.138	0.042	Europe
CBS 708.97	0.066	2.002	0.033	North America

PacBio sequencing results in an improved genome assembly of C. rosea IK726

PacBio sequencing of the C. rosea IK726 genome resulted in 71.2 Mbp sequence data. The PacBio reads were subjected to error correction of the longest reads by subread filtering, mapping and de novo assembly of all reads into long polished contigs using HGAP ver 3.0 and the SMRT pipeline (Chin et al., 2013; Roberts et al., 2013). Combining the polished contigs of the PacBio sequencing with the previously reported Illumina-based genome of C. rosea IK726 (Karlsson et al., 2015), expanded the genome from 58.3 (ver. 1) to 70.7 (ver. 2) Mbp across 767 scaffolds (Supplementary Table 2). N50 increased from 790 kbp to 1.9 Mbp, and the gap length decreased from 924 to 207 kbp (Supplementary Table 2). Comparisons between ver. 1 and ver. 2 using MUMmer revealed a strong correlation between the longer PacBio contigs and the ver. 1 genome (Supplementary Figure 2). Several smaller PacBio contigs were left without a corresponding scaffold of the old version, demonstrating the superior coverage generated by PacBio. A total of 21,246 genes were predicted, although this number dropped to 17,508 when excluding exclusively ab initio predicted genes.

GWAS identifies SNPs associated with growth rate at 10°C

The 52 Illumina-sequenced genomes exhibited an average alignment percentage to the C. rosea genome ver. 2 of 91% (Supplementary Table 1). The genomes had an average of 2.6 Gbp of total reads with a Q value of ≥30, resulting in an average coverage of 35X. Initially, we used vcftools to extract a total of 1,898,673 SNPs. These were further filtered during the GWAS analysis with PLINK, filtering SNPs using a minor allele frequency cut-off of 10% and a Hardy-Weinberg statistic of 10−5, leaving 63,726 SNPs. Based on 8,000 randomly sampled SNPs from the vcftool results, no population structure among the 53 strains was detected using the Structure software. The LD half decay distance was determined to 625 bp (Figure 1). Overall, the unfiltered SNPs were identified across 717 scaffolds (Figure 2, Supplementary Table 3).

Figure 1

The linkage disequilibrium half decay is 625 base pairs in this population of Clonostachys rosea. Linkage disequilibrium of all identified SNPs in a 10,000 base pair window. The green line demonstrates the linkage disequilibrium half decay mark, the blue line highlights an r2 value of 0.2, and the red line describes the overall r2 decay according to the Hill and Weir formula.

Figure 2

Unfiltered SNPs in Clonostachys rosea across all scaffolds, numbered as described in Supplementary Table 3. Manhattan plots of the scaffold positions and P-values of cold growth for all SNPs across the all scaffolds. The red line indicates a P-value of 0.05-log10.

Empirical Bayesian multiple hypothesis testing identified an association between 1,478 SNP markers and normalized growth rate at 10°C with an lfsr lower than 0.05 (Supplementary Table 4). These 1,478 SNPs were located on 227 different scaffolds (Supplementary Table 4). A single SNP (unitig_823:6016) was identified as associated with normalized growth rate at 10°C using FDR correction (Supplementary Table 4).

Gene content in genomic regions associated with growth rate at 10°C

Out of the 1,478 significantly cold growth associated SNP markers, 581 were located within or close (< 1,000 bp) to 265 different genes, mostly predicted to encode proteins with unknown function (Supplementary Table 4). Out of these, 175 SNPs were located in exonic regions, 37 in intronic regions, 2 in splice sites, and 40 in the untranslated regions (UTRs) on the 3′ or 5′ ends. 145 of the 1,478 SNPs were highly significant according to the ashr analysis (lfsr < 10−10). The 20 most significant SNPs based on ashr are depicted in Table 2. The most common gene ontology categories were catalytic activity (GO:0003824), protein binding (GO:0005515), metabolic process (GO:0008152), oxidoreductase activity (GO:0016491), oxidation-reduction process (GO:0055114), ATP binding (GO:0005524), and zink ion binding (GO:0008270). The most highly associated SNP to growth rate at 10°C according to standard PLINK association (P = 1.5*10−7, FDR = 0.01) was located in gene CRV2G00015476. This gene was predicted to encode a beta lactam utilization protein, LamB. Nineteen SNPs were identified in an intergenic region between a gene for a predicted short-chain dehydrogenase (CRV2G00018099) and a gene for a predicted glutathione-dependent formaldehyde-activating enzyme (CRV2G00018100).

Table 2

Top 20 highest significance SNPs reported to be associated with growth rate at 10°C in ashr analysis.

SNP	lfsr	Location	Location	Annotation	NCBI Conserved domain search
unitig_630:6733	0.00	intergenic	CRV2G00019736(dist = 3623)	Protein of unknown function	No conserved domains identified
unitig_878:353	0.00	intergenic	NONE	.	.
unitig_760:385	0.00	intergenic	CRV2G00020684(dist = 6578)	Similar to FRQ Frequency clock protein	FRQ domain
unitig_47:563	0.00	intergenic	NONE	.	.
unitig_47:602	0.00	intergenic	NONE	.	.
unitig_793:2158	0.00	intergenic	CRV2G00016485(dist = 1475)	Protein of unknown function	No conserved domains identified
unitig_793:2167	0.00	intergenic	CRV2G00016485(dist = 1484)	Protein of unknown function	No conserved domains identified
unitig_298:2434	0.00	intergenic	NONE	.	.
unitig_298:2460	0.00	intergenic	NONE	.	.
unitig_280:2479	0.00	intergenic	NONE	.	.
unitig_345:10970	0.00	intergenic	CRV2G00021787(dist = 1241)	Protein of unknown function	FN3-like domain
unitig_246:11272	0.00	exonic	CRV2G00021642	Similar to fmpE Nonribosomal peptide synthetase	AFD class I superfamily, alpha am amid superfamily
unitig_640:12334	0.00	exonic	CRV2G00022020	Protein of unknown function	Fungal_TF_MHR, GAL4, DNA pol3 gamma3 superfamily
unitig_640:12340	0.00	exonic	CRV2G00022020	Protein of unknown function	Fungal_TF_MHR, GAL4, DNA pol3 gamma3 superfamily
unitig_818:85374	0.00	intergenic	CRV2G00017605(dist = 1606)	Mitochondrial chaperone BCS1-B	UhpC superfamily
unitig_601:93320	0.00	intergenic	CRV2G00019351(dist = 7767)	Protein of unknown function	Fasciclin domain
unitig_643:109199	0.00	intergenic	CRV2G00019327(dist = 50251)	Similar to betA Oxygen-dependent choline dehydrogenase	Choline dehydrogenase
unitig_697:153157	0.00	upstream	CRV2G00018410(dist = 872)	Protein of unknown function	No conserved domains identified
unitig_823:628705	0.00	intergenic	CRV2G00015724(dist = 3408)	Similar to FUS6 Efflux pump FUS6	TRI12 Superfamily
scf_014:970797	0.00	intergenic	CRV2G00012293(dist = 2794)	Protein of unknown function	GH18 chitinase-like superfamily

Ten and eight SNPs, respectively, were identified close, exonic or within 1,000 bp up- or downstream, to two genes predicted to encode mitochondrial chaperone proteins (CRV2G00019610 and CRV2G00008672) (Supplementary Table 4). Two additional genes encoding heat shock protein (HSP) chaperones contained exonic SNPs; the small ATP-independent HSP chaperone CRV2G00019134 and the ATP-dependent HSP70 chaperone CRV2G00021296. Four SNPs were identified close, exonic, or within 500 bp upstream, to a gene for a predicted efflux membrane protein (CRV2G00021344). A more detailed analysis revealed that is was a major facilitator superfamily (MFS) transporter from the drug:H+ antiporter-2 family (2.A.1.3), with the highest sequence similarity to members from the 2.A.1.3.33, 2.A.1.3.65, and 2.A.1.3.73 subfamilies that are implicated in drug, fungicide, or fungal secondary metabolite export. We also identified SNPs nearby (< 700 bp) two additional genes putatively encoding MFS efflux pumps from the same drug:H+ antiporter-2 family (CRV2G00015724 and CRV2G00018834), also implicated in drug resistance. Other transport protein-encoding genes associated with significant SNPs were CRV2G00019382 and CRV2G00018067 encoding homologs of Saccharomyces cerevisiae FMP42, possibly associated with mitochondria.

We also identified two different genes encoding putative esterases/lipases (CRV2G00019392, CRV2G00021710) with three and one SNPs, respectively, exonic or within 600 bp downstream, significantly associated with growth rate at cold temperature. We detected seven SNPs close to or within three genes associated with chitin metabolism; CRV2G00008512, CRV2G00019295 (chitin synthases), and CRV2G00019292 (the chiA4 chitinase). Additionally, we identified CRV2G00012293 (chitinase domain containing protein) with a significantly cold growth associated SNP at a distance of 2794 bp (Table 2). A sequence analysis of CRV2G00008512 and CRV2G00019295 showed that they displayed similarity with chitin synthases from other Sordariomycete fungi (Supplementary Table 5). A phylogenetic analysis further revealed that CRV2G00019295 clustered together with chitin synthase 1 from N. crassa, while CRV2G00008512 was more distantly related with chitin synthases from other Sordariomycetes (Figure 3). Besides the two previously mentioned mitochondrial chaperone protein genes and the FMP42 transporters, additional genes in close proximity with significant SNPs were predicted to encode proteins related to mitochondrial functioning; CRV2G00021544 putatively encoding a mitochondrial kinase activator of the bc1 complex required for the biosynthesis of Coenzyme Q and CRV2G00019380 similar to a mitochondrial aminolevulinate synthase involved in heme biosynthesis (Supplementary Table 5). We also detected three SNPs associated (one 6578 bp upstream, two in the 5′ UTR region) with CRV2G00020684, and one SNP associated with CRV2T00007803 (exonic), associated with frequency clock protein (FRQ) encoding gene in C. rosea (Supplementary Table 4, Table 2).

Figure 3

Phylogenetic analysis of Clonostachys rosea chitin synthases. Amino acid sequences of selected chitin synthases were aligned using MUSCLE and a phylogenetic analysis conducted using maximum likelihood methods. Bootstrapped branch statistics values are depicted above the branches. The bar indicates mean number of amino acid substitutions per site. Protein identifiers include NCBI Gene ID number.

Discussion

Fungi exhibit several advantages for applying GWA analysis such as small genomes, a haploid life stage and production of high numbers of clonal propagules that can be used for repeated phenotypic assessments. Especially the possibility to genotype haploid individuals by sequencing makes GWAS particularly powerful, since the phenotype is the results of a single haplotype and not a result of two haplotypes as in the diploid case where heterozygosity can confound the signal. Despite these advantages, GWA analyses have not yet been utilized to study the genetic basis for biocontrol traits in fungal BCAs. In this study, we therefore set out to establish C. rosea as a population genomic model system for BCA fungi and to evaluate its usefulness for GWAS.

A high-quality reference genome with few gaps is an important aspect for a successful GWAS. The PacBio sequencing platform has been shown to improve older assemblies, primarily for eukaryotic organisms, by reconciling the long reads with older second generation short read sequencers (Chakraborty et al., 2016; Jayakumar and Sakakibara, 2017). Combining PacBio polished contigs with Illumina-based sequences substantially increased N50 and decreased gaps in the C. rosea ver. 2 genome. It also increased the total genome size to 70.7 Mbp, which is considerably larger than the genomes of other hypocrealean mycoparasitic species (Karlsson et al., 2018) such as Trichoderma spp. (31.7–39.0 Mbp, Kubicek et al., 2011; Studholme et al., 2013; Xie et al., 2014), Tolypocladium ophioglossoides (31.2 Mbp, Quandt et al., 2015) and Escovopsis weberi (29.5 Mbp, de Man et al., 2016). The draft genome of the closely related species C. chloroleuca (Moreira et al., 2016) is 55.4 Mbp (Sun et al., 2015), although this was based on Illumina-technology only and therefore comparable to the C. rosea ver. 1 genome of 58.3 Mbp (Karlsson et al., 2015). The large genomes of Clonostachys spp. is however comparable to certain related Fusarium spp., such as F. solani (51.1 Mbp, Coleman et al., 2009) and F. oxysporum (61.4 Mbp, Ma et al., 2010). The C. rosea ver. 2 genome was predicted to contain 21,246 genes that were 6,978 more than ver. 1 (Karlsson et al., 2015). However, 3,738 of these genes were solely predicted through ab initio methods and some of these may represent artifacts, suggesting that the total gene number in C. rosea is between 17,508 and 21,246.

We further established a population genomic dataset by sequencing the genomes of 52 different C. rosea strains. The high coverage of the Illumina-based genomes and the high sequence similarity with the strain IK726 ver. 2 genome indicate high quality of the generated genomes. Previous fungal GWASs have used both smaller and larger populations to successfully associate SNP markers with virulence in H. annosum (23 strains), Z. tritici (106 strains), P. nodorum (191 strains), with virulence, deoxynivaleol production and azole sensitivity in F. graminearum (119 strains) and with azole sensitivity in Rhynchosporium commune (120 strains) (Dalman et al., 2013; Gao et al., 2016; Mohd-Assaad et al., 2016; Talas et al., 2016; Hartmann et al., 2017).

Another requirement of successful GWA analyses is the presence of enough recombination events in the population. The homothallic reproductive mode of C. rosea (Schroers et al., 1999) may be an obstacle for this type of analysis if repeated selfings between identical genomes results in an essentially clonal population. However, the high genetic variation and rapid LD decay detected in our worldwide C. rosea population is a strong indication of reoccurring outcrossing between genetically different strains, showing that C. rosea is a facultative homothallic species. Genetically distinct C. rosea individuals have been isolated from the same field in Denmark (Bulat et al., 2000), indicating local genetic variation. Another homothallic species, Aspergillus nidulans, has been shown to occasionally outcross, resulting in genetic recombination within populations (Kronstad, 2007; López-Villavicencio et al., 2013). Recombination and possible outcrossing has also been reported for Aspergillus fumigatus and Letharia species (Kroken and Taylor, 2001; Pringle et al., 2005).

High variation in growth rate at 10°C in C. rosea, with no connection with the geographic origin of the strains and the lack of genetic population structure indicate that a GWAS of the trait is possible. High growth rate at lower temperature may be an important trait for BCAs, especially for rapid colonization of new substrates with limited amounts of competing microorganisms such as wounds, flowers, and developing roots (Jensen et al., 2016). The connection between cold tolerance and interference competition through mycoparasitism or antibiosis (Jensen et al., 2017) is less well established. In fact, C. rosea strain IK726 was among the most slow-growing strains at 10°C in our study, at the same time as it repeatedly has been shown to be a very efficient BCA against several different fungal plant pathogenic species (Jensen et al., 2007).

The normal distribution of the cold tolerance trait in C. rosea further suggests a polygenic inheritance that depends on several different genes. Many medically and agriculturally relevant traits are in fact multigenic and it has been notoriously difficult to achieve a reliable estimation of the true effect size of individual loci to genetic variation contributing to polygenic phenotypic differences (Boyle et al., 2017). Thus, the emerging view of the genome as a uniform distribution of contributing, however small, sequence variants to different phenotypes, the “omnigenic” model, highlights the propensity to underestimate the true significance of SNPs and similar genomic variations on a trait (Boyle et al., 2017). So, in line with the omnigenic trait view of the genome, we used the R package ashr that produces SNPs with a true effect size of zero or non-zero (Stephens, 2016). Recently, ashr has been used in some GWASs as an empirical Bayes large scale hypothesis testing alternative to false discovery rate analysis, by using effect size and its standard error, instead of using only one parameter, such as P-values (Stephens, 2016; Boyle et al., 2017; Petit et al., 2017). The approach by Stephens (2016) has been proposed to be more robust when dealing with large scale hypothesis testing, producing a local false sign rate (lfsr) as an analogous alternative to false discovery rate. The use of ashr allowed us to identify 1,478 SNPs significantly associated with C. rosea growth rate at 10°C, located physically (< 1000 bp) close to about 265 genes.

As C. rosea is found globally, an evolved ability for cold tolerance would be integral for a generalist lifestyle across climates (Knudsen et al., 1995; Robinson, 2001; Abreu et al., 2014; Karlsson et al., 2015; Moreira et al., 2016). One important reason for wanting to study cold tolerance refers to the survival in temperate and arctic regions, as well as, in cold storage spaces, imperative for the long-term practical function and fitness of BCAs. The mechanisms for adapting to cold are numerous, and may include secretion of proteins with high catalytic efficiency active at low temperatures, transporters, lipid metabolism, secretion of secondary metabolites that counteract freezing damage of the membrane and intracellularly (Robinson, 2001; Blanc et al., 2012). Temperature and day length, depending on the geographical location, have been suggested to be major environmental factors impacting local adaptation of diverging N. crassa populations (Ellison et al., 2011). Significant SNPs were identified in or close to three genes encoding putative chaperone proteins. Regulation of mitochondrial chaperone proteins has been identified in association with temperature adaptation in diverse organisms such as the moth Epiblema scudderiana and the fungus S. cerevisiae (Schmitt, 1996; Lyons et al., 2015; Zhang et al., 2018). As chaperone proteins have also been discussed as important for cold tolerance in plants as well, this result demonstrates an interesting basis for further studies of chaperone importance in cold stress in C. rosea (Al-Whaibi, 2011). Several other genes putatively encoding proteins related to the mitochondria were identified by the GWAS, and further emphasize the importance of mitochondrial functioning for abiotic stress tolerance. The importance of mitochondrial activity in response to cold stress has been discussed in both animals and plants (Mollica et al., 2005; Fangue et al., 2009; Karami-Moalem et al., 2018).

Chitin is an important component of the fungal cell wall, and balancing chitin synthesis and degradation provides the appropriate plasticity and rigidity needed for hyphal growth (Specht et al., 1996). Therefore, it was interesting to find two putative chitin synthase genes and the chiA4 chitinase gene among the genes associated with cold tolerance. The chiA4 gene was reported to be constitutively expressed in C. rosea (Tzelepis et al., 2015), which fits well with a role in cell wall modeling during growth. Chitinases have also been implicated in cold tolerance in plants (Nakamura et al., 2008; Kashyap and Deswal, 2017). Certain fungal chitin synthases have also been associated with adaptation to temperature stress (Lenardon et al., 2009; Liu et al., 2017). The chitin synthase encoded by CRV2G00019295 was found to be most related to chs-1 in N. crassa, which has been found to play a major role in cell wall biogenesis (Yarden and Yanofsky, 1991).

Three different putative membrane efflux transporter genes from the MFS drug:H+ antiporter-2 family also contained SNPs associated with cold tolerance. Both the MFS transporter families 2.A.1.3.33 and 2.A.1.3.65 are shown to evolve under selection for increased gene copy numbers in C. rosea, emphasizing the importance for membrane transport in the lifestyle of C. rosea (Nygren et al., 2018). The connection between membrane transport and cold tolerance is less established. Efflux pumps have been associated with cold stress adaptation in bacteria, specifically by Moraxella catarrhalis, where cold stress reportedly activates transcription and efflux activity the AcrAB-OprM system (Spaniol et al., 2015). Furthermore, studies indicate a connection between ions, secondary metabolite production (such as flavonoids) and membrane transport with cold tolerance in animals and plants (Kaplan, 2004; Wei et al., 2006; Storey and Storey, 2012; Wang M. et al., 2013; Wang X. C. et al., 2013). Additionally, secondary metabolites such as melanin have proven to be important in adjustment to temperature for the fungus R. commune and Cryptococcus neoformans (Rosas and Casadevall, 1997; Zhu et al., 2018). Thus, differences in transporter effectiveness due to differing alleles between strains may impact the ability of fungi to adjust their osmotic equilibrium, concentration of important metabolites for protein stability and membrane stability.

Finally, indications of lipid modifications in C. rosea as a response to colder temperatures comes from the association of SNPs situated in several putative lipases. Lipid structure and concentration is important for cold adaptation in many organisms (Arthur and Watson, 1976; Upchurch, 2008; van Dooremalen et al., 2011; Treseder and Lennon, 2015) and further studies of the association between the lipases and cold tolerance C. rosea is warranted. Comparative genomics and a growth assay at 10°C of N. crassa has suggested that the FRQ protein may be of significance in cold resistance, in part due to the link between lower temperatures and less daylight hours across geographical location (Ellison et al., 2011). From the GWAS we detected four significant SNPs associated with two genes in C. rosea encoding predicted FRQ proteins with some homology (~80–90% cover, 33–36% identities) to N. crassa FRQ; CRV2G00020684 (one SNP 6578 bp upstream, two in the UTR5 region) and CRV2G00007803 (one exonic SNP). These results support the link between growth in cold and FRQ (Table 2, Supplementary Table 4).

We here present an improved genome assembly of C. rosea strain IK726 and whole-genome re-sequencing of an additional 52 C. rosea strains, thereby establishing C. rosea as the premier model species for genomic investigations of fungal BCAs. We further show that GWAS is a powerful tool for studying the genetic basis for phenotypes related with biocontrol activity. We identified an association between cold tolerance and SNPs residing in multiple genes, encoding putative mitochondrial chaperone and abiotic stress proteins, membrane efflux transport proteins, and lipases. These genes are of interest for further studies of cold stress and potential biocontrol efficiency in temperate climates.

Author contributions

All authors were involved in planning the study and writing the manuscript. M-HS, S-DL, and H-JS performed identification of the strains, KI and MD performed fungal growth and maintenance, and KI performed DNA extraction and preparation. MB and MK performed the growth assays. MB, MBD, and MK performed the genome sequencing analysis and GWAS analysis. DJ, MBD, and MK secured the funding.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.