Genomic Characterization of the Mycoparasite Pestalotiopsis sp. Strain cr013 from Cronartium ribicola

The aeciospores of a fungus infecting Pinus armandii were collected in Kunming, Yunnan Province, China, in August 2013. The aeciospores were incubated on filter paper soaked with distilled water until colonies or mycelia appeared at 25°C. After seven days of cultivation, a strain was isolated from the aeciospores and identified as Pestalotiopsis sp. (strain cr013). The conidia of the strain cr013 were inoculated into an improved Fries liquid medium (1 g KH₂PO₄, 0.5 g MgSO₄ · 7H₂O, 0.1 g NaCl, 0.13 g CaCl₂ · 2H₂O, 20.0 g sucrose, 5.0 g ammonium tartrate, 1.0 g yeast extract, 1.0 g NH₄NO₃, 1 l distilled water; natural pH). After 60 hours of culture at room temperature at 150 r/min, the strain cr013 mycelia were filtered and collected, washed with sterile water three times, stored in liquid nitrogen, and sent to BGI Technology Co., Ltd. (China) for follow-up work.

Pestalotiopsis sp. strain cr013 DNA was isolated using the bacterial genomic DNA Extraction Kit (Tiangen Biotech(Beijing)Co., Ltd., China) and cut into 50–800 bp fragments with a Covaris® E220 Focused-ultrasonicator (Covaris, LLC., UK). High-quality 150–250 bp DNA fragments were screened using AMPure XP beads (Beckman Coulter, Inc., USA) and repaired with T4 DNA polymerase (Enzymatics, Inc., USA). The ends of selected fragments were ligated with T-tailed adapters and magnified over eight cycles by KAPA HiFi Hot-Start ReadyMix (Kapa Biosystems, Inc., USA). Then, the amplification products were subjected to singlestrand circularization using T4 DNA ligase to generate a single-stranded circular DNA library. Subsequently, T4 DNA ligase was used to cyclize the magnified products to create a single-stranded circular DNA library.

The WGS library was sequenced using the BGISEQ-500 platform (BGI-Qingdao, China), and the original 100 bp double-terminal sequencing data were obtained. The raw reads with a high ratio of N (ambiguous) bases and low-mass bases were removed using SOAPnuke (v1.6.5) (Chen et al. 2018) with the following parameters: “-l 15 -q 0.2 -n 0.05 -Q 2 -c 0”, to obtain clean data. Genomic DNA was quantified by dsDNA BR analysis, and DNA integrity was evaluated by pulsed-field gel electrophoresis (PFGE). To obtain a long Oxford Nanopore fragment, BluePippin (SAGE Science, Inc., USA) was used to select fragments of approximately five μg of genomic DNA. Library preparation and sequencing were performed using Ligation Sequencing Kits SQK-LSK109 (Oxford Nanopore Technologies, UK). DNA repair, end repair, and A-tailing were carried out with NEBNext® FFPE DNA Repair Mix (M6630; New England Biolabs® Inc., USA) and NEBNext® Ultra™ II End Repair/dA-tailing Module (E7546; New England Biolabs^® Inc., USA). The DNA was then purified with AMPure XP (A63882; Beckman Coulter Inc., USA). DNA was linked to the adaptor by NEBNext® Quick T4 DNA Ligase (E6056; New England Biolabs® Inc., USA). According to the instructions, samples were loaded on FLO-MIN 106D R9.4.1 and sequenced on GridION for 48 hours. Prior to assembly, the size of the genome, heterozygosity, and repeatability were estimated by k-mer analysis according to the second-generation data. The nanopore assembly of the clean BGI-seq 500 WGS reads was conducted with Canu (Koren et al. 2017), which is larger than two kb with the following parameters: “useGrid = false maxThreads = 30 maxMemory = 60 g -nanopore-raw *.fastq -p -d”. Then, we used Pilon (Walker et al. 2014) to correct base errors with the second-generation data and obtain the final assembly results. The assembly’s confidence with sordariomyceta_odb9 was assessed using BUSCO (v3.0.1) (Simão et al. 2015).

Multiple tools were used to identify the repetitive elements. Transposable elements (TE) were recognized by alignment against the Repbase (Bao et al. 2015) database using RepeatMasker (v 4.0.5) (Tarailo-Graovac and Chen 2009) with the following parameters: “-nolow -no_is -norna -engine wublast”, and RepeatProteinMasker (v 4.0.5) with the following parameters: “-noLowSimple -pvalue 0.0001”, at the DNA and protein levels, respectively. The de novo repeat library was evaluated using RepeatModeler (v1.0.8) and LTR-FINDER (v1.0.6) (Xu and Wang 2007) with the default parameters. For the de novo-identified repeats, repeat sequences were classified using RepeatMasker (v4.0.5) with the same parameters. In addition, the tandem repeats were identified using Tandem Repeat Finder (v4.04) (Benson 1999) with the following parameters: “2 7 7 80 10 50 000 -d -h”.

For the prediction of homologous genes, the protein sequences of P. fici W106-1 (GCA_000516985.1) and Pestalotiopsis sp. JCM 9685 (GCA_001599175.1) were downloaded from NCBI by homology-based prediction. The genomes were annotated by homologue using Gene-Wise (http://www.ebi.ac.uk/Tools/psa/genewise) with default parameters. For de novo prediction, Augustus (http://bioinf.uni-greifswald.de/augustus) and GeneMark (http://exon.gatech.edu/GeneMark) softwares were used with the parameters “--ES – fungus”. Finally, EVM (http://evidencemodeler.sourceforge.net) integration software was used to integrate different annotation results. Among them, the weight of homologous annotation results was set as 10, and the weight of prediction results of two de novo software programs was set as 1. To evaluate the integrity of the predicted gene set, BUSCO software and the sordariomyceta_Odb9 database were used to evaluate the integrity of our assembled protein-coding genes.

The predicted genes were aligned to the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al. 2004), SwissProt (Magrane and UniProt Consortium 2011), Cluster of Orthologous Groups of proteins (COG) (Tatusov et al. 2003), CAZyme (Cantarel et al. 2009), NR database, Pathogen Host Interactions (PHI) and NR using BLAST (v2.2.26) with the parameters “-p BLASTP -e 1e-5 -F F -a 4 -m 8”, E-value < 1 × 10^-5. Additionally, the motifs and domains of genes were identified by protein databases, including Pfam, SMART, PANTHER, PRINTS, ProSite, and ProDom using InterProScan (5.16–55.0) (http://www.ebi.ac.uk/Tools/InterProScan), and Gene Ontology (GO) (Ashburner et al. 2000) databases using BLAST (v2.2.26) with the parameters “-p BLASTP -e 1e-5 -F F -a 4 -m 8”, E-value < 1 × 10^-5. The cr013 strain assembly was uploaded to the antiSMASH (v5.0) (Blin et al. 2019) website to identify the secondary metabolite gene clusters.

The mechanism of mycoparasites mainly includes the production of cell wall degrading enzymes and small molecular toxins, and these processes involve CAZyme, cytochrome P450, protease, and Zn2Cys6 transcription factor. Therefore, to comparatively analyze these genes in strain cr013, T. atroviride, T. virens, T. reesei and T. harzianum, BLASTP and Hmm (http://hmmer.janelia.org) software were used to identify candidate gene families. First, the protein sequences of four other Trichoderma species were downloaded from the NCBI database and aligned against the protein sequences of five species using BLASTP v2.2.31 with an E-value of 1 × 10^-5 to search for significant hits. Additionally, Hmmsearch was used to predict gene families in conjunction with the Hmm model from the Pfam database. The protein sequences of each gene family were aligned using MUSCLE v3.8.1551 (Edgar 2004) if the gene family lacked the Hmm model, and hmm models were built using Hmmbuild. The significant hits (E-value of 1 × 10^-5) from BLASTP and Hmmsearch were dereplicated, filtered, and merged to generate final hits. The sequences without domains of the resulting hits were removed using the Conserved Domain Database (CDD) (Lu et al. 2020) and Simple Modular Architecture Research Tool (SMART) (Letunic and Bork 2018; Letunic et al. 2020) databases. Furthermore, the sequences coding fewer than 80 amino acids were discarded, and the longest transcript was reserved because multiple transcripts were identified.

The concentration and purity of the strain cr013 were determined using a Qubit™ 4 Fluorometer (Invitrogen™, USA) and a NanoDrop™ Spectrophotometer (Thermo Fisher Scientific Inc., USA) and then assessed by 1% agarose gel electrophoresis. DNA integrity was assessed by gel electrophoresis (Fig. 1). According to the Quality Standard of Nanopore Genome Sequencing Samples, the quality of the samples met the requirements for database building, and subsequent genome sequencing could be carried out.

Fig. 1.

The Nanopore platform was used to sequence the long fragment of the strain cr013, and a total of 11.10 G of raw data was generated. Prior to assembly, the K-mer was selected as 15, and k-mer analysis was performed according to the second-generation data to evaluate the genome size, heterozygosity, and repeatability. The filtered data were processed using Jellyfish software. The results indicate that the genome size of strain cr013 was 47.4 Mb, the genome repeatability was 18.5%, and the heterozygosity was 0.06%. Nanopore data were assembled using Canu, and the second-generation data were used for basic error correction using Pilon to obtain a fine genome map. A total of 13 scaffolds were assembled into the genome with a size of 46.69 Mb. BUSCO was used to evaluate the integrity of the genome assembly. The results showed that 97.2% of the core genes could be annotated in the genome. Complete and single-copy BUSCOs were 96.80%, complete and duplicated BUSCOs accounted for 0.4%, fragmented BUSCOs accounted for 1%, missing BUSCOs accounted for 1.8%, and 3,725 BUSCO groups, which reflected the high integrity of the assembly results, were searched in total.

The Pestalotiopsis sp. strain cr013 mycoparasite genome (GenBank accession number: JACFXT010000000, BioProject ID: PRJNA647543, BioSample ID: SAMN15589943) was assembled into 13 scaffolds with an N50 of 6.39 Mb, total size of 46.69 Mb, and a 52.1% GC content. A total of 14,893 genes were predicted. Through BUSCO and sordariomyceta_Odb9 database analysis, 96.5% of the genes were predicted to be complete, reflecting the high annotation integrity of the annotated gene set. Regarding repetitive sequences, the cr013 strain showed a 2.07% repeat content and 1.68% transposable element (TE) content (Table I).

After obtaining the genome of the cr013 strain, we compared and annotated the gene database to determine the functions and related descriptions of the genes. The results reflect the functional classification of the gene set of the strain as a whole and facilitate subsequent research to identify target functional genes. In this study, the predicted gene set was aligned to NR, KEGG, COG, SwissProt, GO, ARDB, PHI, and CAZyme using BLASTP (E-value ≤ 1 × 10^-5). The results are shown in Table II. Finally, approximately 94.10% of the protein sequences were aligned to eight databases.

In KEGG annotations, genes are usually classified into five branches based on the following KEGG metabolic pathway categories: cell processes, environmental information processing, genetic information processing, metabolism, and organic systems. Among the five pathways annotated through the KEGG analysis of the strain cr013, metabolic pathways were the most enriched category, mainly related to the metabolism of carbohydrates, amino acids, lipids, coenzymes, and vitamins (Fig. 2). The analysis of KEGG enrichment showed that 9,768 genes corresponding to KEGG pathways were enriched in 129 metabolic pathways. Five branches of metabolic pathways involved the most genes (Table III).

Fig. 2.

According to the COG classification of the strain cr013 genes obtained by sequencing, 2,812 genes were classified into 24 categories (Fig. 3), except for general function prediction. The most involved genes were posttranslational modification, protein turnover, chaperones (195 genes, 6.93%), signal transduction mechanisms (171 genes, 6.08%), translation, ribosomal structure and biogenesis (228 genes, 8.1%), amino acid transport and metabolism (271 genes, 9.63%), carbohydrate transport and metabolism (306 genes, 10.9%), coenzyme transport and metabolism (206 genes, 7.32%), energy production and conversion (327 genes, 11.62%), lipid transport and metabolism (376 genes, 13.37%), and the biosynthesis, transport and catabolism of secondary metabolites (319 genes, 11.34%).

Fig. 3.

A total of 8,924 genes were annotated in the GO database, which can be divided into three categories with 52 components (Fig. 4): biological processes (25 branches), cellular components (14 branches), and molecular functions (13 branches). The most significant number of genes in the biological process category were associated with metallic processes and cellular processes; the greatest number of genes in the cellular component category were associated with membranes, membrane parts, cells, cell parts, and organelles; and the greatest number of genes in the molecular function category were related to quantitative activity and binding.

Fig. 4.

CAZyme is a classic database relevant to carbohydrate-active enzymes, which include many enzymes that can catalyze the degradation, modification, and biosynthesis of carbohydrates. Among the 709 genes annotated in strain cr013, 255 were glycoside hydrolases (GHs), 92 were glycosyl transferases (GTs), 17 were polysaccharide lyases (PLs), 41 were carbohydrate esterases (CEs), 170 were carbohydrate-binding modules (CBMs), and 134 were auxiliary activities (AASs).

Fungi can produce many carbohydrate-active enzymes, among which chitinase and 1, 3-β-glucan enzymes are important hydrolase enzymes (Zhao et al. 2013). Therefore, the GH families of Trichoderma spp. and the strain cr013 were further analyzed. There were three GH18 and four GH19 genes in the chitinase family of the cr013 strain. The number of chitinase genes in the cr013 strain was significantly lower than that in Trichoderma sp. mycoparasite. The number of 1, 3-β-glucans (GH17, GH55, GH64 and GH81 families) in the strain cr013 was 20, more significant than the number in the selected Trichoderma sp. mycoparasite. Cuomo et al. (2007) found that the number of GH75 family genes (chitosanase) increased significantly during cell wall degradation. The number of GH75 family genes in the strain cr013was similar to that in the Trichoderma sp. mycoparasite (Table IV). The results showed that there were significantly fewer hydrolytic enzymes related to cell wall degradation in the strain cr013 than in the selected Trichoderma sp. mycoparasite, which could explain the difference in mycoparasitism mechanism between Trichoderma and Pestalotiopsis.

The secondary metabolite gene clusters of strain cr013 and other selected fungi were analyzed using antiSMASH with default parameters, and the results are shown in Table V. Two ambuic acids and four reducing polyketides with novel structures that had strong antitumor activity were isolated from the mycoparasitic strain cr013. Polyketide compounds were found for the first time in the mycoparasitic Pestalotiopsis sp. Therefore, the biological synthesis mechanism of these compounds was explored by genome sequencing, and more secondary metabolites were expected to be found. Further analysis results by antiSMASH showed that the genes of the strain cr013 were distributed in 64 gene clusters, including six terpenes, nine NRPSs, 15 similar NRPSs, 19 T1 PKSs, one T3 PKS, two indoles, six hybrid NRPS-PKSs, one betalactone, two hybrid NRPS-terpenes, one RiPP, and two hybrid indole-PKSs. Compared to the selected fungi, the number of gene clusters of secondary metabolites of the strain cr013 was much greater than Trichoderma sp. mycoparasite. The number of gene clusters of secondary metabolites was less than endophytic Pestalotiopsis sp. However, the strain cr013 contained hybrid NRPS-terpenes, which was not predicted in the Pestalotiopsis sp. reference genome. According to the results analysed by antiSMASH, the Pestalotiopsis strain cr013 had the potential to produce abundant secondary metabolites.

Cytochrome P450 participates in many important cellular processes (van den Brink et al. 1998) and complex biotransformation processes in fungi, and the enzymes of this group catalyze the conversion of hydrophobic intermediates of primary and secondary metabolic pathways (Crešnar and Petrič 2011). There are 234 cytochrome P450-coding bases and 74 protease-coding bases in the strain cr013 genome, which is a significantly higher number than in the selected Trichoderma sp. genome. A total of to 152 transcription factors were found in the genome sequence analysis, among which the 13 genes encoding C2H2 transcription factors, and three genes encoding Zn2cys6 transcription factors represented a significantly smaller number than in the Trichoderma sp. genome (Table VI).