Genome, genetic evolution, and environmental adaptation mechanisms of Schizophyllum commune in deep subseafloor coal-bearing sediments.

To understand the genomic evolution and adaptation strategies of fungi to subseafloor sedimentary environments, we de novo assembled the genome of Schizophyllum commune strain 20R-7-F01 isolated from ∼2.0 km-deep, ∼20-millionyearsago (Mya) coal-bearing sediments. Phylogenomics study revealed a differentiation time of 28-73 Mya between this strain and the terrestrial type-strain H4-8, in line with sediment age records. Comparative genome analyses showed that FunK1 protein kinase, NmrA family, and transposons in this strain are significantly expanded, possibly linking to the environmental adaptation and persistence in sediment for over millions of years. Re-sequencing study of 14 S. commune strains sampled from different habitats revealed that subseafloor strains have much lower nucleotide diversity, substitution rate, and homologous recombination rate than other strains, reflecting that the growth and/or reproduction of subseafloor strains are extremely slow. Our data provide new insights into the adaptation and long-term survival of the fungi in the subseafloor sedimentary biosphere.

Introduction

Subseafloor sediment is one of the largest and taxonomically diverse microbial habitats and organic carbon pools, with 2.9×1029 microbial cells, accounting for 0.18–3.6% of the Earth’s total living biomass (Hoshino et al., 2020; Kallmeyer et al., 2012; Parkes et al., 2014; Wormer et al., 2019). Besides the relatively well-studied prokaryotes, a growing number of studies have shown that fungi are also an important component of the subseafloor biosphere, with the deepest distribution down to ∼2.5 km (Liu et al., 2017; Pachiadaki et al., 2016). Most of these fungi are affiliated with Ascomycota, Basidiomycota, and Chytridiomycota, whose numbers and distribution are not closely related to depth, but are significantly related to the total organic carbon, electron accepters, temperature, and salinity (Orsi et al., 2013; Rédou et al., 2015). Phylogenetic analysis and the fungal-host co-deposition hypothesis suggest that subseafloor fungi are closely related to terrestrial species, and almost all subseafloor fungi exist in terrestrial or aquatic habitats (Liu et al., 2017; Orsi et al., 2015; Quemener et al., 2020). Previous studies have shown that subseafloor fungi may be alive and metabolically active under the extremely energy-limited conditions of the deep subseafloor (Orsi et al., 2013,2015; Sohlberg et al., 2015). Our previous laboratory investigation also revealed that the subseafloor S. commune possess normal growth ability under anaerobic conditions and employ various strategies, such as ethanol fermentation and production of amino acids, increase in mitochondrial number, and autophagosome formation, to adapt to an anaerobic environment (Zain Ul Arifeen et al., 2021a, 2021b). However, the fungus could not produce fruit bodies in the absence of oxygen, suggesting that fungi can only reproduce asexually in deep subseafloor sediments (Zain Ul Arifeen et al., 2021b). All these data indicate that subseafloor fungal species may have unique metabolic activities that helped them sustain life in the deep subseafloor environments for up to ∼81 to 48 Mya (Ivarsson, 2012; Ivarsson et al., 2013). Although a large number of studies highlight the diversity and activity of fungi in the deep subseafloor biosphere (Liu et al., 2017; Orsi et al., 2015; Pachiadaki et al., 2016; Quemener et al., 2020; Rédou et al., 2015), how fungi evolved and adapted to the energetically challenging subseafloor environment and have persisted over millions of years remains largely unknown.

S. commune is one of the most widely distributed mushrooms in the forests on Earth and an effective wood-degrading basidiomycete, which can produce a large number of hydrolases (xylanase, pectinase, cellulase, and endoglucanase) and oxidoreductase like laccases (Tovar-Herrera et al., 2018). It has also served as a model mushroom to study cell wall biogenesis, hyphal fusion and development, matingtype, heterologous expression of genes, and gene deletions (Ohm et al., 2010). Sequence search of NCBI and JGI databases and previous reports revealed that 58 strains of S. commune had been completely or partially sequenced, of which type strain H4-8 had been fully studied (Ohm et al., 2009). Based on culture study of fungal diversity in subseafloor lignite coal-bearing sediment samples collected during the Integrated Ocean Drilling Program (IODP) Expedition 337 (Inagaki et al., 2015), the distribution of S. commune has been extended to the sediment down to 2,457 meters below seafloor (mbsf) (Liu et al., 2017). S. commune, as one of the predominant subseafloor fungi in those particular geologic horizons, not only has special characteristics of anaerobic growth, but also has a unique anaerobic energy generation pathway, which might endow it to subsist energy-limited subseafloor sedimentary environments (Zain Ul Arifeen et al., 2021a). This long-term physical isolation of S. commune in deep subseafloor sediments provides a unique opportunity to investigate how genomic diversity, nucleotide substitution, gene expansion/contraction affect the genetic evolution of S. commune over millions of years since they were buried from the ancient forest to lignite-bearing sediment in the subseafloor.

Here, we report a high-quality genome sequencing of S. commune 20R-7-F01 that has been buried in a 1,966 m-deep lignite coal-bearing sediment to advance our insight of S. commune evolution and adaptation to the deep subseafloor extremes. A large-scale single nucleotide polymorphisms (SNP) data from 14 high-coverage genomes of S. commune strains derived from different habitats, including nine (6 subseafloor strains; 1 marine strains, 2 terrestrial strains) newly re-sequenced accessions from this study and five terrestrial accessions from a previous study, was used to investigate the genetic diversity and evolutionary history of S. commune and the mechanism of its adaptation to the depositional environments in the deep subseafloor biosphere.

Results

Genome assembly and scaffolding and annotation

Whole genome of S. commune 20R-7-F01 (CGMCC 11604) was sequenced using a shotgun sequencing strategy and PacBio RSII long-read sequencing platform. A total of 376,337 reads representing a cumulative size of 4.58 Gb were generated with PacBio RSII sequencing (Table S1). The total sequence coverage of the whole genome assembly was ∼110✕. After filtering out the low-quality reads, we obtained a composite 40.79 Mb genome from the retained 3.11 Gb clean reads. The final assembly comprised 162 scaffolds with N50 of 1.83 Mb, with the longest scaffold being 4.5 Mb (Table 1 and Figure 1). The BUSCO analysis suggested that the annotation set was well completed, with 92.8% complete BUSCOs and 3.7% missing BUSCOs. The Quast analysis indicated that the quality of S. commune 20R-7-F01 assembly was similar to that of previously sequenced strain H4-8, but superior to that of the strains LoenenD and TattoneD (Tables S2 and S3).

Table 1

Genome assembly and annotation summary of S. commune strains

Assembly feature	Strain 20R-7-F01	Strain LoenenD	Strain TattoneD	Strain H4-8
Genome size (Mbp)	40.79	35.88	36.46	38.48
Coverage (X)	113.1X	112.3X	109.4X	8.29X
Number of Scaffold	162	1,774	1,707	36
N50 (bp) of contigs	1,826,793	54,148	54,683	2,548,518
GC content (%)	57.32	57.50	57.45	56.67
Gene number	10,765	13,827	15,199	13,210
Average gene length (bp)	1,725	1,708	1,692	1,795
Average exon length (bp)	213	247	264	249
Average intron length (bp)	91	76	72	79
Average number of exons per genes	5.9	5.55	5.27	5.7

Figure 1

The circos diagram of S. commune 20R-7-F01 genome

The outermost layer is the chromosome and its size. The second and third layers are CDS on the positive and negative chains, and the different colors indicate the functional classification of different COGs of the CDS. The fourth and fifth layers are gene density on the positive and negative chains. The sixth and seventh layers are GC content, the green part indicates that the GC content in this area is higher than the whole genome average GC content, and the blue part indicates that the GC content in this area is lower than the whole genome average GC content. Links between all genes represents inparalogs and the bold line indicates the top five genes with the most copies.

Gene annotation on the genome assembly of S. commune 20R-7-F01 was performed using genomic and transcriptomic data. Results showed that the genome comprised 10,765 protein-coding genes, 220 tRNA, 78 rRNA, and 69 other ncRNAs fragments. The average length of gene was 1,725.2 bp, and each gene had 5.9 exons on average (Tables S4 and S5). The average length of exon and intron was 213 bp and 91 bp, respectively. The size of the coding region was 18.57 Mb, accounting for 45.53% of the whole genome. 11.06% of the genome sequences contained repeats and the most frequent repetitive elements were LTRs (1.66%) (Table S6). The content of GC in the genome was 57.32%, and GC in the coding region was 59.31%. The genome size is similar to those observed in the existing S. commune strains isolated from the Earth’s surface, while the number of protein-encoding genes is lower than other terrestrial strains (Table 1).

Genomic comparison between strains 20R-7-F01 and H4-8

To investigate whether subseafloor strain S. commune 20R-7-F01 has unique genetic characteristics, we performed comparative genomic analysis of this strain with the most well-studied terrestrial type strain H4-8. First, a total of 92 gene clusters, including 237 genes, were found only in strain 20R-7-F01, but not in strain H4-8. Except 184 unknown functional genes, other genes are mainly involved in DNA repair, transposon protein, alpha-1,3-glucan synthase, aryl-alcohol dehydrogenase (Table S7). In particular, the number of genes related to DNA replication and repair in strain 20R-7-F01 was 175, 10 more than that in strain H4-8. Second, compared with strain H4-8 (79 bp), strain 20R-7-F01 has larger size of introns (91 bp). Third, a total of 29 gene families, containing 170 genes, mainly encoding FunK1 protein kinase, NADH flavin oxidoreductase, Cytochrome p450 were identified in strain 20R-7-F01, but the copy number of these clusters in strain 20R-7-F01 was at least two more than that of strain H4-8 (Table S8). Fourth, both strains shared 152 conserved syntenic blocks (30.82 Mb), which was accounting for 76 and 80% of the genome size in strain 20R-7-F01 and strain H4-8, respectively. However, a 0.16 Mb inversion occurred on unitig_432 of strain 20R-7-F01, corresponding to the scaffold_7 of strain H4-8 (Figures S1 and S2). In addition, strain H4-8 had 11,343 inserted duplications and 1,021 translocations, while strain 20R-7-F01 had only 9,640 inserted duplications and 958 translocations (Table S9). Fifth, according to Ka/Ks ratio of 16,474 collinear genes, we identified a total of 525 genes that provide evidence of positive selection in strain 20R-7-F01 relative to strain H4-8, which are mainly involved in FunK1 protein kinase, glycosyl hydrolase, as well as regulation of Rho protein signal transduction (Table S10).

Genome evolution and divergence time estimates

To determine the evolutionary relationships of S. commune 20R-7-F01 with other closely related fungal isolates, the protein sequences of 11 species/strains were analyzed with a gene family clustering method. A total of 957 single-copy orthologous genes were selected to construct a phylogenomic tree with the maximum likelihood method. Phylogenomic analysis showed that the species in Agaricales were clustered into a branch, of which S. commune 20R-7-F01 and S. commune H4-8 were a clade. Based on the phylogeny and fossil records, we suggest that S. commune 20R-7-F01 diverged from the reference S. commune H4-8 genome about 28–73 Mya (Figure 2). A Mitochondrial phylogenetic tree also confirmed that the two strains differentiated about 4.8–22 Mya (Figure S3), which was consistent with the estimated ages of sediment (approx. 20 Mya) from which strain 20R-7-F01 were isolated (Gross et al., 2015). In addition, compared with other fungal species, the subseafloor S. commune 20R-7-F01 contained 14 significantly expanded gene families encoding FunK1 kinase, NmrA family, transposon protein, and carboxy-cis, cis-muconate cyclase, as well as two contracted gene families encoding pepsin-like aspartate proteases and integrase (Figure 2 and Table S11).

Figure 2

Phylogenetic tree and divergence time of S. commune 20R-7-F01

The branch lengths of the phylogenetic tree are scaled to estimated divergence time. The blue bars on the nodes indicate the 95% credibility intervals of the estimated posterior distributions of the divergence times. The overall timeline is shown below the phylogenetic tree.

Genetic relationship and population structure

To reveal the evolutionary process of S. commune using a large-scale SNP and insertion/deletion (InDel) dataset generated from 14 high-coverage genomes of S. commune, including nine newly re-sequenced accessions and five previously studied terrestrial accessions. Approximately 1.77Gb of high-quality reads were collected using Illumina whole genome resequencing of the nine strains, which was aligned against the reference genome with a mapping rate ranging from 73.39 to 96.55%, resulting in coverage depths ranging from 35 to 71✕(Table S12). Combining our data with the five previously sequenced strains (all with >50✕ coverage depth), we identified a total of 19,755,267 SNPs and 3,753,725 InDels in the 14 accessions. Among them, 5,595,289 SNPs and 1,009,532 InDels were found in subseafloor strains, accounting for 39.53 and 36.79% of terrestrial and marine strains, respectively (Table S13).

Using whole-genome SNP data, we performed a neighbor-joining (NJ) phylogenetic analysis and principal component analysis (PCA) on the 31 samples (15 strains, 16 basidiospores) from different habitats. The NJ tree clearly showed that strains collected from subseafloor sediments notably differ from those obtained from terrestrial habits (Figure 3A). The ITS sequence analysis of 127 strains from NCBI databases further confirmed the differentiation between subseafloor and terrestrial strains of S. commune (Figure S4). All these results showed that subseafloor strains of S. commune are a distinct subset that do not overlap with other sampled diversity. In addition, the marine strain MCCC 3A00233 was clustered in the terrestrial group, suggesting that this strain may have originated from the land. In addition, the results of PCA analysis are consistent to those of NJ tree analysis. Although strains MF and Hom2-8 were relatively far from other terrestrial strains, the first two principal components, PC1 and PC2, could clearly distinguish the subseafloor and terrestrial strains (Figure 3B).

Figure 3

Phylogenetic, population stratification and principal component analyses of S. commune

(A) Maximum likelihood phylogenetic tree of subseafloor, marine, and terrestrial S. commune population.

(B) PCA analysis for 30 S. commune samples.

Genetic differentiation and diversity

To understand whether loss-of-function polymorphisms (LoF) such as frameshift and premature stop codon play a selective role in the environmental adaptation of S. commune, we calculated these polymorphisms based on SNP datasets. A total of 684 and 2,874 premature stop codons that were caused by SNVs (stop-gain), 18,429 and 66,327 insertion/deletion-induced frameshift variants (frameshift) leading to the disruption of a transcript reading frame, and 101,371 and 392,083 nonsynonymous SNVs were identified in subseafloor and terrestrial strains of S. commune, respectively (Table S13). Taking all these stop-gain, frameshift, and nonsynonymous SNVs variants together, we identified 120,484 and 361,284 LoF variants across 9,362 and 10,453 protein-coding genes in subseafloor and terrestrial strains of S. commune, respectively. This indicates more unexpected redundancy in terrestrial strains than in subseafloor strains. In addition, 29 LoF genes were only found in the subseafloor strains of S. commune, among which two genes were related to transposon protein (Table S14), and the other 1,141 LoF genes were not found in the subseafloor population, involving in structural constituent of ribosome, pyrophosphatase activity, hydrolase activity and ATPase activity (Table S15).

Moreover, genes with LoF mutations showed a correlation with nucleotide diversity (π). The average nucleotide diversity (π) in LoF genes (subseafloor population = 0.0201; terrestrial population = 0.0196) was higher than that in non-LoF genes (subseafloor population = 0.0141; terrestrial population = 0.0112) (Table S16).

Roles of mutation and recombination

To explore the relative effect of nucleotide diversity originating from mutations versus homologous recombination (r/m) on the genetic diversification of populations, we calculate the relative rate of recombination to mutation (R/θ), the mean length of recombined DNA (δ), and the mean divergence of imported DNA (v) for branch tips in S. commune genome phylogeny, which allows for a calculation of r/m (r/m = (R/θ) ✕δ✕ v). The result showed that the r/m of subseafloor strains was 0.0526, which was approximately 5,000 times lower than that of terrestrial strains (r/m = 266.66) (Table 2). To avoid the impact of different geographical locations on the analysis results, we choose the small-scale geographical and clonal terrestrial strains (5 strains from Moscow and 4 strains from Florida) (Baranova et al., 2015) to compare with our subseafloor strains by using ClonalFrame software. The result showed that the r/m of Moscow strains was 145.75, and the r/m of Florida strains was 71.10, which were approximately 2,770 and 1,351 times higher than that of subseafloor strains (r/m = 0.0526), respectively (Table S17). To investigate whether the heterozygosity of dikaryotic strains 6R-2-ZF01 (2.96%), 15R-5-ZF01 (3.05%) and 24R-3-ZF01 (2.98%) (Table S13) had an impact on the calculation of r/m value, we compared the r/m value of subseafloor homozygous samples (0.0139) with those of both homozygous and heterozygous samples (0.0526) (Table S18), and found that the r/m value of subseafloor strains was significantly lower than that of terrestrial strains (r/m = 266.66). In addition, we also calculated the index of association (IA) of subseafloor strains and Moscow strains, and found that the frequency of IA in subseafloor strains was concentrated around 0.55, which was higher than that of Moscow strains (0.03) (Figure S5). These consistent results indicated that homologous recombination played a minimal role in diversification of the subseafloor strains.

Table 2

The contribution of recombination and mutation to nucleotide diversity of subseafloor and terrestrial S. commune populations

Group	R/θ	δ	v	r/m
Subseafloor	0.00238761	2084.436	0.0105791	0.0527
Terrestrial	0.122791	11831.35	0.183557	266.67

To estimate genomic diversity of subseafloor and terrestrial strains, we determined the nucleotide diversity (π) and Tajima’s D based on SNP variants of 10 re-sequenced strains (Table S19). The results showed that the average nucleotide diversity (π = 0.01758) of the subseafloor strains was approximately ½ that (π = 0.0333) of terrestrial strains, and the nucleotide diversity between the subseafloor strain 20R-7-F01 and other terrestrial strains ranged from 0.068 to 0.104 after Jukes-Cantor corrections (Table S20). Compared with the average π value of homozygous samples (0.01356), the average π value of samples containing both homozygous and heterozygous strains was 0.01758, which was much lower than that (π = 0.0333) of terrestrial strains (Table S21), indicating that the heterozygosity of dikaryotic strains did not affect the results of genomic diversity analysis. Moreover, the subseafloor strains also had a lower ratio of nonsynonymous/synonymous mutation (nonsyn/syn = 0.22) and transitions to transversions (Ts/Tv = 1.93) than that of terrestrial strains (nonsyn/syn = 0.24, Ts/Tv = 1.98), probably due to the extremely slow growth or reproduction of subseafloor strains in situ (Table S13).

Selective sweep analysis

Since selective sweep regions usually contain loci related to environment adaptation, we performed a whole-genome screening of the overlapping selective sweep regions by combining Tajima’s D with paired F-statistics (Fst) analysis. A total of 104 and 244 selective sweep regions with 0.25 Mb (153 gens) and 0.68 Mb (325 genes) genome size were identified in the subseafloor (top/bottom 10% for Tajima’s D2.82/0.88) and terrestrial groups (top/bottom 10% for Tajima’s D = 1.55/-0.49), Fst = 0.4540, respectively. There were 97 genes in the selective sweep regions shared by the subseafloor and terrestrial groups, and 56 and 228 unique genes, respectively. The specific genes in subseafloor group are associated with the activities of selective transposase (4 genes) and hydrolases (2 genes) (Table S22).

Discussion

S. commune is a ubiquitous white rot fungus with a worldwide distribution and has been utilized as a model system for studying mating-type gene function and mushroom development (Ohm et al., 2010). Although the genome of 58 terrestrial strains of S. commune has been sequenced to date, only strain H4-8 has been fully annotated (Ohm et al., 2009). Here, we presented for the first time a high-quality genome of S. commune 20R-7-F01, which can be used as a model strain of S. commune subsisting deep subseafloor environment over 20 Mya for the study of its origination, evolution, and environmental adaptation mechanism by comparison with terrestrial strains. The subseafloor strain contained 10,765 protein coding genes, of which 13.86% had no ortholog in the strain S. commune H4-8, indicating that subseafloor S. commune are genetically different from their terrestrial counterparts.

Compared with terrestrial strains H4-8, the subseafloor strain 20R-7-F01 has a similar-sized genome, but longer introns and more numbers of genes associated with DNA repair and transposon proteins. The genome size of the subseafloor strain 20R-7-F01 was similar to that of the terrestrial strains, suggesting that hetero recombination or lateral gene transfer played a small role in driving the evolution of the S. commune during ∼20 million years of depositional history (Nelson-Sathi et al., 2015). It has been reported that the high fidelity of DNA replication and repair mechanisms plays an important role in genome conservation of subsurface microbes (Becraft et al., 2021). As the in situ temperature, where strain 20R-7-F01 was isolated, was about 45–50°C, it generally causes damage to DNA, proteins and other biomolecules (Inagaki et al., 2015; Lever et al., 2015; Steen et al., 2013). Therefore, the subseafloor strain possess extra genes related to DNA repair as well as large introns that buffer and regulate fluctuations in protein concentration in fungal cells and protect coding DNA from mutations (Lynch and Conery, 2003; Stival Sena et al., 2014), which might have helped the fungus to survive there and play an essential role in the genome conservation. In addition, five genes encoding transposon proteins in subseafloor S. commune 20R-7-F01 were notably expanded compared with other fungi, which indicated that transposon proteins may play an important role in genome plasticity and effective means for organisms to adapt to stressful or extreme environments (Vigil-Stenman et al., 2017; Li et al., 2014).

Meta-omics analyses revealed that subseafloor fungi are not only metabolically active but may play important ecological roles in the subseafloor environments (Quemener et al., 2020; Ivarsson et al., 2018). Little is known about the mechanisms by which fungi adapt to extreme conditions of the subseafloor; however, limited research suggested that fungi sustain life by hydrolyzing their own cell walls (autolysis) under carbon starvation conditions (Emri et al., 2005; Kim et al., 2011). Our previous work showed that the subseafloor strains were able to utilize lignite/lignin as a sole carbon and energy source (Zain Ul Arifeen et al., 2020), possessed normal growth ability under anaerobic conditions by employing various strategies, including ethanol fermentation and production of amino acids, increase in mitochondrial number, and autophagosome formation (Zain Ul Arifeen et al., 2021a, 2021b). In this study, we also found that subseafloor strain 20R-7-F01 contained more CAZymes coding genes related to cell wall hydrolysis than terrestrial strain H4-8, and these genes may help the fungus obtain energy by degrading the buried organic matter (from terrestrial forest soil to lignite coal) in the sedimentary system (Inagaki et al., 2015). In addition, genes encoding FunK1 kinase, NmrA family, and carboxy-cis, cis-muconate cyclase were also found to be expanded in S. commune 20R-7-F01 compared to the other fungus, suggesting that NmrA family genes play a vital role in regulating the activity of the GATA transcription factor AreA during nitrogen metabolism in various fungi (Stammers et al., 2001), and indicating that carboxy-cis, cis-muconate cyclase is a major enzyme for the dissimilation of aromatic compounds derived from lignin (Thatcher and Cain, 1975). On the contrary, genes encoding retropepsins, pepsin-like aspartate proteases, and integrase were significantly contracted in S. commune 20R-7-F01, indicating that S. commune 20R-7-F01 could be barely infected by a virus since it isolates from the subseafloor sedimentary environment.

LoF mutations are known to play important roles in adaptation and phenotypic diversification as well as rewire the cell’s metabolism (Hottes et al., 2013; Xu et al., 2019), also found that the presence of LoF mutations is correlated with the level of nucleotide diversity, the density of transposable elements (TEs), and gene family size (Xu et al., 2019). However, very little is known about the evolutionary patterns of LoF mutations at the genomic level and the adaptive effects of LoF variants at the population level (Albalat and Cañestro, 2016) in the deep biosphere. Our study first revealed that the frequency and quantity of LoF mutations in subseafloor sedimentary strains of S. commune were less than that of terrestrial strains, possibly due to the extremely slow growth and/or proliferation of fungi in the lignite coal-bearing strata, with minimal estimates of microbial turnover time ranging from a few months to over 100 years under the laboratory experimental conditions (Trembath-Reichert et al., 2017). In addition, two LoF genes, related to transposon proteins, were only found in the population of subseafloor strains, which indicated that transposon proteins may play an important role in genome plasticity and effective means for organisms to adapt to stressful or extreme environments (Vigil-Stenman et al., 2017; Li et al., 2014).

The ratio of homologous recombination and nonsynonymous/synonymous mutation can be used to characterize genome evolution (Vos and Didelot, 2009). In this study, we found that the recombination rate and nonsynonymous to synonymous substitution rate of subseafloor strains of S. commune were notably reduced compared with that of terrestrial strains. Moreover, similar results were suggested in the subseafloor bacterial population, which had a lower homologous recombination rate than the terrestrial bacterial populations due to the physical isolation of individual cells, reduced cell concentrations, and reduced availability of extracellular DNA for recombination in subseafloor sediments (Lever et al., 2015). Thus, the limited availability of nutrients (Quemener et al., 2020), slow growth/proliferation of hyphae/spores (Pachiadaki et al., 2016; Raghukumar et al., 2004), and few cell numbers (Lever et al., 2015; Quemener et al., 2020) may be the main reasons for the decrease in the genetic diversity of S. commune strains in deep subseafloor sediments compared with those on land of this modern world.

Although the strain 20R-7-F01 was isolated from core sediment Unit III and the geological age of this sediment has been determined to be ∼20 Mya (Gross et al., 2015), the actual age of the fungus is however difficult to determine. There are great differences in the age of subseafloor fungi reported; for instance, the fossilized fungi isolated from the deep subseafloor basalts (Ivarsson, 2012) and the continental igneous crust (Drake et al., 2019, 2021) had been dated as ∼81-48 and ∼39 Mya, respectively. However, the age of these fungi might theoretically be much younger than the age of their habitats where active hydrothermal circulation may occur (Ivarsson, 2012; Drake et al., 2019, 2021). This suggests that fungal age based on geological age of the sediments might not represent the accurate age of these anciently deposited fungi to the subseafloor sediment. Interestingly, based on genomic and mitochondrial phylogenetic analysis we found that the divergence time of 20R-7-F01 and H4-8 was 28–73 Mya and 4.8–22 Mya, respectively, which is consistent with the geological age of the sediment (Gross et al., 2015; Inagaki et al., 2012). Although, we could not rule out all possibilities about the origin of fungi in subseafloor sediments, multiple evidences from structural geology, sedimentology, (bio)geochemistry, and microbiology consistently support that subseafloor fungi are indigenous and have persisted in subseafloor sediments over millions of years (Inagaki et al., 2015; Liu et al., 2017; Zain Ul Arifeen et al., 2021a). However, the arrival of these fungi to subseafloor sediments through underground circulation is not feasible because the occurrence of active fluid circulation is not physically possible under the 2 km-deep subseafloor sedimentary system. In particular, the reverse advection from shallow to deep below the unconformity layer (marine-terrestrial transition) associated with carbonates and conglomerates will not be possible (Inagaki et al., 2012).

Hypervariable species are particularly interested for population genetics (Baranova et al., 2015; Seplyarskiy et al., 2014), and population genetics can examine the distribution of genetic variation and levels of genetic diversity within and between populations (Keats and Sherman, 2014). Genetic diversity of a population is both a factor and an outcome of a wide range of evolutionary processes. In the vast majority of the species, nucleotide diversity (π), the probability that two alleles randomly chosen from the population at a nucleotide site differ from each other, is below 0.03–0.05 (Leffler et al., 2012). However, our results show that S. commune exceeds this value, for examples, the nucleotide diversity of the subseafloor strain 20R-7-F01 and other 24 terrestrial haploid strains (Baranova et al., 2015; Seplyarskiy et al., 2014) ranged from 0.07 to 0.11, also the π of the subseafloor strain 20R-7-F01 and the terrestrial strain H4-8 reached about 0.104. These results are consistent with other studies where the diversity at synonymous sites is 0.2 in the American population of S. commune and 0.13 in the Russian population, thus confirming that S. commune is the most polymorphic among all studied eukaryotic species (Baranova et al., 2015; Seplyarskiy et al., 2014).

Conclusions

In this study, we present the first de novo whole genome sequencing and assembly of S. commune 20R-7-F01, isolated from a ∼20 millionyearsago coal-bearing sediment at ∼2.0 km below the ocean floor. Phylogenomic analysis reveals that the divergence between strain 20R-7-F01 and terrestrial type-strain H4-8 is approximately 28–73 Mya. The expansion of the genes encoding FunK1 kinase, NmrA family, transposon protein, and Carboxy-cis, cis-muconate cyclase possibly contributed to the extreme environmental adaptation of S. commune from the coastal to the deep subseafloor. Nucleotide diversity, substitution rate, and the homologous recombination rate were lower in subseafloor strains than those in terrestrial strains, reflecting that the growth or reproduction of subseafloor strains is extremely slowly. Our data contribute to the study of S. commune evolutionary trajectory and phylogenomic diversity among the species and provide novel targets in understanding the molecular adaptation mechanisms of the fungi underlying the geological perturbations associated with the deep subseafloor sedimentary biosphere.

Limitations of the study

Although we report here new insights into the evolutionary history and adaptation mechanism of S. commune through comparative genomic analyses, more detailed studies, such as large-scale population genome analysis, will be needed in the future to identify vital gene loci and physiological characteristics associated with environmental adaptation. Also, the adaptation and survival mechanisms of deep subseafloor microbial communities along with increasing depth, pressure, and temperature during depositional regimes (e.g., Heuer et al., 2020; Köster et al., 2021; Beulig et al., 2022) could be addressed by scientific ocean drilling and mimicking cultivation-based laboratory experimentations in the future.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests should be directed to the Lead Contact, Changhong Liu (chliu@nju.edu.cn).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Fungal strains and culture conditions

A total of nine S. commune strains were used for sequencing in this study. Six strains, including three dikaryotic (6R-2-ZF01, 15R-5-ZF01, 24R-3-ZF01) and three corresponding monokaryotic (6R-2-F01, 20R-7-01, 24R-3-F01) strains, isolated from 1496 ∼1993 mbsf sediments during IODP Expedition 337 (Liu et al., 2017). Different subseafloor strains come from different cores, different distances, and different geological ages; for example, subseafloor strain 6R-2-F01 was isolated from sediment core Unit 2 at a depth of 1,496 mbsf with an estimated geological sediment age of Miocene, while strains 15R-5-F01/20R-7-F01/24R-3-F01 were derived from sediment core Unit 3 at depths of 1,924 mbsf, 1966 mbsf and 1,993mbsf; respectively, having geological age of early to late Miocene (Gross et al., 2015). Two terrestrial strains CFCC 7252 and CFCC 86625 were purchased from China Forestry Culture Collection Center, which were isolated from Populus wood in Songshan, Beijing (altitude:700m) and Jurong, Jiangsu (altitude: 300m) of China, respectively, and the other five terrestrial strains (225DK, 227DK, MF, Hom2-8, 207) were obtained from NCBI and JGI database. Moreover, one strain MCCC 3A00233, collected from marine sediment of the Atlantic Ocean, was purchased from Third Institute of Oceanography, State Oceanic Administration, People’s Republic of China. Details of habitat of marine and terrestrial strains are shown in Table S12. Details of the habitat and culture conditions of S. commune strains have been described previously (Zain Ul Arifeen et al., 2020). For DNA and RNA isolation, mycelium of S. commune was grown in 250 mL conical flask containing 150 mL PD (200 g/L potato and 20 g/L glucose) and incubated at 30°C and 200 rpm for 3–5 days.

Method details

De novo genome sequencing and assembly

Genomic DNA of S. commune 20R-7-F01 (CGMCC 11604) was extracted using a CTAB-PEG method (Huang et al., 2018). An Illumina paired-end DNA library of 20 kb insert size was generated using the TruSeq DNA Sample Preparation Kit (Illumina, USA) and the Template Prep Kit (Pacific Biosciences, USA). Genome sequencing was performed using the PacBio RSII long-read sequencing platform at BGI Genomic (Shenzhen, China). A total of 4.58 Gb of raw PacBio data (∼113.1 ✕ genome coverage) was obtained, of which 1.47 Gb was removed due to low complexity (length <100 bp), low quality (score <0.8), and adapter and duplication contamination (subreads length <1000 bp; score <0.8). The resultant clean PacBio data were assembled using SMRT Analysis (v2.3.0) HGAP3 workflow to construct scaffold-level contigs (Chin et al., 2013). Integrity of the genome assembly was assessed using BUSCO software (version v3.0.2) (Waterhouse et al., 2018). Quast software was used to evaluate the assembly quality of S. commune 20R-7-F01 genome (Gurevich et al., 2013).

Transcriptome sequencing and assembly

Total RNA was extracted from mycelium of S. commune 20R-7-F01 using an RNAprep Pure Plant Kit (TIANGEN, Cat. #DP441, China) following the manufacturer’s instructions. Magnetic oligo (DT) beads were used to enrich polyA mRNA tails of three independent RNA samples. The cDNA libraries were sequenced using the Illumina HiSeq™ 2000 platform by Personal Biotechnology Company (Shanghai, China). After trimming of low quality reads (Q<20) and adapter contamination, ∼72.5 Gb high-quality transcriptome data reads were obtained, which was mapped to genome assembly of strain 20R-7-F01 using TopHat (Trapnell et al., 2009), and performed for gene prediction using Cufflinks (Roberts et al., 2011).

Gene annotation

The repetitive DNA in S. commune 20R-7-F01 genome was identified using RepeatModeler and RepeatMasker (http://www.repeatmasker.org/). Mitochondrial sequences were identified by BLASTN search using published mitochondrial genome of S. commune (GenBank ID: AF402141.1) (Forget et al., 2002) and were excluded when annotating genes in the nuclear genome. The rRNA, tRNA, and ncRNA were predicted with RNAmmer (Lagesen et al., 2007), tRNAscan-SE (Lowe and Eddy, 1997), and Rfam (Gardner et al., 2009), respectively. Open reading frames were predicted with Fgenesh (Salamov and Solovyev, 2000), Exonerate (Slater and Birney, 2005), Cufflinks, EVidenceModeler (Haas et al., 2008), based on both assembly and RNA-seq data. Gene functions were annotated using NCBI non-redundant protein database (Nr), Eukaryotic Clusters of Orthologous Groups (KOG), SwissProt, Gene Ontology (GO), KEGG, TIGRFAM, and PFAM databases (Li et al., 2021). Mitochondrial genes were predicted according to the annotation of S. commune (GenBank ID: AF402141.1) after aligning the local mitochondrial sequences to the published mitochondrial genome using BLASTN.

Gene family identification and phylogenetic evolution analysis

A total of 10 fungal species (i.e., 11 strains) assigned to Basidiomycota including S. commune (two strains 20R-7-F01 and H4-8), Coprinopsis cinerea, Laccaria bicolor, Postia placenta, Phanerochaete chrysosporium, Cryptococcus neoformans and Ustilago maydis, and Ascomycota such as Saccharomyces cerevisiae, Aspergillus nidulans and Neurospora crassa were used to perform gene family identification. Protein sequences of these species except S. commune 20R-7-F01 were downloaded from NCBI and JGI databases, and analyzed by OrthoFinder with default parameters to find orthologous genes (Emms and Kelly, 2015). The orthologous genes between S. commune 20R-7-F01 and S. commune H4-8 were analyzed using Orthomcl software (Li et al., 2003). The single-copy genes of the 11 analyzed strains were aligned using MAFFT to create a super-gene for each species (Katoh and Standley, 2013). After Gblocks alignment optimization (Castresana, 2000), the conserved blocks of super-genes were used for phylogenetic tree construction using RAxML (Stamatakis, 2014). The divergence time was estimated using the “Independent rates (clock = 2)” model and the “HKY85” model of MCMCTREE program in the PAML package (Yang, 2007), based on the fossil calibration of the most recent common ancestors (MRCA) of N. crassa and A. nidulans (divergence at 290-470 Mya), and L. bicolor and C. cinerea (divergence at 120-180 Mya) (Kumar et al., 2017). The MCMC process was run for 1,000,000 iterations, after a burn-in of 2,000 iterations, and the convergence was evaluated by ESS, which was calculated with Tracer software. The mtDNA protein sequences of A. nidulans, C. cinerea, C. neoformans, L. bicolor, N. crassa, S. cerevisiae, S. commune and U. maydis were downloaded from GenBank (accession numbers JQ435097.1, AACS02000068.1, CP003834.1, MK697670.1, KC683708.1, KP263414.1, AF402141.1, and DQ157700.1, respectively) and also analyzed by OrthoFinder with default parameters to find orthologous genes. Finally, three single copy homologous genes (cytochrome c oxidase subunit 1, cytochrome c oxidase subunit 2, cytochrome c oxidase subunit 3) were identified. The cytochrome c oxidase gene families, as molecular markers, always widely used in molecular evolutionary and phylogenetic relationship studies within and between fungus (Belloch et al., 2000; Okane et al., 2021). Additionally, the same method was used to construct the mitochondrial phylogenetic tree. The expansion and contraction of orthologous gene family were calculated with CAFE and defined as “significantly expanded or contracted” if p-value ≤0.01 (De Bie et al., 2006). The positive selected genes (PSGs) between S. commune 20R-7-F01 and S. commune H4-8 were detected using the MCScanX (i.e., add_ka_and_ks_to_collinearity.pl) (Wang et al., 2012). In addition, synteny analysis of strains 20R-7-F01 and H4-8 was performed using MUMmer sequence alignment package (Delcher et al., 2002). Conserved syntenic blocks between the two strains were identified using the MCScanX package with default parameters based on a minimum requirement of five collinear orthologous genes

Genome resequencing, alignment and variant calling

The genome DNA of S. commune strains (Table S12) was extracted and fragmented to generate approximately 300 bp library insert size. The paired-end reads were sequenced on an Illumina HiSeq 2500 platform at BGI Genomic (Shenzhen, China). Raw data were filtered by removing the adapters and low-quality reads (Q<20). All the filtered re-sequencing reads were mapped to the reference genome of S. commune 20R-7-F01 using BWA-MEM with default parameters (Li and Durbin, 2009). After removing reads with low mapping quality (MQ < 40), both paired-end and single-end mapped reads were used for SNP detection across entire samples of S. commune strains using the GATK toolkit (version 4.0.8.1). SNPs and small InDels (insertion/deletions, 1–50 bp) were called using GATK software (DeSumma et al., 2017), filtered using SelectVariants according to the standard hard filtering criteria: “QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0” (DePristo et al., 2011), and assigned to specific genomic regions and genes using Annovar based on 20R-7-F01 genome annotations (Wang et al., 2010). Additionally, bcftools (Danecek et al., 2021) was used to calculate homozygous and heterozygous sites with “snp.vcf” file as input, and estimated the heterozygosity as “number of heterozygous SNPs”/”genome size”. Loss-of-Function (LoF) variants were identified including mutations of stopgain, nonsynonymous SNV, and frameshift (Xu et al., 2019).

Phylogenetic and population structure analyses

A subset of 2,704,270 SNPs obtained from resequencing of 14 S commune strains was used for phylogenetic and population structure analysis. Vcf files were converted to hapmap format with custom perl scripts and to PLINK format file by PLINK (Purcell et al., 2007). Under the p-distances model, a neighbor-joining tree of all samples was constructed with SNPhylo software (Lee et al., 2014). GCTA was used to carry out principal component analysis (PCA), and the first 2 eigenvectors were extracted (Yang et al., 2011).

VCFtools (Danecek et al., 2011) was used to estimate the Fst, nucleotide diversity and Tajima’s D for the whole genome of both subseafloor and terrestrial populations with non-overlapping 10 kb sliding windows. The top 10% of Fst with top/bottom 10% of Tajima’s D, defined by a sliding window, were selected as candidate selective sweep regions. We further integrated the adjacent selective sweep regions and annotated the corresponding genes across genome of strain 20R-7-F01.

Detection of homologous recombination

The contributions of mutations and recombination to the genomic diversity in the concatenated core genome alignment, the number of recombination events (imports) per genome, and the positions of recombination hot spots, was investigated using ClonalFrameML (Didelot and Wilson, 2015). Nucleotides unaffected by recombination are referred to as unimported and nucleotides subject to recombination are referred to as imported. ClonalFrameML provides the relative rate of recombination to mutation (R/Theta), the mean length of recombined DNA (Delta), and the mean divergence of imported DNA (Nu). These results were used to calculate the relative contribution of recombination versus mutation to the overall genomic diversity (r/m), using the following formula r/m = (R/Theta) ×Delta × Nu. ClonalFrameML was performed in two separate runs containing a core genome alignment that contained (1) only the sediment genomes, and (2) only the terrestrial strains. The resulting r/m values from these two groups (presented in Table 2) were then used to interpret the relative importance of mutations compared to recombination, in the separate groups (e.g., terrestrial strains versus subseafloor strains). Additionally, the same method was used to calculate the homologous recombination of small-scale geographical terrestrial strains (5 strains from Moscow and 4 strains from Florida) (Baranova et al., 2015), which could avoid the impact of different geographical locations on the analysis results. We acknowledge that because the dataset contains genomes covering the diversity of a single fungal genus (Schizophyllum), the only detectable recombination most events are from donors from the species under study, so that the main source of recombination is not external. Moreover, we also calculated the index of association for subseafloor population and Russia population by using a function in R package “Poppr” (Kamvar et al., 2014, 2015). The function is named “samp.ia”, which can randomly sample sites to calculate the index of association. The specific parameter settings are as follows: the number of snps to be used to calculate standardized index of association was set 100L; the number of times to perform the calculation was set 100L.

Quantification and statistical analysis

Quantification and statistical analysis used in the genome sequencing and assembly, genome quality assessment, evolutionary analysis and comparative genome analysis can be found in the relevant sections of the method details.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Total DNA and RNA of Schiophyllum commune	This study	NA
Deposited data
original code	This study	https://github.com/liuxuan-425lab/genome_code
Resequencing data	This study	GeneBank: PRJNA738972
RNA sequence data	This study	GeneBank: PRJNA543698
Genome assembly data	This study	GeneBank: PRJNA544166
Other genome data	Baranova et al. (2015)	GeneBank: PRJNA236351
Other resequencing data	Seplyarskiy et al. (2014)	GeneBank: PRJNA234274
Software and algorithms
BUSCO	Waterhouse et al. (2018)	https://busco.ezlab.org/
HGAP3	Chin et al. (2013)	http://www.pacb.com/devnet/
Quast	Gurevich et al. (2013)	http://bioinf.spbau.ru/quast
TopHat	Trapnell et al. (2009)	http://ccb.jhu.edu/software/tophat/index.shtml
Cufflinks	Roberts et al. (2011)	http://cole-trapnell-lab.github.io/cufflinks/
RNAmmer	Lagesen et al. (2007)	http://www.cbs.dtu.dk/cgi-bin/nph-sw_request?rnammer
tRNAscan-SE	Lowe and Eddy (1997)	https://github.com/UCSC-LoweLab/tRNAscan-SE
Fgenesh	Salamov and Solovyev (2000)	http://linux1.softberry.com/berry.phtml?topic=fgenesh&group=programs&subgroup=gfind
EVidenceModeler	Haas et al. (2008)	https://evidencemodeler.github.io/
OrthoFinder	Emms and Kelly (2015)	https://github.com/davidemms/OrthoFinder
Orthomcl	Li et al. (2003)	https://orthomcl.org/orthomcl/app/downloads/software/
MAFFT	Katoh and Standley (2013)	https://github.com/GSLBiotech/mafft
Gblocks	Castresana (2000)	https://github.com/atmaivancevic/Gblocks
RAxML	Stamatakis (2014)	http://evomics.org/learning/phylogenetics/raxml/
PAML	Yang (2007)	http://abacus.gene.ucl.ac.uk/software/paml.html
CAFE	DeBie et al. (2006)	https://github.com/hahnlab/CAFE
MUMmer	Delcher et al. (2002)	http://mummer.sourceforge.net/
BWA-MEM	Li and Durbin (2009)	https://github.com/lh3/bwa
GATK	De Summa et al. (2017)	https://gatk.broadinstitute.org/hc/en-us
Annovar	Wang et al. (2010)	http://www.openbioinformatics.org/annovar/annovar_download_form.php
PLINK	Purcell et al. (2007)	http://zzz.bwh.harvard.edu/plink/
SNPhylo	Lee et al. (2014)	https://github.com/thlee/SNPhylo
GCTA	Yang et al. (2011)	https://yanglab.westlake.edu.cn/software/gcta/
VCFtools	Danecek et al. (2011)	http://vcftools.sourceforge.net/
ClonalFrameML	Didelot and Wilson (2015)	https://github.com/xavierdidelot/ClonalFrameML
MCScanX	Wang et al. (2012)	https://github.com/wyp1125/MCScanX
bcftools	Danecek et al. (2021)	http://github.com/samtools/bcftools
Poppr	Kamvar et al. (2014)	https://github.com/grunwaldlab/poppr