Genome and transcriptome reveal lithophilic adaptation of Cladophialophora brunneola, a new rock-inhabiting fungus

ABSTRACT

Rock-inhabiting fungi (RIF) are slow-growing microorganisms that inhabit rocks and exhibit exceptional stress tolerance owing to their thick melanised cell walls. This study reports the identification of a novel rock-inhabiting fungus, Cladophialophora brunneola sp. nov. which was isolated from a karst landform in Guizhou, China, using a combination of morphological and phylogenetic analyses. The genome of C. brunneola was sequenced and assembled, with a total size of approximately 33.8 Mb, encoding 14,168 proteins and yielding an N50 length of 1.88 Mb. C. brunneola possessed a larger proportion of species-specific genes, and phylogenomic analysis positioned it in an early diverged lineage within Chaetothyriales. In comparison to non-RIF, C. brunneola displayed reduction in carbohydrate-active enzyme families (CAZymes) and secondary metabolite biosynthetic gene clusters (BGCs). Transcriptome analysis conducted under PEG-induced drought stress revealed elevated expression levels of genes associated with melanin synthesis pathways, cell wall biosynthesis, and lipid metabolism. This study contributes to our understanding of the genomic evolution and polyextremotolerance exhibited by rock-inhabiting fungi.

KEYWORDS:

• Rock-inhabiting fungi
• Cladophialophora brunneola
• genomic characteristics
• stress-responsive transcriptome
• lithophilic adaptation

1. Introduction

Fungi exhibit a global distribution, inhabiting diverse ecosystems such as tropical, polar, terrestrial, aquatic, and even endosymbiotic environments, acquiring nutrients through various ecological strategies, including saprophytic, parasitic, and predatory lifestyles (de Boer et al. 2005; Ohm et al. 2012; Yang et al. 2012, 2023; Pusztahelyi et al. 2015; Gostinčar et al. 2018; Coleine et al. 2022). Bare rock surface, characterised by limited moisture and nutrients, high UV radiation, and fluctuating temperatures, represents one of the most inhospitable habitats on earth (Gorbushina 2007). Rock-inhabiting fungi (RIF) have developed adaptations to thrive in this challenging environment, benefiting from the reduced competition with other microbes and their extremotolerance (Gorbushina et al. 1993; Wollenzien et al. 1995; Chertov et al. 2004). Despite their worldwide distribution, RIF have often been overlooked due to their small size, extremely slow growth, and lack of diagnostic morphological features. These fungi are characterised by the absence of sexual reproductive structures and form compact, melanised colonies on rocks. They are considered obligate rock dwellers with limited growth rates, often referred to as “lithophilic” fungi (Liu et al. 2021). RIF are a polyphyletic group mainly affiliated with Dothideomycetes, Eurotiomycetes, and Arthoniomycetes in Ascomycota, especially in the orders Dothideales, Capnodiales in Dothideomycetes, and Chaetothyriales in Eurotiomycetes (Gueidan et al. 2008; Ruibal et al. 2008, 2009; Egidi et al. 2014). Currently, limited genomic data is available for RIF. The first RIF genome to be sequenced was that of Cryomyces antarcticus MA 5682 (Sterflinger et al. 2014), followed by Coniosporium apollinis CBS 100218 (Teixeira et al. 2017) and Knufia petricola MA 5789 (Tesei et al. 2017). Proteomic data subsequently revealed changes in protein expression patterns under stress conditions, including a reduction in the number of expressed proteins and the upregulation of unidentified proteins (Tesei et al. 2012; Zakharova et al. 2013, 2014). It appears that RIF reduces metabolic levels as a mechanism to withstand adverse environmental factors. However, further transcriptomic and proteomic data are necessary to fully understand the mechanisms of stress adaption and metabolic activity during active growth and dormant phases.

Cladophialophora Borelli is classified as an asexual genus in the family Herpotrichiellaceae (Badali et al. 2008). It is composed of black yeast-like fungi species. The genus is considered to be simple hyphomycetes, characterised by brown hyphae that give rise to branched chains of pale brown conidia. Currently, Cladophialophora comprises over 40 species, which frequently act as opportunistic pathogens, causing a range of human infections from mild cutaneous lesions to fatal encephalitis (Horre and de Hoog 1999; Badali et al. 2008, 2009, 2010). These fungi can also be found as phytopathogens (Crous et al. 2007) and in environmental samples (Usui et al. 2016; Kiyuna et al. 2018; Das et al. 2019; Sun et al. 2020). Numerous molecular phylogenetic studies have shown that Cladophialophora is polyphyletic within the order Chaetothyriales. Additionally, this genus is closely related to various anamorphic genera, such as Exophiala, Cyphellophora, Fonsecaea, Knufia, Phialophora, and the teleomorphic genus Capronia (Kiyuna et al. 2018; Quan et al. 2020). Only five species of black yeast-like fungi have been sequenced in this genus, including three opportunistic pathogens (C. bantiana, C. carrionii, and C. yegresii) and two species from the environment (C. psammophila and C. immunda) with the ability to degrade polyaromatic hydrocarbons (Sterflinger et al. 2015; Teixeira et al. 2017).

In a survey of resources and diversity of RIF in Guizhou Province, a previously unrecognised species within the genus Cladophialophora was discovered. Through phylogenetic analysis of multiple gene sequences, including the internal transcribed spacer (ITS), small subunit of ribosomal RNA (SSU), ribosomal large subunit (LSU), translation elongation factor (TEF) exon genes, and β-tubulin (TUB) exon genes, the isolate was identified as Cladophialophora brunneola. Although RIF face numerous simultaneous stress factors in their natural habitat, the mechanisms underlying their adaptation to such niches remain poorly understood. Given that C. brunneola represents one of the typical lineages of RIF in the family Herpotrichiellaceae of Eurotiomycetes, we propose that analysing genomic and stress-response transcriptome data can provide new insights into the genome evolution and the development of polyextremotolerance in RIF.

2. Materials and methods

2.1. Sample collection on rocks and strain isolation

Rock samples were collected in a karst landform in Fanjingshan National Nature Reserve in Tongren City, Guizhou Province, China (27°55′ N, 108°42′ E, 2,220 m altitude). Vegetation was abundant around the circumstances of sampling sites (Figure 1a). Samples containing black colonies were obtained by using a sterile chisel to detach them from the stock rocks (Figure 1b). The collected samples were promptly sealed in plastic bags, stored in an ice box, and transported to the laboratory. Rock pieces (approximately 0.5–1 cm³) with black colonies were obtained from the rock samples using an industrial stone splitter (Model CM-10, Hydrasplit, Park Industries, Inc., St. Cloud, MN, USA). These rock pieces were subsequently surface-disinfected with 95% (v/v) ethanol for 3–5 seconds, transferred to physiological saline containing 0.001% (v/v) tween-20, washed with sterilised distilled water, and dried on sterilised filter papers. Approximately 1 gram of each sample was pulverised in a sterilised mortar, and the resulting powder was suspended in 2–3 mL of sterile water. Aliquots of suspensions (200 μL) were evenly spread onto Petri plates containing dichloran rose Bengal chlortetracycline (DRBC) agar with the addition of 100 mg/L streptomycin after autoclaving. The plates were then incubated at 10 °C for 4 weeks. Slow-developing colonies exhibiting dark pigmentation on the plates were transferred to MEA (malt extract 2%, peptone 2‰, and agar 1.5%) for further purification and identification. The type specimens were deposited in the Herbarium of Mycology at the Institute of Microbiology, Chinese Academy of Sciences in Beijing, China (HMAS), and ex-type living cultures were deposited in the China General Microbiological Culture Collection Center (CGMCC).

Figure 1. Sampling site and culture of Cladophialophora brunneola (CGMCC 3.18770). (a) Landscape of sampling site. (b) Rocks bearing black colonies. (c) Colony on 2% MEA after 4 weeks. (d–h) Conidial chains. (i, j) Conidia. Scale bars = 10 µm.

2.2. Morphological characterisation

The morphological characteristics of the strain were evaluated by incubating on 2% MEA for a duration of 4 weeks at 23 °C. After incubation, measurements were taken for colony diameter and colour. To examine the mycological characteristics, such as the presence of conidia, aerial mycelium, density, and pigment production, a light microscope (Nikon Eclipse 80i, Tokyo, Japan) was utilised.

2.3. DNA extraction, PCR amplification, and sequencing

The genomic DNA was extracted from fungal mycelia grown on MEA plates using the rapid “thermolysis” protocol (Zhang et al. 2010). The amplification of gene fragments, including ITS (internal transcribed spacer gene region), SSU (small subunit ribosomal RNA gene), LSU (large subunit ribosomal RNA gene), TEF (translation elongation factor), and TUB (the partial β-tubulin gene), was performed using the following primer pairs: ITS4/ITS5 for ITS (White et al. 1990), NS1/NS4 for partial SSU (White et al. 1990), LROR/LR5 for partial LSU (Rehner and Samuels 1994), 728 F/986 R for partial TEF (Carbone and Kohn 1999), and Bt2a/Bt2b for partial TUB (Glass and Donaldson 1995). The PCR programmes were performed as described by Su et al. (2015). The PCR products were subsequently sequenced by Beijing Tianyihuiyuan Bioscience and Technology after evaluation through electrophoresis.

2.4. Phylogenetic analyses

The newly obtained sequences were queried against the nuclear database of NCBI. Related available sequences from GenBank were obtained and included in the analysis (Table 1). Multiple sequence alignments of five genes were generated using MAFFT v.7.407 with Q-INS-I strategy (Katoh and Standley 2013). Conserved blocks were selected from the initial alignments using Gblocks 0.91b (Castresana 2000). For phylogenetic reconstructions, both maximum likelihood (ML) and Bayesian inference (BI) methods were employed. ML analyses were conducted using RAxML v.8.2.12 (Stamatakis 2014) with the GTR + GAMMA model. Maximum likelihood bootstrap proportions (MLBP) were determined by performing 100 replicates. Models of evolution for BI analyses were estimated using PartitionFinder 2 (Lanfear et al. 2017). Bayesian inference posterior probabilities (BIPP) were calculated using Markov Chain Monte Carlo sampling (MCMC) in MrBayes v3.2.6 (Ronquist et al. 2012) under the estimated model of evolution. Two independent analyses with four Markov chains were run for 2,000,000 generations until the average standard deviation of the split frequencies dropped below 0.01. Trees were sampled every 100th generation. The initial 25% of MCMC sampling generations were discarded as burn-in. The trees were visualised using FigTree v1.4.4.

Table 1. Characteristics of genome assembly of Cladophialophora brunneola..

Features	C. brunneola
Size (bp)	33,808,635
Coverage (fold)	171 ×
Scaffold No. (> 1 kb)	182
Scaffold No. (> 100 kb)	25
The longest scaffold length (bp)	5,633,944
The shortest scaffold length (bp)	1,948
Scaffold N50 (bp)	1,884,340
Scaffold N90 (bp)	410,579
G + C content (%)	51.03
Repeat rate (%)	4.53
Protein-coding genes	14,168
Gene density (gene per Mb)	419
Mean gene length (bp)	2,299.38
Mean cDNA length (bp)	2,183.07
Exons per gene	2.27
Mean exons length (bp)	1,348.39
Mean intron length (bp)	109.1
tRNA	185

2.5. Genome sequencing, assembly, and annotation

We sequenced the genome of C. brunneola strain CGMCC 3.18770 using the PacBio Sequel and Illumina HiSeq platform at Allwegenes Tech (Tianjin, China). Library construction and genome sequencing were carried out using the company’s standardised pipeline. The PacBio Sequel reads were obtained and subjected to trimming, correction, and assembly using the Canu assembler version 1.8 (Koren et al. 2017). To minimise error rate, the read correction parameter “correctedErrorRate = 0.035” was employed. Moreover, Illumina short reads were used to further polish the assembly genome with Pilon v1.23 (Walker et al. 2014). The completeness of genome annotation was accessed using the ascomycota_odb10 database in BUSCO v5.2.2 (Seppey et al. 2019). Knufia petricola MA 5789, Coniosporium apollinis CBS 100218 from previous research were selected for RIF group. Other fungal genomes used for comparative analysis were downloaded from NCBI-GenBank database.

Tandem Repeat Finder v.4.04 (Benson 1999), LTR_FINDER (Xu and Wang 2007), RECON v.1.08, RepeatScout v.1.0.5, and RepeatModeler v.1.0.11 (http://www.repeatmasker.org; last accessed on 13 October 2022) were used to construct species-specific repeat sequence databases through identification of structural features and de novo prediction. Subsequently, these databases were merged with Repbase (Bao et al. 2015) to form a comprehensive reference repeat database. Finally, RepeatMasker v.4.0.8 (Chen 2004) was utilised to predict types and quantities of repeat sequences in the genome.

The GETA pipeline (https://github.com/chenlianfu/geta; last accessed on 13 October 2022) was used to generate a comprehensive set of protein-coding genes. RNA-seq reads obtained from different treatment tissues were mapped to the assembly using Hisat2 (Kim et al. 2019). Ab initio gene predictions were performed using Augustus (Stanke et al. 2004), which was trained with the RNA-seq data. The predicted gene sets were functionally annotated using eggNOG-mapper (Huerta-Cepas et al. 2019) and InterProScan (Jones et al. 2014) with default parameters. Noncoding RNA was predicted using Infernal v.1.1.2 (Nawrocki and Eddy 2013) with the Rfam 14.1 database (Kalvari et al. 2018). The presence of secondary metabolite BGCs was searched on the antiSMASH website (fungal version 5.2.0) (Blin et al. 2019). For the prediction of CAZymes, we employed blastp and HMMER v.2.2 (Finn et al. 2011) against the dbCAN2 (Zhang et al. 2018) database.

2.6. Phylogenomic reconstruction

Single-copy orthologs were identified in Orthofinder v.2.5.4 (Emms and Kelly 2015). Sequence alignments were performed using MAFFT v.7.407 (Katoh and Standley 2013), and gaps were eliminated using Gblocks 0.91b (Castresana 2000). Subsequently, the genes were concatenated to generate a supergene sequence. The best-fitting model was estimated with ModelTest-NG (Darriba et al. 2020). A maximum likelihood (ML) phylogenomic tree was built using RAxML v.8.2.12 (Stamatakis 2014) with 1,000 bootstraps to infer branch support.

2.7. Culture under PEG-induced drought stress and transcriptome analysis

Polyethylene-glycol (PEG), a non-permeable osmolyte known to induce significant water deficit (Lagerwerff et al. 1961), was used to assess the drought resistance of C. brunneola. Shake culture medium was prepared by adding 0%, 5%, 10%, 20%, and 40% PEG to 2% MEB solution, and was incubated at 23 °C for a month of observation. For RNA-seq analysis, samples treated with 0% (CK) and 20% (P) PEG for 7 days of cultivation were selected. Each sample was conducted in triplicate, and RNA was extracted using Trizol (Invitrogen^TM, Carlsbad, CA, USA) (Meng and Feldman 2010). Sequencing was performed on the Illumina NovaSeq platform at Annoroad Gene Technology Co., Ltd. Library construction and genome sequencing followed the standardised pipeline of the company. TrimGalore-0.6.6 was used to remove adapters and low-quality bases, and the resulting clean reads were aligned to the genome of C. brunneola using hisat2 (Kim et al. 2015). Alignment files were then ordered using Samtools 1.6, and quantification of the mapped reads to exons was conducted with stringtie v1.3.5 (Pertea et al. 2015). The TPM (transcripts per million) method was used to calculate the expression level. Differentially expressed genes (DEGs) were identified using DEseq2 R package (Wang et al. 2010) with adjusted P-value < 0.05 and |log2Foldchange| > 1 as thresholds. Further Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs was performed using the combination of Fisher’s exact test.

2.8. Melanin coarse extraction and spectrophotometric measurements

Melanin was extracted following previously established methods with minor modifications (Pal et al. 2013; Liu et al. 2018; Pralea et al. 2019). Approximately 1 mg of fungal biomass from the culture filtrate was mixed with 1 mol/L NaOH and refluxed for 2 hours. Subsequently, the solution was centrifuged at 10,000 r/min for 5 minutes, and the resulting supernatant was acidified overnight at room temperature using concentrated HCl until pH reached to 2–3. The resulting precipitate was collected by centrifugation at 10,000 r/min for 5 minutes, and washed with 2 mL of ddH₂O and 3 mL of organic solvents (chloroform:ethyl acetate:ethanol = 1:1:1). The precipitate was then dissolved in 1 mol/L NaOH, followed by acidification using concentrated HCl again as described above. The reacquired precipitate was freeze-dried to obtain crude melanin.

Extracted crude melanin samples were suspended in 1 mL of 1 mol/L NaOH to measure absorbance in the range of 200–800, utilising the ThermoGenesys-10S UV-Vis spectrophotometer. The blank control was prepared using 1 mol/L NaOH solution.

3. Results

3.1. Taxonomy

Cladophialophora brunneola W. Sun, L. Su, M.C. Xiang and X.Z. Liu, sp. nov. (Figures 1, 2)

Figure 2. Phylogenetic tree of the Cladophialophora genus generated using maximum likelihood analysis with combined of ITS, SSU, and LSU sequences. Node values (ML/BI) indicating bootstrap support (left) and Bayesian posterior probability (right) are shown. Sequences generated in this study are highlighted in orange. The outgroup used in this analysis is Cyphellophora europaea strain CBS 129.96.

MycoBank: MB824202.

Etymology: brunneola (Lat., masc. adj.), referring to the colour of conidia and conidial chains.

Type: China, Guizhou Province, isolated from rock on 12 November 2014, by Meichun Xiang and Wenjiao Cai. The holotype is deposited at HMAS 247720 (dried culture), and the ex-type culture is preserved at CGMCC 3.18770. Another culture is also preserved at CGMCC 3.18769.

Description: Following a 4-week incubation period at 25 °C on MEA, the colonies reach a diameter of 12 mm. The central areas exhibit a light brown colouration and possess short aerial hyphae, whereas the periphery displays a dark brown, pale, and glossy appearance with pigment inside the medium, forming a zonate pattern (Figure 1c). The reverse of colonies appears olive brown. The mycelium is composed of less branched, smooth, thin-walled, septate hyphae that contain lipid droplets. These hyphae vary in colour from light to dark brown and exhibit a width of 1.1–1.9 μm. The conidia are unicellular, thick-walled, and range from hyaline to pale olivaceous in colour, from spherical to ellipsoid in shape. The dimensions of conidia are 3.5–5.5 × 4.2–4.9 μm, and they always form long chains (Figure 1d-j). No teleomorph form is observed.

Notes: Phylogenically, C. brunneola is closely related to C. humicola and C. minutissima. However, the conidia of C. brunneola are much smaller. Morphologically, C. humicola produces subcylindrical to narrowly ellipsoid conidia measuring 8–17 × 1.5–2.5 μm, characterised by an unthickened cell wall (Crous et al. 2007). C. minutissima produces cylindrical to fusiform conidia with truncated and thickened hila, measuring 8–22 × 1–2 μm (Davey and Currah 2007).

3.2. Phylogenetic analyses

The alignment results showed that the ITS sequence of C. brunneola was less than 97% similar to other Cladophialophora species, displaying 96.03% and 95.44% similarities to C. minutissima UAMH 10709T and C. humicola CBS 117536, respectively. Consequently, a 3-locus sequences alignment was analysed, incorporating ITS (1–300 characters), LSU (301–603 characters), and SSU (604–1512 characters), and involved a total of 27 taxa, with Cyphellophora europaea selected as the outgroup (Table S1). Bayesian inference (BI) and maximum likelihood (ML) generated concordant topologies. For the Bayesian analyses, the HKY + I + G model was the best fit for ITS and LSU regions, while the GTR + G model was suitable for SSU region. The maximum likelihood tree, presented in Figure 2, depicted a strongly supported grouping (MLBP/BIPP = 100/1) of CGMCC 3.18770 and CGMCC 3.18769, distinct from C. minutissima and C. humicola (Figure 2).

3.3. Genomic features of C. brunneola

The genome of C. brunneola (strain CGMCC 3.18770) was sequenced using Illumina paired-end and PacBio Sequel platforms, resulting in a total genome size of 33.8 Mb. The N50 length of assembly was 1.88 Mb. By searching approximately 100 bp in length with TTAGGG repeats, three contigs with telomeric repeats on both ends were identified. In addition, 21 scaffolds contained these characteristic telomeric repeats on either the 5′ or 3′ end (Table S2). To evaluate the completeness of the assembled genome, BUSCO analysis was performed, which detected 99.2% of the 758 expected genes, indicating a high-quality assembly.

Approximately 1.53 Mb (representing 4.53% of the genome) in C. brunneola genome were repetitive DNA sequences and transposable elements (TEs). TEs comprised 0.97 Mb (2.87%) in total, consisting of 0.47 Mb (1.39%) retroelements, 0.09 Mb (0.28%) DNA transposons, and 0.41 Mb (1.21%) unknown transposons. The most prevalent transposable elements were LTRs, occupying 0.4 Mb (1.20%) of the genome (Table S3). Using the GETA pipeline, we obtained a final dataset comprising 14,168 protein-coding gene models. The average gene length was 2,299 bp, resulting in a gene density of 419 genes per Mb. The average CDS length was 2,183.07 bp, with an average of 2.27 exons and a mean exon length of 1,348.39 bp (Table 1, Figure 3).

Figure 3. Circos graph illustrating genome characteristics of Cladophialophora brunneola. (a) 25 scaffolds (>100 kb) representing over 93.83% of the assembled genome, with large fragment duplications marked. Transposons are shown in green and retrotransposons are shown in purple. (b) GC ratio per 100 kb. (c) Secondary metabolite BGCs indicated by red arrows. (d) Gene density per 100 kb.

3.4. Phylogenomic analysis and functional gene annotation

A maximum likelihood phylogenomic tree was conducted using 2,402 single-copy orthologs from 31 fungi within the order of Chaetothyriales (Eurotiomycetes), with representative species from Dothideomycetes serving as the outgroup (Table S4). The results showed that C. brunneola was distantly related to its relatives in Herpotrichiellaceae and represented an early lineage in the Chaetothyriales (Figure 4). The large number of genes encoded in the C. brunneola genome can be attributed to the abundance of species-specific genes (2,137 genes). Similarly, the other two RIF, K. petricola and C. apollinis, which encode 9,308 and 10,469 genes respectively, exhibited a similar pattern in terms of orthologs. To determine the allocation of cellular functions within the proteome, protein functions were assigned based on sequence similarity using KOG classification. C. brunneola possessed a substantial number of genes categorised under “general function prediction only”. Additionally, in terms of the total count, C. brunneola had the highest number of proteins in the categories of “signal transduction mechanisms” (567 genes), followed by “lipid transport and metabolism” (542 genes) and “posttranslational modification, protein turnover, chaperones” (522 genes), as depicted in Figure 5. Regarding the 4,906 annotated genes with KOG function, C. brunneola had the highest proportion of genes associated with “energy production and conversion” (429 genes, 5.69%), but the smallest proportion of genes related to “coenzyme transport and metabolism” (104 genes, 1.38%) compared to other fungal genomes in this study (Table S5). Of 2,137 species-specific genes in C. brunneola, only 45 were annotated with KOG, with “general function prediction only” and “signal transduction mechanisms” having the highest number of annotated genes (Figure 5).

Figure 4. Phylogenomic relationships among representatives from Herpotrichiellaceae, trichomeriaceae, and cyphellophoraceae species. The tree was rooted using species in Dothideomycetes as outgroup. Three rock-inhabiting fungi (RIF) are highlighted in red. Genes were categorised as core (present in all species), shared (present in at least two species), and species-specific.

Figure 5. KOG functions of genes in whole genome and species-specific in Cladophialophora brunneola.

We performed a detailed analysis of the CAZymes in C. brunneola and compared them to those of other fungi, including two RIF species (K. petricola and C. apollinis) and three saprophytic fungi (Aspergillus nidulans, Hysterium pulicare, Penicillium expansum). C. brunneola exhibited a larger number of CAZymes (366 genes) (Table 2), primarily attributed to its glycoside hydrolases (GHs) (176 genes), which exceeded those of K. petricola and C. apollinis by 26 and 40 genes, respectively. In terms of proportion of total CAZymes, the profiles of 109 glycosyltransferases (GTs), 18 carbohydrate esterases (CEs), 15 carbohydrate-binding modules (CBMs), and 55 auxiliary activities (AAs) in C. brunneola resembled those of the other two RIF species. However, no protein with polysaccharide lyases (PLs) activity was identified in the genome of C. brunneola. In contrast to saprophytes, all three RIF species had significantly fewer CAZymes, particularly in GHs, PLs, CBMs, and AAs. Nonetheless, GTs constituted a higher proportion of these RIF species compared to saprophytes.

Table 2. Comparison of CAZymes between RIF and saprophytes.

	RIF			Saprophytes
	Coniosporium apollinis	Cladophialophora brunneola	Knufia petricola	Aspergillus nidulans	Hysterium pulicare	Penicillium expansum
Total	286	366	297	512	552	496
GHs	136 (47.5%)	176 (48.1%)	150 (50.5%)	256 (50%)	257 (46.5%)	247 (49.8%)
GTs	77 (26.9%)	109 (29.7%)	88 (29.6%)	90 (17.5%)	105 (19%)	105 (21.2%)
PLs	2 (0.7%)	0 (0%)	2 (0.6%)	23 (4.5%)	4 (0.7%)	9 (1.8%)
CEs	14 (4.9%)	18 (4.9%)	15 (5.1%)	33 (6.4%)	36 (6.5%)	26 (5.2%)
CBMs	11 (3.8%)	15 (4.1%)	14 (4.7%)	43 (8.4%)	35 (6.3%)	40 (8.6%)
AAs	50 (17.4%)	55 (15%)	36 (12.1%)	90 (17.6%)	137 (24.8%)	95 (19.1%)

A total of 26 putative BGCs responsible for secondary metabolite were identified in C. brunneola, encompassing 12 nonribosomal peptide synthetases (NRPSs), 6 polyketide synthases (PKSs), 3 terpene synthases, 2 beta lactone synthetases, 2 ribosomal synthesised and post-translational modified peptides (fungal RiPPs), and 1 NRPS-PKS hybrid enzyme (Table 3). In comparison to saprophytic species, the secondary metabolite biosynthetic gene clusters in RIF species exhibited a significant reduction, mainly due to the absence of PKS and NRPS clusters.

Table 3. Comparison of secondary metabolite BGCs between RIF and saprophytes.

	RIF			Saprophytes
	Coniosporium apollinis	Cladophialophora brunneola	Knufia petricola	Aspergillus nidulans	Hysterium pulicare	Penicillium expansum
Total	15	26	11	53	26	65
NRPS	4	12	6	17	8	24
PKS	3	6	2	15	13	18
Beta lactone synthetase	2	2	0	3	3	1
Terpene synthase	5	3	3	9	0	5
Fungal-RiPP	0	2	0	0	0	0
Hybrid	1	1	0	7	2	15
Indole synthase	0	0	0	2	0	2

3.5. Cultured features and transcriptome analysis of C. brunneola under PEG-induced drought stress

The accumulation of melanin in the culture filtrate and fungal biomass, as observed under microscopic examination (Figure 6a), showed a positive correlation with the increase in PEG concentration. Moreover, crude melanin extracts were obtained from the mycelia using an alkaline dissolution and acid precipitation procedure, with the exception of the 40% PEG culture due to its stickiness. UV absorption analysis of crude melanin extracts within the range of 200–800 nm further corroborated the increase in melanin accumulation with higher PEG concentrations, as depicted in Figure 6b, 6c.

Figure 6. (a) Melanin accumulation and hyphal morphology of Cladophialophora brunneola in liquid culture, scale bars = 10 µm. (b) UV-Vis absorption curve of crude melanin extraction (200–800 nm). (c) Absorption value at 220 nm of crude melanin extraction from samples treated with varying PEG concentrations.

Simultaneously, we conducted RNA-seq analysis on mycelia harvested from liquid cultures supplemented with 0% and 20% PEG over a period of 7 days. Alignment of clean reads to the C. brunneola genome revealed expression in 80.24% (11,368 genes) of the predicted genes across at least one sample (Table S6). Prominent distinctions were observed between PEG-treated and control groups (Figure 7a), with the identification of 1,759 differentially expressed genes (DEGs). Among these DEGs, 1,008 genes were upregulated, while 751 genes showed downregulation in the PEG-treated samples compared to control samples without PEG-induced stress (Figure 7b).

Figure 7. (a) Pairwise Pearson correlation coefficients comparing the transcriptome data of control (CK) and PEG-treated (P) samples in Cladophialophora brunneola. (b) Volcano plots illustrating differentially expressed genes (DEGs) between control and PEG-treated samples.

In the genome of C. brunneola, we identified homologs participating in three different melanin biosynthesis pathways (Table S7) as well as upregulated genes that may contribute to melanin production (Figure 8). Under PEG-induced stress, 8 out of the 14 genes in the DHN-melanin pathway were up-regulated in three replicates (Figure 8a). In the DOPA-melanin pathway, tyrosinases and laccases, which play vital roles, were activated and expressed (Figure 8b). CLBR09244, a homolog of tyrosinase, exhibited more than a 4-fold increase in expression compared to the control condition (TPM 783.8 vs. TPM 193.7), significantly surpassing the average expression level of all genes (TPM 71.2). Additionally, C. brunneola conserves the alternative pathway of producing a third type of melanin through the utilisation of L-tyrosine degradation pathway, which results in the production of pyomelanin as a side-product (Schmaler-Ripcke et al. 2009; Keller et al. 2011) (Figure 8c). In response to drought stress, tyrosine aminotransferase (Tat) and hydroxyphenylpyruvate dioxygenase (hppD) exhibited slight upregulation, while hmgA, maiA, and fahA were downregulated, leading to pyomelanin accumulation.

Figure 8. Gene expression levels in melanin synthesis pathways of Cladophialophora brunneola. (a) DHN-melanin pathway. (b) DOPA-melanin pathway. (c) L-tyrosine degradation pathway.

Under PEG-induced stress, the upregulated genes (z-score > 0) showed enrichment in GO terms related to organisation or biogenesis of the cell wall (GO:0071554; P-value = 5.4e-10), metabolic process involving macromolecules in the cell wall (GO:0044036; P-value = 1.1e-08), metabolic process involving polysaccharides in the cell wall (GO:0010383; P-value = 2.34e-08), and metabolic process involving chitin in the cell wall (GO:0006037; P-value = 3.20e-07) in biology process (BP) category. Additionally, in the cellular component (CC) category, the enriched GO terms include the cell periphery (GO:0071944; P-value = 3.41e-05), the cell septum (GO:0030428; P-value = 2.31E–04), cell tip (GO:0051286, P-value = 2.81E–04) and the cell wall (GO:0005618, P-value = 9.43E–04) (see Figure 9a and Table S8).

Figure 9. (a) GO enrichment plot of different expression genes (DEGs) in Cladophialophora brunneola, with z-score calculated as (number of upregulated genes – number of downregulated genes) divided by the square root of the number of genes in a specific GO term. (b) KEGG pathway enrichment of DEGs. (c) Expression level of genes involved in fatty acid degradation and biosynthesis.

The results of KEGG pathway enrichment analysis indicated the upregulation of genes associated with lipid metabolism (Table S9). These included fatty acid degradation (ko00071, P-value = 1.85E–08), biosynthesis of unsaturated fatty acids (ko01040, P-value = 1.46E–03), metabolism of alpha-linolenic acid (ko00592, P-value = 4.55E–03), and degradation of valine, leucine and isoleucine (ko00280, P-value = 1.18E–02) and transport and catabolism in peroxisome (ko04146, P-value = 2.55E–04) (Figure 9b). Additionally, the upregulated genes were enriched in amino sugar and nucleotide sugar metabolism (ko00520, P-value = 1.36E–03), pentose and glucuronate interconversions (ko00040, P-value = 3.63E–03) related to carbohydrate metabolism, as well as tyrosine metabolism (ko00350, P-value = 1.44E–02) and arginine and proline metabolism (ko00330, P-value = 3.32E–02) related to amino acid metabolism.

Twenty-six genes involved in mitochondrial and peroxisomal fatty acid degradation were upregulated in DEGs of lipid metabolism (Table S10, Figure 9c), likely providing energy for cellular maintenance. Two acyl-CoA synthetases (ACSL and ACSS3) among these DEGs play a critical role in catalysing the activation of free fatty acids (FAs) by forming acyl-CoA thioesters. Carnitine O-acetyltransferase (CPT, CLBR05220) and mitochondrial carnitine/acylcarnitine transporter (CACT, CLBR00227, and CLBR01704), which are components of the carnitine system responsible for transferring long chains of fatty acyl-CoA across the mitochondrial inner membrane, were also upregulated. Five acyl-CoA oxidases or dehydrogenases (ACOX, ACADM, and ACADSB), six enoyl-CoA hydratases (ECI1, ECI2, ECH1, and ECHS1), thiolases (ACAA1, ACAT, and SCPX), and thioesterases (PTE, CLBR11962) were also upregulated in the subsequent cycle of β-oxidation. Alcohol dehydrogenase (CLBR07439), the first and rate-limiting enzyme in the oxidation of fatty acids in peroxisomes, was also slightly upregulated. Additionally, the production of unsaturated fatty acids by C. brunneola as a means to maintain membrane fluidity was indicated by the upregulation of 3-oxoacyl-[acyl-carrier-protein] reductase (FabG, CLBR07084) and fatty acid desaturase (FAD2, CLBR04329).

Compatible osmolytes, such as trehalose, polyols, proline, and malondialdehyde are essential for cells to resist osmotic stress (Jennings 1984; Jennings and Burke 1990; Pfyffer et al. 1990; Hounsa et al. 1998). CLBR10229 and CLBR08760, which encode pyrroline-5-carboxylate reductase and trehalose 6-phosphate synthase, respectively, were significantly upregulated after PEG treatment (Table S6), resulting in the accumulation of pyrroline and trehalose.

4. Discussion

CGMCC 3.18770 and CGMCC 3.18769 are two isolates that represent the novel fungus, C. brunneola sp. nov., from the genus Cladophialophora. Phylogenetic comparisons suggest that C. minutissima, C. humicola, C. sylvestris, and C. modesta are the closest relatives of C. brunneola. However, extensive phylogenetic analyses demonstrate that C. brunneola is a distinct species, separate from these related fungi in the family Herpotrichiellaceae, and it appears in the early lineage of the order Chaetothyriales. These results suggest that C. brunneola exhibits significant differences in single-copy ortholog sequence compared to its relatives and has undergone specific environmental selection during genome evolution.

The influence of natural selection in the evolution of RIF genomes may be much deeper than phylogeny. Comparative genomic analysis of C. brunneola and other two published RIF species reveals genomic characteristics consistent with a rock-inhabiting lifestyle. Despite having a relatively small genome size, RIF displayed higher gene density (374.4 vs. 339.9 genes per Mb) compared to saprophytic fungi. Genes encoding enzymes involved in the degradation of glucans, polysaccharides, and chitin were reduced. For example, chitinase (GH18), chitosanase (GH75), CBM50, and CBM18, which bind to chitinous components, as well as gluco/chitooligosaccharide oxidase (AA7) and family GT31, which encode a diverse set of transferases, were less abundant in RIF. Fungal secondary metabolites (SMs) are known to play a role in self-protection, defence, competition, and communication with other organisms (Keller 2015). The average number of SMs biosynthetic gene clusters in RIF was almost one-third of saprophytes fungi (17.3 vs. 48). Clearly, RIF have reduced primary and secondary metabolic requirements, enabling themselves to be compatible with oligotrophic rock surfaces with fewer competitors.

The rock surface and its surroundings are inherently stressful environments, characterised by extreme temperature, UV exposure, and osmotic pressure. In this study, drought and osmotic stress were simulated using PEG treatment. We observed that C. brunneola produces significant quantities of melanin, either excreted into surroundings or deposited on the cell wall. Melanin, a pigment widely distributed in both Prokaryotes and Eukaryotes, serves as an ultraviolet absorber (Plonka and Grabacka 2006), provides structural rigidity in cell walls, and can store water and ions, making it advantageous in preventing desiccation (Bell and Wheeler 1986; Butler and Day 1998). Furthermore, genes involved in cell wall biogenesis and melanin synthesis pathways were found to be upregulated under PEG-induced stress, suggesting that the thickening cell wall creates more space for the increased deposition of melanin. The phenomenon of upregulated genes significantly enriching lipid metabolism is intriguing. It is speculated that the abundance of lipid droplets possessed in the cells of RIF (Gorbushina 2003; Onofri et al. 2007) is critical for desiccation tolerance, as the oxidation of fatty acids can supply cells with energy and water.

In conclusion, this study identified a novel rock-inhabiting fungus, Cladophialophora brunneola sp. nov., which was isolated from a karst landform in Guizhou, China, using a combination of morphological and phylogenetic analyses. The genome characteristics and transcriptional response of C. brunneola CGMCC 3.18770 clearly illustrated its adaptation to diverse stresses and harsh environments. The integrity and rigidity of the cell wall play a crucial role in the survival of C. brunneola in its niche. The constitutive production and upregulation of melanin biosynthesis throughout its life cycle, during drought conditions, enable melanised C. brunneola to survive and adapt to a multi-stress environment. Furthermore, C. brunneola shares genomic characteristics in common with rock-inhabiting fungi from various taxa, indicating a convergence in genome evolution. Therefore, the availability of genome sequences for other rock-inhabiting fungal species will facilitate more comprehensive investigations of the convergent mechanisms underlying their adaptation to extreme environments.

Acknowledgments

The author thank Dr. Baosong Chen (State Key Laboratory of Mycology, Institute of Microbiology, Chinese Academy of Sciences, Beijing, China) for kindly help in melanin extraction.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The sequencing data are available at NCBI BioProject ID PRJNA952342.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/21501203.2023.2256764

Additional information

Funding

The work was supported by the National Natural Science Foundation of China (31970011, 31670014).