Penicillium rubens strain 212 and other P. rubens strains, isolated from diverse agricultural soils and plant samples in Spain, are listed in Table 1. Conidia from these strains were stored in 20% glycerol for a long term at −20°C except for PO212, which was stored at 4°C as dried conidia. Dried conidia of PO212 were produced in a solid-state fermentation system and dried as previously described by Larena et al. (2003a). Penicillium strains were grown on potato dextrose agar (PDA) or minimal medium (MM; Espeso et al., 2019) with 5 mM ammonium tartrate and D-glucose 1% (w/v) as nitrogen and carbon source, respectively, and incubated at 25°C for 5 days. For short-term storage, strains were kept at 4°C on solid media.

The biocontrol activity of P. rubens strains was tested using the pathogenic isolate 1A of FOL, provided by Dr. Cristina Moyano from the Laboratory for Assessment of Variety, Seed and Nursery Plants, INIA-CSIC (Madrid, Spain). FOL was stored at 4°C in tubes containing sterile sand. For mycelial production, conidia from FOL stored in sterile sand were germinated on Czapek Dox agar (CDA; Difco Laboratories, Detroit, MI, United States) and cultivated in darkness at 25°C for 7 days. Microconidial inoculum of FOL was produced in 250 ml flasks containing 150 ml of sterile Czapek Dox broth (Difco). Each flask was inoculated with three mycelial plugs (1 cm diameter) from the 7-day-old cultures on CDA (De Cal et al., 1995) and incubated for 5 days at 25°C in a rotary shaker (model 3527; Lab-Line Instruments, Inc.) at 150 rpm. Microconidia were separated from the mycelial mass by filtration through glass wool. The conidial concentration was determined using a hematocytometer and adjusted to 10⁶ microconidia/ml.

At least two growth chamber experiments were carried out on tomato plants to evaluate the biocontrol efficacy of S27 against FOL as described by Villarino et al. (2018). Seeds of tomato cultivar “San Pedro”, which is susceptible to races 1 and 2 of FOL, were used in all experiments. Tomato seeds were sown in sterile trays (27 cm × 42 cm × 7 cm) that contained an autoclaved (for 1 h at pressure of 1 kg cm⁻² and temperature of 121°C, during 3 consecutive days) mixture of vermiculite (Termita; Asfaltex, S.A., Barcelona, Spain) and peat (Gebr. BRILL substrate GmbH & Co. KG; 1:1, v:v). The trays were maintained in a growth chamber at 25°C with fluorescent lighting (100 μE m⁻² s⁻¹, 16-h photoperiod) and 80–100% relative humidity for 3–4 weeks. Tomato seedlings (with at least two true leaves) were treated 7 days before transplanting with an aqueous conidial suspension (6 × 10⁶ conidia per gram of substrate) of S27 or PO212. Conidial suspensions of PO212 and S27 were prepared as follows: Dried conidia of PO212 were rehydrated in sterile distilled water (SDW) using a rotatory shaker at 150 rpm for 2 h (CERTOMAT® RM). Conidia of S27 were harvested from colonies grown on PDA and incubated in the dark at 25°C for 7 days. The day before treatment, the viability of PO212 and S27 conidia were estimated by measuring their germination as previously described (Larena et al., 2003a). For each replicate (three by sample type), the germination of 50 randomly selected conidia was counted and viability was calculated and expressed as a percentage (De Cal et al., 1990). Seven days after treatment, tomato seedlings were transplanted from seedbeds into 100-ml flasks containing 100 ml of Hoagland solution (Hoagland and Arnon, 1950) so that the roots were in contact with the solution, as described by De Cal et al. (1997a). An aliquot of an aqueous (SDW) conidial suspension of FOL was added to the flasks just before transplanting so that the final conidial concentration in the flasks was 1 × 10⁵ conidia/ml. Plants that had been inoculated with the pathogen but not treated with any strain of P. rubens, were used as the control. Five replicate flasks, each containing four plants, were used per treatment. The flasks were placed in a randomized complete block design in a growth chamber for 4 weeks under the conditions described earlier in this subsection. The complete experiment was done two times. Disease severity was graded on days 7, 14, 21, and 28 after the transplant. Disease severity followed a 1–5 index scale: 1, healthy plants (0–24%); 2, yellow lower leaves (25–49%); 3, dead lower leaves and some yellow upper leaves (50–74%); 4, dead lower leaves and wilted upper leaves (75–99%); and 5, dead plants (100%; De Cal et al., 1995). All plant roots were placed in humidity chambers at the end of each experiment, and the presence or absence of the pathogen in the crown after 5 days of incubation at 25°C was recorded.

Data of disease severity and incidence were analyzed by ANOVA with the STATGRAPHICS program (XVII Centurion. v. 17.2.00). When the F test was significant at a value of p of < 0.05, means were compared using the Student–Newman–Keul’s multiple range test (Snedecor and Cochran, 1980).

Genomic DNA (gDNA) was extracted from the mycelia of Penicillium strains (Table 1) grown in liquid MM at 25°C for 2 days. Mycelia were harvested by filtration using Miracloth (Calbiochem, Merck-Millipore, Darmstadt, Germany). Samples were lyophilized for at least 6 h. For each sample, mycelium was pulverized using a ceramic bead in a FastPrep-24 homogenizer (MP Biomedicals™), one pulse for 20 s at minimum speed. A sample of 100 mg of powdered mycelium was mixed with 1 ml of DNA extraction solution (25 mM Tris–HCl pH 8.0, 250 mM sucrose, and 20 mM EDTA pH 8.0). Before incubation for 15 min at 65°C, 100 μl of 10% SDS were added to each sample. Next, proteins and cellular debris were removed by adding 1 ml of Phenol/Chloroform/Isoamyl alcohol mixture per sample and further mixing on a rotary shaker (Rotator Multi Bio RS-24) for 15 min at room temperature (RT). Organic and aqueous phases were separated by centrifugation on a benchtop centrifuge at maximum speed for 5 min. gDNA was precipitated from the aqueous phase by the addition of 1/10 vol of sodium acetate pH 6 and 0.6 vol of 2-propanol, followed by incubation for 15 min at RT and centrifugation for 5 min at RT. To wash the gDNA pellet, 1 ml ethanol (80%) was added, followed by centrifugation at max. Speed for 5 min at RT. After drying the ethanol in samples, pellets were dissolved in 500 μl sterile water, and samples were treated with DNase-free RNase A (5 mg/ml; 37°C, 60 min). gDNA was newly precipitated and washed as described before. Finally, gDNA was dissolved in 200 μl nuclease-free water and stored at −20°C until use.

Sequencing of PO212 gDNA was performed on an Illumina MiSeq 500 cycles at the “FPCM, Fundación Parque Científico de Madrid” sequencing facility using 150 bp and 250 bp paired-end sequencing reads. For S27 gDNA, Stab Vida (Portugal) performed the construction and sequencing of DNA libraries. DNA libraries were sequenced on the lllumina Hiseq 4000 platform, using 150 bp paired-end sequencing reads.

Standard PCR protocols were used: initial denaturation at 95°C for 2 min, 30 cycles of denaturation at 95°C for 10 s, annealing at 55°C for 30 s, and extension at 72°C for 1 min/kb. After 30 cycles, an extension at 72°C for 10 min and storage at 4°C. The polymerase Takara (TaKaRa Taq™) was used for amplification. Oligonucleotides were designed using Vector NTI™ Suite 8 and are listed in Table 2. The oligonucleotides for the amplification of MAT1-1 and MAT1-2 were previously described in Espeso et al. (2019). PCR products were analyzed in 0.8% agarose-TAE electrophoresis and when required, they were purified using the PCR clean-Up kit (Macherey-Nagel), following the manufacturer’s instructions. Sequencing of DNA fragments was done by Stab Vida’s (Portugal) sequencing service.

For sequencing of the PO212 genome, two libraries of 150 and 250 bp paired-end fragments were produced and raw sequencing reads were subjected to quality control using the FastQC program.¹ A Q > 28 along the read length and the k-mers and nucleotide distribution were homogeneous in both libraries. An A5-miseq pipeline was used (Coil et al., 2015) to assemble the PO212 genome. The nuclear reads were extended using FLASH (Magoč and Salzberg, 2011) and were assembled using the A5-miseq pipeline, until the output, named final.scaffolds.fasta, was obtained. The previously aligned mitochondrial reads were assembled to obtain a single mitochondrial contig of ~28 kbp.

De novo assembly of the S27 genome was performed using the A5-miseq pipeline (Coil et al., 2015) and further revised with Gapcloser from the Soapdenovo assembler (Luo et al., 2015). The aligned S27 mitochondrial reads were assembled apart from the genome in 13 contigs of ~34 kbp.

Completeness of the genome assemblies was assessed by the Benchmarking Universal Single-Copy Orthologs (BUSCO) v.2.0.1 software tool (Simão et al., 2015). The assembled genomes were annotated using the MAKER (v.2.31.9) pipeline (Campbell et al., 2014). Before annotation, a species-specific repeat library was constructed using RepeatModeler (v.1.0.8) to mask repeats (Tempel, 2012). Gene models were predicted with AUGUSTUS (Stanke et al., 2006) using two strategies, ab initio gene predictors and training with Aspergillus nidulans as a model genome.

Gene models of the PO212 and S27 genomes were manually curated using the P2niaD18 (GCA_000710275) and PrWis (GCA_000226395) strain proteins available in the EnsemblFungi and the National Center of Biotechnology Information (NCBI) database as evidence for gene prediction using Apollo (Dunn et al., 2019). Predicted proteins were functionally annotated using BLASTp (Altschul et al., 1990) against the non-redundant database of the NCBI and classified using InterProScan and Pfam analysis (Jones et al., 2014).

The comparison of the PO212 genome with genomes from other strains of P. rubens and P. chrysogenum, which are available in the database and listed in Table 3, was performed through the Quast bioinformatics tool (Gurevich et al., 2013). To visualize the genomes of the P. rubens and P. chrysogenum strains, Icarus, a genome visualizer based on the Quast tool, was used (Mikheenko et al., 2016). These genomes are divided into two groups according to the literature: P. chrysogenum (P2niaD18, KF-25, HKF42, NCPC10086, and v1.0) and P. rubens (PrWis and Biourge 1923; Table 3).

Due to the current controversy in the classification of P. chrysogenum and P. rubens, in this study, all these strains have been classified as belonging to the P. chrysogenum or P. rubens clade, based on the barcodes of three genes encoding proteins: β-tubulin (BenA), RNA polymerase II second largest subunit (RPB2), and calmodulin (CaM). The barcodes used were JF909949 (BenA), JX996658 (RPB2), and JX996263 (CaM) to P. rubens and AY495981 (BenA), JN121487 (RPB2), and JX996273 (CaM) to P. chrysogenum (Visagie et al., 2014). The sequences of P. rubens and P. chrysogenum were aligned to the barcodes with the Multiple Sequence Alignment tool of Clustal Omega (EMBL-EBI).

MUMmer4 software package (Marçais et al., 2018) was used for pairwise alignment of PO212 and S27 assemblies by setting PO212 as the reference; the minimum length of a matched group to 20 bp and the distance an alignment extension will attempt to extend to poor scoring regions before yielding to 100 bp. Comparative analysis was performed using Circos tools (Krzywinski et al., 2009). Moreover, CLC Genomics Workbench 12 (QIAGEN Bioinformatics; https://digitalinsights.qiagen.com/products-overview/discovery-insightsportfolio/analysis-and-visualization/qiagen-clc-genomics-workbench/) was used to align the reads of S27 to the assembled genome of PO212. The CLC alignment workflow was used to map the raw reads to a reference genome and to detect variants. Default values were used, except for the minimum coverage and minimum count, which were set to 10 and 2, respectively, to avoid loss of information. Of all the variations found, homozygous variations were selected. Only variations with 100% of frequency were taken into account for manual verification by PCR amplification using specific primers and sequencing. To visualize amino acid changes in genes between PO212 and S27 strains, we used Integrative Genome Viewer (IGV version 2.9.4; Robinson et al., 2011).

A phylogenetic tree was constructed based on nucleotide sequence from assemblies of genomes from PO212, S27, and other strains of P. rubens and P. chrysogenum available in the database and listed in Table 3. For this purpose, Aspergillus nidulans strain FGSC-A4 genome (Wortman et al., 2009) and P. oxalicum strain HP7-1 genome (Zhao et al., 2016) were used as outgroups. The phylogenetic tree was generated using the software JolyTree (Criscuolo, 2019). The algorithm uses an alignment-free distance-based procedure for inferring phylogenetic trees from genome contig sequences. The pairwise dissimilarity between each pair of genomes was estimated using the Mash method (Ondov et al., 2016). Balanced Minimum-Evolution (BME) was used to optimize the evolutionary distances calculated in the first step. Finally, the rate of elementary quartets (REQ) was used for assessing the branch confidence values to construct the final tree (Criscuolo, 2019).

In addition, a tree obtained with the concatenated sequences benA, caM, and RPB2 barcodes was generated using the maximum likelihood method with 1,000 bootstrap replications. The tree was drawn to scale, with branch lengths representing the inferred evolutionary distances. Between 1,879 and 1,931 bp were presented in the final data set. Sequences of P. oxalicum barcodes were used as an outgroup. Multiple sequence alignments were performed using the ClustalW algorithm (Thompson et al., 1994). The Tamura-Nei model and MEGA11 software were used for the phylogenetic analysis of sequence data (Tamura and Nei, 1993; Tamura et al., 2021).

The genome of PO212 was sequenced via paired-end Illumina MiSeq technology, providing two DNA-seq libraries: one paired-end small fragment library and a second obtained from long-range DNA fragments. The resulting assembly was performed de novo and evaluated in terms of N50 and L50 using Quast. Automated assembly and manual sequence verification yielded an estimated PO212 genome size of 29.82 Mb at ~200× coverage. PO212 genome was organized into 65 scaffolds with an N50 scaffold length of 1.88 Mb, comprising the first six longest scaffolds and reaching an N90 scaffold length of 0.38 Mb, corresponding to the first 18 scaffolds (Table 4 and Supplementary Table S1). PO212 genome had a GC content of 49.07%, similar to that of other Penicillium genomes. After the manual curation of automated ORF predictions by AUGUSTUS, we determined at least 10,164 ORFs in the PO212 genome.

The data presented in the study are deposited in the GenBank repository, accession number JAPDLE000000000.

To perform a comparison of PO212 with other genomes, we chose those genomes deposited in databases of Penicillium species classified as P. rubens or P. chrysogenum. Given the controversy in the current classification of P. rubens and P. chrysogenum, we first analyzed the barcode sequences corresponding to benA, caM, and RPB2 in these deposited genomes and the newly sequenced PO212 to determine which strains corresponded to each species.

Using the Clustal Omega Multiple Sequence Alignment tool (EMBL-EBI), we generated multiple alignments with each of the three barcode sequences (Supplementary Figure S1). These alignments indicated that all strains share the same sequence changes characteristic of P. rubens barcodes (indicated as pink boxes in Supplementary Figure S1). In the genomes of strains KF-25 and v1.0, benA, caM, and RPB2 barcodes displayed a mixture of sequence changes between those of P. rubens and P. chrysogenum. This analysis indicated that most probably all strains, including PO212 and with the exception of KF-25 and v1.0, may belong to the P. rubens clade, as the nucleotide sequences of the barcodes were identical to those described as P. rubens as compared to the sequence described as P. chrysogenum (Supplementary Figure S1). Supplementary Figure S2 shows the tree generated for all studied strains of Penicillium spp. by concatenating the sequences of benA, caM, and RPB2 gene fragments (between 1,879 and 1,931 bp). This tree clustered the multilocus sequence of P. rubens strain DTO 98E8 and all multilocus sequences from previously described P. rubens strains and differentiated those from P. chrysogenum strains KF-25 and v1.0. This cluster of P. rubens strains included PO212 multilocus. Consistent with the multilocus analysis shown in Supplementary Figure S2, a genome-based phylogenetic tree constructed using nucleotide sequence from assemblies of all studied strains showed the identical distribution of strains into two clusters (Figure 1).

For a more detailed comparison, multiple comparative analyses of these genomes were conducted using the Quast tool and the PO212 genome as a reference genome (Supplementary Table S1). The genomes of strains KF-25 and v1.0 were included as distant members of the P. rubens and P. chrysogenum clades. All genomes used in this comparison contained similar GC content (Supplementary Figure S3). The percentage of the genome represented relative to the PO212 genome ranged from 92.91 (KF-25) to 98.579% (HKF42). The duplication ratio was low (1.001–1.007), indicating high sequence conservation and the absence of major duplication events in the BCA PO212. In terms of mismatches, the greatest differences were found, as expected, with those that presented the greatest differences in barcode alignments (KF-25 and v1.0). However, the most similar strains to PO212 in the number of mismatches were HKF42 and Biourge 1923. The best-assembled genome, organized into five scaffolds, is P2niaD18 and the comparison with PO212 assembly showed the presence of 98.46 mismatches per 100 kbp, probably indicating the distant origins of both strains.

Quast analysis also showed the differences between the PO212 assembly and the other selected Penicillium genomes, suggesting the presence of numerous variations that could prevent immediate identification of those causing the BCA phenotype. Figure 2 provides an overview of the distinctive organization of the assemblies used in this comparison with respect to PO212. The excessive number of variations and the difficulties in testing biocontrol characteristics of the strains with a sequenced genome prompted us to choose for sequencing and detailed genomic analysis an alternative P. rubens strain that lacked the BA but was close to the PO212 strain.

The S27 strain was isolated from the soil and was primarily classified as a non-BA strain (Villarino et al., 2016). Notably, S27 displayed similar colonial morphology to PO212 (Figure 3). Sequencing of ITS1-5.8S-ITS2 regions classified S27 as a P. rubens strain, and a dendrogram based on BOX and repetitive extragenic palindromic (REP) DNA fingerprints placed S27 close to the isolate PO212 (Villarino et al., 2016).

The efficacy of S27 against Fusarium wilt in tomato plants was recently determined. In contrast to PO212, strain S27 did not significantly (p ≤ 0.05) reduce either the disease severity or incidence caused by FOL in tomato plants. FOL-inoculated and S27-treated plants showed similar symptoms to FOL-inoculated and untreated control plants (control +, Figure 3). Thus, we chose S27 for genomic sequencing as the best tool for a comparative genomic analysis with PO212 to target genes potentially related to PO212 BA.

The genome of the S27 strain was sequenced using the Illumina Hiseq 4000 platform by 150 bp paired-end sequencing reads that generated 2,345 Mbp (15,532,752 sequence reads). These reads showed coverage of ~82.3× over PO212 assembly. De novo assembly of the S27 genome was performed and evaluated in terms of N50 and L50 using Quast. Automated assembly and manual sequence verification allowed estimating for the S27 strain a genome size of 29.89 Mb. S27 genomic assembly was organized into 414 scaffolds with an N50 scaffold length of 0.42 Mb. At least 10,164 coding sequences were estimated after manual curation of ORF prediction by AUGUSTUS.² Barcode analyses of P. rubens and P. chrysogenum showed that S27 was also a P. rubens strain (Supplementary Figure S1). Table 4 shows a summary of the two genomes. The data presented in the study are deposited in the GenBank repository, accession number JAPDLD000000000.

Genome wide alignment showed the strong conservation of PO212 and S27 genomic sequences (Figure 4). MUMmer plot showed that 99.95% of the PO212 genome assembly matched that of the S27 strain with identity equal to or greater than 75%. Only 0.05% of the genomic sequence seemed to be unique to PO212 but did not contain any predicted gene model.

Quast analysis (Supplementary Table S1) also revealed the high similarity between PO212 and S27 genomes. Comparisons of the S27 and PO212 assemblies showed that up to 99.81% of the S27 genome was represented in PO212 and a duplication ratio of 1.002. In this comparison 1,890 mismatches were found, a rate of 6.34 changes per 100 kbp of genomic sequence, which was significantly lower than the previous genomic comparison (Supplementary Table S1). The phylogenetic trees using either multilocus or genomic analysis also confirmed the proximity of the two strains PO212 and S27 and other P. rubens strains by joining them in a clade (Figure 1 and Supplementary Figure S2).

Figure 5 shows a circos ideogram representing a comparison between the assembly of PO212 and S27 raw reads. S27 sequencing data present full coverage of the PO212 genome. Figure 6A shows that 97.86% of S27 reads mapped over the 65 scaffolds PO212 genomic assembly. The remaining 2.14% of S27 reads mapped to the mitochondrial genome. With the aid of Circos analysis and CLC software, we searched for those Single Nucleotide Variants (SNVs) displaying 100% allelic variation between S27 and PO212 genome. We detected at least 104 variations between both genomic sequences. Most of this low number of variations was classified as SNVs (Figures 5, 6B). In addition to 71 nucleotide changes found in coding or intergenic regions, we found seven Multi-Nucleotide Variants (MNVs), 14 insertions, nine deletions, and three replacements (Figure 6B). Only 15 of these variations were found in coding regions (Figures 5, 6C). However, manual curation still reduced the number of variations to six with amino acid change (see Figure 7 and Table 5) and two silent mutations (Supplementary Table S2).

We focused on six genes that carried an SNV causing an amino acid change. Only one SNV caused an early stop in the coding sequence. In g4471.t1 gene (orthologue of Pc16g08360 in PrWis), the CGA codon for arginine 327 was changed to a TGA stop codon. This SNV caused the truncation of the hypothetical protein S27g3229.t1 at amino acid 326, missing 3/4 of the ORF compared to the predicted PO212 protein (g4471.t1; Table 5). The remaining SNVs described in Table 5 caused punctual amino acid substitutions (Figure 7). Notably, these SNVs were found only in S27. The orthologues of PO212 showed the same nucleotide sequences as in the PrWis and P2niaD18 reference genomes.

We then confirmed the presence of these variations between S27 and PO212 genomes by PCR amplification of these regions and subsequent sequencing (Table 6). To investigate whether any of these variations were specific to the biocontrol phenotype, we chose eight strains from our collection of P. rubens strains classified accordingly to their BA. We sequenced the regions where those SNVs were mapped. We found that SNVs are present in g196.t1, g1339.t1, and g3169.t1 genes (nomenclature as in PO212 genome) specific to S27 (Table 6). For the remaining genes, it was feasible to establish two well-differentiated groups based on specific changes in g3741.t1, g3975.t1, and g4471.t1 (Table 6).

Given the absence of specific mutations or genes that could explain the BA of PO212, we performed a search for putative homologs of 13 genes involved in biocontrol in Trichoderma species (Table 7; Sharma et al., 2011). Twelve were found in the genomes of PO212 and S27 without any difference in nucleotide sequences.