Whole-genome sequencing was completed by BGI Genomics Co., Ltd. The genomic DNA was extracted from single spore mycelia of the Penicillium janthinellum P1. Extracted genomic DNA (Grijseels et al., 2016) then was sequenced using high-throughput Illumina Hiseq 2500 sequencing platform for mate paired-end sequencing together with the PacBio RSII single-molecule real-time sequencing platform, respectively.

The Illumina raw reads were evaluated using FastQC v0.11.6¹, and processed by utilizing the NGS QC Toolkit v2.3.3 (Patel and Jain, 2012) with Phred 20 as the cutoff to delete the low-quality reads. The resulting reads with linker and adaptors, with >10% unknown bases (N) were filtered out for further analysis. FastUniq v1.1 (Xu et al., 2012) was employed to remove PCR repeats and then the high quality paired-end reads were obtained for following assembly. An in-house bash5tools.py command of SMRT Analysis was implemented to filter low-quality reads and obtain clean reads of PacBio sequencing.

The clean reads of the PacBio sequencing were de novo assembled to obtain the original assembly contigs by FALCON, then were optimized and upgraded by FinisherSC (Lam et al., 2015) using the clean reads of the PacBio sequencing. Bowtie2 was conducted to map the clean reads of the Illumina sequence back to the PacBio assembly contigs to obtain the sam file. Use Samtools to convert sam files to bam files, and sort and index them. The upgraded assembly contigs and indexed bam file were used as the input file of Pilon to obtain the final assembly. The quality of the assembly was assessed by QUAST (Niegowski and Eshaghi, 2007).

Protein-coding genes in P1 genome were predicted using four softwares including GeneWise² (Madeira et al., 2019), Augustus³ (Keller et al., 2011), GeneMark-ES⁴ (Borodovsky and Lomsadze, 2011), and MAKER⁵ (Cantarel et al., 2008). In this paper, Augustus was trained on gene models for P1 by combining the Coprinus gene model and RNA-Seq data of P1. GeneWise and GeneMark-ES with default parameters achieved different gene prediction sets, respectively. Finally, all prediction gene results were integrated by MAKER.

The putative proteins of P1 were aligned against the RefSeq non-redundant proteins (NCBI Nr) and the SWISS-PROT Protein Knowledgebase (Swiss-Prot) using BLAST and the HMMER v3.3 (Klingenberg and Marugan-Lobon, 2013) search against the Pfam protein families database (Pfam⁶) (El-Gebali et al., 2019; Mistry et al., 2021) to assign general protein function profiles. The predicted genes were then aligned against the Conserved Domain Database (CDD⁷) using RPS-BLAST. Metabolic pathways of the P1 genome were classified based on the KEGG Ortholog results against the Kyoto Encyclopedia of Genes and Genomes (KEGG⁸). KOG enrichment analysis was assigned based on the best match of alignments against The Clusters of Orthologous Groups of proteins (COG⁹) through BLASTP. RepeatMasker v4.0.7 and RepeatModeler v1.0.5 software were used to filtrate and annotate the retroelements, DNA transposons, simple repeats and low-complexity DNA sequences. Simple sequence repeats (SSRs) were annotated by MISA. The tRNA and anticodons of the P1 genome were predicted by tRNAscan-SE v1.23 (Schattner et al., 2005). rRNA genes were predicted by Barrnap v0.9. Other non-coding RNAs such as sRNA and snRNA were identified using Rfam v9.1. AntiSMASH was used to predict the putative secondary metabolism gene clusters by using genome files and annotation files.

24 fungal species related to Ascomycota together with P1 were used for phylogenetic analysis (Supplementary Table 7) (Wang et al., 2019). The protein sequences of 25 fungi were compared by BLASTP (e-value < le-5). Then, OrthoMCL under default parameters was used to analyze the blastp results. Batch multiple sequence alignment of these genes were calculated by MAFFT v7.221 software (Katoh and Standley, 2013), subsequently combined into long sequences for each species. The conserved sites of multiple sequence alignment results were extracted using Gblocks v0.91b (parameters -t = p), the best amino acid replacement model was determined by ProTest v3.2. The phylogenetic tree was constructed using RAxML-8.1.24 (bootstraps 1000) by the input of the final alignment sequence. Three pairs of species’ divergence times were estimated by the web¹⁰ and added three fossil calibration points to the molecular clock analysis (Floudas et al., 2012). The most recent common ancestor (MRCA) of Aspergillus oryzae and Aspergillus terreus were diverged at 65 MYA. Gloeophyllum trabeum and Saccharomyces cerevisiae were diverged at 628 MYA, while Trichoderma parareesei and Aspergillus clavatus were diverged at 343 MYA (Chen et al., 2016). The divergence time of other nodes was calculated by r8s v1.80 (PL method, TN algorithm and the smoothing parameter value set to through cross-validation). The tree was idealized in iTOL¹¹. The orthologous gene family expansion was analyzed by CAFE v3.1 based on the ultrametric tree (Han et al., 2013). The OCGs undergone expansion/contraction in P1 were predicted by Swiss-Prot. Genome annotation files and protein sequences of Penicillium janthinellum P1, Penicillium oxalicum 114-2 and Penicillium janthinellum ATCC¹² were used to perform collinearity analysis via MCScanX software (Wang et al., 2012). Genome comparison was investigated by Mauve v2.3.1¹³.

CAZymes were annotated by dbCAN (covered fraction ratio > 0.3, e-value < 1e-5) and BLASTP against CAZyme database (e-value < 1e-5, maximum hit number set 500). The blastp results are integrated with the HMM results to achieve the final CAZyme annotation. Genes related to heavy metal chromium resistance, which are normally regulated by heavy metal stress, were found and analyzed from the proteins encoded by the P1 genome.

Through swiss-prot annotation of the whole protein group and gene annotation file, the locations of target encoding genes on the scaffold were determined. Analyses of other key genes adjacent to the target genes were conducted to try to interpret the characteristics of the strain. Meanwhile, the target gene’s locus of other strains with similar genes was analyzed to observe the diversity of genes.

As shown in Table 1, a total of 34,520,396 reads were produced using Illumina sequencing technology, thus yielding 5,178 Mb of raw sequencing data. The results of PacBio Sequel platform were shown in Table 2. The total number of bases measured is about 13 GB.

The pre-processed Illumina and Pacbio clean data were assembled using SOAPdenovo and FALCON, respectively. The resulting genome size was approximately 31 Mb, resulting in 1,289 contigs with Contig N50 of 50.4 kbp using SOAPdenovo. FALCON pipeline was used for the assembly of PacBio sequencing data, generating 143 contigs with Contig N50 of 3077.5 kbp (Table 3).

With third-generation data for assembling and second-generation data for evaluation and correction, the de novo genome assembly of P. janthinellum P1 with a total length of 30.8 Mb (Table 3) was yielded. A total of 29 contigs were generated with Contig N50 of 3,088.6 kbp (Table 3). This size was consistent with the estimated genome size of 29–34 Mb for three Penicillium species (Yang et al., 2016; Peng et al., 2020; Zhang et al., 2020).

10,955 CDS were identified in strain P1 and the average gene length was 1591 bp using GeneWise, Augustus, GeneMark-ES, and MAKER software for gene prediction.

Through repeated sequence analysis, there were about 1.90% of the repeated sequences (Supplementary Table 3). Among DNA transposons, Tc1-IS630-Pogo accounts for the largest proportion, amounting to 0.36% of the genome. Among the retrotransposons, LTR transposons are the majority, accounting for 0.18% of the genome.

Although simple sequence repeats are usually more likely to occur in non-coding regions, tandem repeats in the coding region of fungal genes cannot be ignored. Research by relevant scholars has shown that SSR occurs in the coding region, which will cause repeated patterns in protein sequences, thus quickly increasing the speed of protein evolution and make fungi adapt to environmental changes faster (Clayton et al., 2017). Therefore, MISA was employed to annotate simple sequence repeats of the coding region of the P1 genome, and the SRR annotation results were shown in Supplementary Table 4.

Through analyzing the distribution of tRNA genes, codon usage patterns and amino acid composition at the genomic level is an effective method to decipher the biologically important characteristics of intracellular tRNA and the bias of codon and amino acid usage. A total of 274 tRNAs were annotated, of which there were no false tRNAs (Supplementary Table 5). The number of tRNAs recognizing the same amino acid varied. There were six tRNAs over 20, which were tRNA-Val (27) carrying valine, tRNA-Arg (26) carrying arginine, and glutamic acid. TRNA-Glu (25), tRNA-Ile carrying isoleucine (24), tRNA-Phe carrying phenylalanine (22), tRNA-His carrying glutamic acid (21). The frequency of use of tRNAs carrying anticodons was also more scattered, with 39 groups of tRNAs that recognize the same amino acid carrying the same anticodons.

rRNA is the structural and functional core of the ribosome, which has peptidyl transferase activity. rRNA mainly serves as a binding site for tRNA and a variety of protein synthesis factors, and often binds to mRNA in the extension of peptide chains (Supplementary Table 6).

Total 19150 Ortholog Cluster Groups (OCGs) were constructed, 8764 of which contained 9776 P1 proteins through OrthoMCL analysis. Approximately 20.2% of predicted proteins in P1 had orthotics in all species, whereas 10.8% of proteins were not orthologous to other fungi, and 19.9% of proteins had at least ten paralogs (Figure 2 and Supplementary Table 8), indicating that the P1 strain shares high genetic similarity with the selected Aspergillus and Penicillium strains.

The phylogenetic tree was constructed based on 1423 selected single-copy orthologous genes presented in ≥70% of species, revealing the evolutionary history of the P1 strain. The topological structure of the tree was consistent with the taxonomic classification of species (Figure 3). The molecular clock analysis showed that Penicillium oxalicum 114-2 had the closest divergence time, which was estimated to be 74 million years ago (MYA). The long divergence time with other strain may be the isolation effect brought by the living environment (Nayfa et al., 2020).

A total of 1183 expansion/contractions family numbers were identified and 102 OCGs were found to have expansion/contraction in P1. After annotation and analysis of those OCGs, it was confirmed that the largest expanded OCG, which encodes Malformin synthetase mlfA, harbors seven genes in P1. The expanded mlfA may indicate that P1 has the potential for biosynthesis of secondary metabolites. Other 17 OCGs including Dynein assembly factor, Aconitate hydratase, Autoinducer 2 sensor kinase/phosphatase, iron transporter, Tol-Pal system protein, oxygen-dependent choline dehydrogenase, high-affinity glucose transporter, L-threonine 3-dehydrogenase, serine/threonine-protein kinase were also expanded in P1. Some enzymes related to respiration have also expanded, which may facilitate strains to resist heavy metal toxicity. Additionally, 306 and 431 syntenic blocks were explored on the basis of the conserved gene order between P1 and Penicillium oxalicum 114-2, P1 and Penicillium janthinellum ATCC, respectively (Supplementary Figure 1). Average protein similarity between the proteome predicted by P1 with Penicillium oxalicum 114-2 and Penicillium janthinellum ATCC was evident (Supplementary Figure 2).

A wealth of varied secondary metabolites were found in P1 using the AntiSMASH. A total of 33 secondary metabolism gene clusters were identified. The predicted putative products of 24 secondary metabolism gene clusters were classified into eight different types: seven T1 polyketide synthase (T1pks) gene clusters, seven non-ribosomal peptide synthase (Nrps) gene clusters, two terpene gene clusters, three T1pks-Nrps gene clusters, two fatty acid gene cluster, one Indole-Nrps gene cluster, one saccharide gene cluster, and one Lantipeptide gene cluster; the remaining nine gene clusters synthesized other unknown secondary metabolites (Supplementary Table 12).

Some microorganisms can secrete extracellular polymeric substances (EPS) such as glycoproteins and soluble amino acids which have the effect of complexing or precipitating metal ions (Naveed et al., 2020). The release of EPS by microorganisms needs the aids of secretory systems and transporters, therefore signal recognition particle (SRP) plays a critical role in the targeting of secretory proteins to cellular membranes (Gowda et al., 1998), and Sec secretion system is responsible for the transport of proteins on the plasma membrane (Fagan and Fairweather, 2011). Within P1 genomic information, 22 proteins were predicted to be components of the eukaryotic Sec-SRP secretion systems (Supplementary Table 9).

Carbohydrate and energy metabolism acted as key attribution in response to heavy metal adaptation (Kan et al., 2019). Comparing carbohydrate active enzymes of similar fungi may reveal the anti-heavy metal properties of P1. The CAZy database¹⁴ contains a series of enzymatic catalytic modules and sugar binding modules that degrade, modify or create glycosidic bonds. A total of 525 encoding CAZymes were annotated using dbCAN. The putative CAZymes included 92 Auxiliary Activities (AAs), 44 Carbohydrate-Binding Modules (CBMs), 27 Carbohydrate Esterases (CEs), 250 Glycoside Hydrolases (GHs), 104 GlycosylTransferases (GTs), 8 Polysaccharide Lyases (PLs), of which carbohydrate active enzymes such as glycoside hydrolase were most annotated (Supplementary Table 10). Compared with P. janthinellum ATCC 10455 and P. oxalicum 114-2, families of AAs in P1 expanded while families of CBMs contracted. AAs acts as redox enzymes that act in conjunction with CAZymes (Supplementary Table 11). This reveals a high oxidoreductase-coding activity in P1, which may be more conducive to the reduction of toxic heavy metals. Cytochrome-c peroxidase of AA2 (EC 1.11.1.5) and p-benzoquinone reductase (NADPH) of AA6 (EC 1.6.5.6) were identified within the P1 genome which were anti-chromium associated proteins. CBM1 as an important cellulose-binding domain family in cellulose degradation was rare in the P1 fungal genome, inferring that it may be related to marine living environment.

Further analyses revealed that 190 CAZymes were predicted to be sugar metabolism enzymes, involving 68 polysaccharide degrading enzymes, 100 glucose degrading enzymes and 12 polysaccharide synthases. The P1 genome was enriched in glucose degrading enzymes assigned to 10 predicted CAZyme families. This finding was compared to 99 genes in P. janthinellum ATCC 10455 and 54 genes in P. oxalicum 114-2. Compared with the above two strains, larger numbers of genes encoding the important glucose oxidase were classified into AA3 in the P1 genome. This result supports the high sugar metabolism activity of P1, and provides a certain insight into anti-metal activity.

The maximum Cr (VI) biosorption capacity for living fungus P1 pellets was about 87% at the optimum condition (Chen et al., 2019), demonstrating strong chromium removal ability (Sriharsha et al., 2020). By analyzing the genes annotated in the P1 genome, the genes that may be related to their high resistance to heavy metal chromium are listed below (Table 4). As we can see, cytochrome P450s (CYPs) were heavily annotated in P1. As a multicomponent electron transport chain system, CYPs were critical in degradation, detoxification, and syntheses of life-critical compounds in organisms (Deng et al., 2007). In addition, CYPs were also found in Phanerochaete chrysosporium and the unique P450s had a potential for the production of useful flavonoids (Kasai et al., 2010). Superoxide dismutase can eliminate radicals (Liu et al., 2020), which may be involved in the process of reducing Cr (VI) by fungi.

The encoding genes flanking to the ChrA and ChrR gene loci were studied in detail, seen in Figure 4. ChrA gene located on scaffold 14 can encode chromate transport protein. OpS5 acted as an oxidoreductase is attributed to secondary metabolite biosynthesis. Its physical distance from ChrA gene is about 5.6 Kb. ChrA gene is also present in other similar fungi, and a number of oxidoreductases are present around it. ChrAs in P. aeruginosa (Cervantes et al., 1990), A. eutrophus (Nies et al., 1990), Shewanella sp. (Aguilar-Barajas et al., 2008) have been studied in detail. These studies examined the resistance of direct homologous microorganisms carrying the ChrA gene to different Cr (VI) concentrations, and found that the ChrA gene alone could not explain the different Cr (VI) resistance of different microorganisms.

The ChrR gene is located on scaffold 16 and can encode quinone reductase. Adjacent to the ChrR, the Ornithine aminotransferase RocD and ABC Transporter C family (ABCC) were identified. Interestingly, the same genetic arrangement was found in Penicillium oxalicum 114-2, Penicillium expansum ATCC 24692, and Penicillium chrysogenum Wisconsin 54-1255, and these three genes were arranged in the same order. Part of Penicillium species are resistant to chromium (Fukuda et al., 2008; Zhang et al., 2011; Long et al., 2020), which may be explained by the special gene arrangement to a certain extent. In addition, Ornithine aminotransferase catalyzes the conversion of Ornithine into glutamic semialdehyde which subsequently can catalyze the reduction of NADP⁺ to NADPH. The generated NADP⁺ and NADPH will then serve as cofactors of ChrR to catalyze the reduction of Cr (VI). Quinone reductase acted as Cr (VI) reductase was confirmed to be an NADH-dependent reductase and also found in L. fusiformis ZC1 (He et al., 2011). ABCC mediates the transport of glutathione (GSH), also be known as an important antioxidant, which can eliminate the ROS produced by the ChrR enzyme during the reduction of Cr ions. GSH is crucial for chromium resistance to microorganisms. Glutathione synthase and glutathione s-transferase were up-regulated under chromium pressure in P. aeruginosa (Kılıç et al., 2010). Also the increase of Cr (VI) toxicity in E. coli was reported due to the absence of GSH (Helbig et al., 2008). ChrR gene, RocD and ABCC also exist in Trichoderma parareesei CBS 125925 and Saccharomyces cerevisiae YB210, but they are not on the same scaffold and far apart. We hypothesis this potentially may be a reason for differences in chromium resistance. ABCC alone was not related to the presence of heavy metal-related genes. In summary, it is noteworthy that whether the target resistance gene has such a genetic arrangement may be used as a basis for judging whether the strain is resistant to heavy-metal chromium at the genetic level.