Anaerobic gut fungi (phylum Neocallimastigomycota) are native to the rumen and hindgut of large herbivores where they are well known for their role in plant cell wall degradation. Anaerobic gut fungi are of significant biotechnological interest due to their production of abundant carbohydrate-active enzymes (CAZymes) and multienzyme complexes, termed fungal cellulosomes (Haitjema et al., 2017), which break down Abbreviations: CAZyme, carbohydrate-active enzyme; N. californiae, Neocallimastix californiae; RINe, RNA Integrity Number equivalent; RQN, RNA Quality Number; JGI, Joint Genome Institute; TPM, transcripts per million; KOG, Eukaryotic Orthologous Groups of proteins; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; IPR, InterPro; SMURF, Secondary Metabolite Unknown Regions Finder; SWEET, Sugars Will Eventually be Exported Transporters; TCDB, Transporter Classification Database; PFL, pyruvate formate lyase; q, adjusted p-value from DESeq2 analysis; p, p-values from Fisher's Exact Tests; GH, glycoside hydrolase; CBM, carbohydrate-binding module; GT, glycosyltransferase; EC, number Enzyme Commission number; SucCoA, succinyl coA; PKS, polyketide synthase; NRPS, non-ribosomal peptide synthetase. lignocellulose, enabling applications in the industrial conversion of abundant lignocellulosic biomass (1.3 billion tons per year) to fuels and chemicals (Baruah et al., 2018). These early-branching fungi also possess a wealth of untapped biosynthetic gene clusters encoding predicted secondary metabolites, which have potential applications in pharmaceuticals and bio-based products such as biopolymers and biofuels (Bhattarai et al., 2021;Moeini et al., 2021;Swift et al., 2021;Zhang and Keasling, 2012).
However, anaerobic fungi feature a zoosporic life cycle that has an unknown influence on their more well known degradative and productive properties (Fig. 1). Free-swimming, flagellated fungal zoospores encyst into fibrous substrate and develop into fungal sporangia, which are vegetative, sac-like structures that form new zoospores. The sporangia of some strains grow rhizoids, root-like structures that branch into substrate and form mats. Ultimately, sporangia develop to maturation and release new zoospores (Bauchop, 1981;Elshahed et al., 2022;Liggenstoffer et al., 2010). Anaerobic fungi likely also possess an understudied aerotolerant stage that may mature through either fungal zoospores or sporangia (Wubah et al., 1991). To date, nearly all studies conducted on the anaerobic fungi are performed in "bulk", meaning that any sequencing efforts are reflective of an unknown mix of these populations, which likely varies in time and are influenced by cultivation conditions.
Until recently, genetic engineering tools have been non-existent for anaerobic gut fungi, hindering efforts to test gene function, construct mutant strains, or test hypotheses related to metabolic function. Some genetic engineering efforts targeted towards the anaerobic gut fungi aim to selectively transform zoospores, due to the expected greater ease of DNA transformation into zoospores vs. mature cells with thick cell walls (Calkins et al., 2018;Durand et al., 1997), while other efforts aim to genetically engineer multiple life stages (Hooker et al., 2023). Better understanding of the life stages of anaerobic gut fungi can improve our ability to design effective genetic engineering tools and mine useful enzymes and metabolites for biotechnological applications. Moreover, separation of these life stages can help aid in our basic understanding of early-branching fungal development. In anaerobic gut fungal genomes, a significant portion (~70 %) of open reading frames lack functional annotation, and predicted functions of many other genes are not yet verified (Haitjema et al., 2014;Henske et al., 2018a;Solomon et al., 2016a). In this study, we develop hypotheses for the roles of some of these genes, particularly those involved in fungal maturation and lifestyle.
In this study, we generated global gene expression profiles for mature fungal mat and zoospore-enriched samples for the model anaerobic fungus Neocallimastix californiae G1, a monocentric, rhizoidal, polyflagellate strain that has previously been genomically sequenced (Haitjema et al., 2017) and subject to other bulk omics characterization (Solomon et al., 2016;Swift et al., 2021). To generate global gene expression profiles, we extracted RNA from two sample types of cultured N. californiae: (1) cell pellets enriched in free-swimming zoospores with some germlings present and (2) fungal mats with mixed life stages. We grew cultures in serum bottles with a rumen fluid-based medium formulation, Medium C (Theodorou and Brookman, 2005), with simple, soluble sugar glucose as the substrate. Although both sample types contained some mixture of fungal cell types, the cell pellets were enriched for zoospores by a Miracloth filtration method, while the extracted mat samples, which consist of thick rhizoidal biofilms, included mature sporangia.
This study marks the first endeavor to transcriptomically compare the two most extreme life stages of anaerobic gut fungi to better characterize their basic cell biology, as well as to direct future efforts towards enzyme and metabolite production and ease of genetic manipulation. We find that the fungal zoospore-enriched cells, which were only exposed to glucose, nevertheless primed to encounter more complex substrates. Chiefly, catabolic CAZymes and putative carbohydrate transporters are upregulated in the zoospore-enriched samples. Additionally, differential regulation of transcription factors and putative secondary metabolite genes across life stages inform their physiological roles. The transcriptomic dataset presented here serves to inform hypotheses for future experimentation with anaerobic gut fungi, with applications in understanding their fundamental biology and harnessing their biotechnological potential.
The N. californiae G1 strain was previously isolated via reed canary grass enrichment from goat feces from the Santa Barbara Zoo (Haitjema et al., 2017;Solomon et al., 2016a). Recent evidence supports phylogenetic synonymization of N. californiae, Neocallimastix lanati, and Neocallimastix cameroonii (Stabel et al., 2021). Fungal cultures were incubated at 39 • C and routinely passaged every 3-4 days by adding 1.0 mL of fungal inoculum to a Hungate tube with 9 mL of fresh, autoclaved, anaerobic, reduced formulation of Medium C (M. K. Theodorou, J. Brookman, 2005), which contains 0.1 g of 4 mm milled reed canary grass as growth substrate, 0.25 g of yeast extract (Thermo Fisher Scientific), 0.5 g of Bacto™ Casitone, and 7.5 % clarified rumen fluid. Growth rates for N. californiae were approximated by changes in culture headspace pressure due to accumulation of fermentation gases, and typical growth rates have been previously reported for this growth condition (Solomon et al., 2016b). As necessary, cultures were preserved and revived cryogenically using established methods (Solomon et al., 2016).
Fig. 1. Anaerobic gut fungi have a unique life cycle. (A) In the anaerobic gut fungal life cycle, motile zoospores encyst into substrate and develop into sporangia, which can grow rhizoids, reproduce asexually, and release new zoospores upon maturation (Gruninger et al., 2014;Heath et al., 1986). The understudied, potentially aerotolerant life stage remains undescribed in this study. (B) A zoospore pellet sample and (C) an extracted fungal mat sample for N. californiae G1 were imaged with 20× magnification on a light microscope. Figure was created with BioRender.com.
All cultures were prepared anaerobically in 100 mL serum bottles with 80 mL autoclaved Medium C and 4 mL of 0.1 g/mL sterile-filtered glucose solution, yielding 84 mL Medium C cultures with 0.48 % (w/v) glucose (M. K. Theodorou, J. Brookman, 2005). All inoculations were performed with 1.0 mL of N. californiae inoculum from 2-day seed cultures. Per zoospore-enriched RNA replicate, a cell pellet was generated by consolidating and filtering batches of 16 cultures (or 15 cultures for zoospore replicate HHCCZ, due to one culture not growing). Batches of cultures were inoculated from two seed cultures, incubated at 39 • C, and harvested after 2 days of growth.
To harvest samples, all 85 mL cultures in a batch were filtered through autoclaved Miracloth (Millipore Sigma) set in a funnel over a media flask, with all parts sprayed with RNase AWAY (Thermo Scientific). Per batch, forceps were used to extract several whole fungal mats as cultures were poured over the Miracloth. These mats were placed in separate 15 mL Falcon tubes with 1 mL of RNAlater, flash-frozen, and stored at -80 • C until RNA extraction. The flow-through with zoospores was aliquoted into multiple 50 mL Falcon tubes that were then successively centrifuged at 2000 xg for 10 min with a swinging bucket rotor at 4 • C. After each round of centrifugation, the supernatant was decanted and pooled. Cell pellets were each resuspended in spent media and consolidated in a single Falcon tube. The process of centrifugation and consolidation was repeated until all pellets were combined, after which the final cell pellet was resuspended in 1 mL of spent media, transferred to a 1.5 mL Eppendorf tube, and centrifuged in a microcentrifuge at 20,000 x g for 5 min at 4 • C. The remaining supernatant was removed by pipetting, 0.5 mL RNAlater was pipetted over the pellet, and the sample was flash-frozen and stored at -80 • C until RNA extraction. For future application of the described method, we recommend filtering through several layers of Miracloth, introducing RNAlater to the zoosporeenriched sample as early as possible to minimize RNA degradation, following all RNAlater guidelines, and inspecting each zoosporeenriched cell pellet under the microscope for proper filtration.
To approximate zoospore yields, the harvest method was performed in smaller batches of cultures (3 × 85 mL cultures per batch) and without RNA-preserving measures (RNase AWAY and RNAlater). Harvested zoospore pellets were resuspended to an amount yielding 100-200 cells per grid on an Improved Neubauer Bright-Line Hemacytometer (Hausser Scientific) to count and classify cells.
Total RNA was extracted and purified from harvested fungal mat and zoospore-enriched samples using the Qiagen RNeasy Mini Kit (Qiagen, Valencia, CA), via the plant and filamentous fungi protocol with liquid nitrogen grinding cell lysis and on-column DNase Digest, as previously described for extracting RNA from anaerobic gut fungi (Solomon et al., 2016). For each sample, RNA quality was measured by RINe (RNA Integrity Number equivalent) score using an Agilent TapeStation 2200 (Agilent, Santa Clara, CA), and total RNA quantity was measured using Qubit Fluorometric Quantitation (Qubit, New York, NY) (data not shown). Samples were then transferred to the Joint Genome Institute and underwent a second quality control check for RQN (RNA Quality Number) score using an (Supplementary Table 1). Ultimately, RNA samples from three zoospore-enriched pellets and 15 fungal mats possessed sufficient quality for sequencing. The discrepancy in the number of generated replicates (3 vs. 15) for zoospore-enriched pellets vs. fungal mats is due to the nature of combining groups of cultures to acquire a single zoospore-enriched replicate, whereas each individual culture can generate one mat replicate. A few fungal mats were arbitrarily selected from each attempted culture batch to check for mat variation between batches.
RNA-Seq was performed at the Joint Genome Institute (JGI) using a 96-Capillary Fragment Analyzer (Agilent), and stranded RNA-Seq libraries were created and quantified by qPCR. Plate-based RNA sample prep was performed on the PerkinElmer Sciclone NGS robotic liquid handling system using the Illumina TruSeq Stranded mRNA HT sample prep kit utilizing poly-A selection of mRNA. Amplification libraries were prepared by amplifying 100 ng of total RNA per sample over 10 cycles of PCR. The prepared libraries were then quantified using the KAPA Illumina library quantification kit (Roche) and run on a LightCycler 480 real-time PCR instrument (Roche). The quantified libraries were then multiplexed, and the pool of libraries was prepared for sequencing on the Illumina NovaSeq 6000 sequencing platform using NovaSeq XP v1.5 reagent kits (Illumina), S4 flow cell, following a 2 × 150 indexed run recipe.
Raw fastq files were filtered and trimmed using the JGI quality control pipeline. Using BBDuk (https://sourceforge.net/projec ts/bbmap/), raw reads were evaluated for artifact sequence by kmer matching (kmer = 25), allowing 1 mismatch, and detected artifacts were trimmed from the 3′ end of the reads. RNA spike-in reads, PhiX reads, and reads containing any Ns were removed. Quality trimming was performed using the phred trimming method set at Q6. Finally, following trimming, reads under the length threshold were removed (minimum length 25 bases or 1/3 of the original read length -whichever was longer).
Filtered reads from each library were aligned to the reference genome (Haitjema et al., 2017) using HISAT2 version 2.2.0 (Kim et al., 2015). The average genome mapping was 93.67 % with 14-48 M genome-mapped reads per library. The dataset suffered high rRNA contamination, averaging rRNA contamination in 14.14 % of raw reads. Strand-specific coverage bigWig files were generated using deepTools v3.1 (Ramírez et al., 2014). The tool featureCounts (Liao et al., 2014) was used to generate the raw gene counts (Supplementary Material 2) file using gff3 annotations. Only primary hits assigned to the reverse strand were included in the raw gene counts. Pearson's correlation was calculated for raw gene counts to evaluate the level of correlation between biological replicates, and ultimately, 15 mat samples and two zoospore-enriched samples were selected for subsequent analysis, with the zoospore-enriched replicate HHCCZ being removed for suspected mat contamination (Supplementary Fig. 1). All gene counts, transcript per million (TPM) normalized gene counts, mapping statistics, and strandedness estimations are provided in Supplementary Material 2.
Prior to performing differential gene expression analysis, raw counts for transcripts encoding identical amino acid sequence were combined under single proteinIDs, with the rationale that protein-encoding genes for identical proteins represent one functional annotation. This consolidation step was also performed for the TPM counts, and the schema for how proteinIDs were consolidated is described in Supplementary Material 2. DESeq2 (v1.36.0) (Love et al., 2014) was used to determine which genes were differentially expressed between mat and zoosporeenriched samples.
Gene annotations were aligned to each N. californiae gene in Supplementary Material 3. Annotations from the Joint Genome Institute (JGI) MycoCosm portal (Grigoriev et al., 2014) were predominantly used to annotate genes detected in the transcriptomic data. The following gene annotation types were considered: KOG (Eukaryotic Orthologous Groups of proteins) (Koonin et al., 2004), GO (Gene Ontology) (Ashburner et al., 2000;Carbon et al., 2021), KEGG (Kyoto Encyclopedia of Genes and Genomes) (Kanehisa, 2019;Kanehisa et al., 2023;Kanehisa and Goto, 2000), and IPR (InterPro) (Paysan-Lafosse L.V. Butkovich et al. Fungal Genetics and Biology 176 103958 et al., 2023). Additional, more specific annotations were also considered from MycoCosm: 1450 carbohydrate-active enzymes (CAZymes) and 39 putative secondary metabolite biosynthetic genes, generated using a Secondary Metabolite Unknown Regions Finder (SMURF) algorithm (Grigoriev et al., 2014;Khaldi et al., 2010;Swift et al., 2021). The dbCAN2 tool (Zhang et al., 2018) was used to putatively assign CAZyme functions to 14 additional genes, based on the following requirements:
(1) at least two out of three dbCAN2 sub-tools classify the gene as a CAZyme and ( 2) the gene has a KOG or InterPro functional annotation. 418 genes with dockerin domains, DOC2 or carbohydrate-binding domain family 10 (CBM10), were annotated. 37 scaffoldin genes were annotated by local BLASTp comparison (e-value <0.01) to previous N. californiae scaffoldin annotations (Haitjema et al., 2017). 9 Sugars
Will Eventually be Exported Transporters (SWEETs) were annotated based on previous work with N. californiae SWEETs (Podolsky et al., 2021). Instead of mitochondria for energy production, anaerobic gut fungi possess hydrogenosomes. Putative hydrogenosome functions were assigned to 69 N. californiae genes based on their homology to hydrogenosome-associated genes previously identified for Caecomyces churrovis (Leggieri et al., 2022). DeepLoc (Almagro Armenteros et al., 2017) was used to predict localization of transcripts to subcellular locations: the cytosol, mitochondria, plastid, or other cellular compartments (lenient confidence cutoff of 50 %). 14 out of 67 genes with putative hydrogenosome function were DeepLoc-localized to the hydrogenosome, based on localization to mitochondria or plastid. Ten pyruvate formate lyases (PFLs) were localized to the cytoplasm, and ten other PFLs were localized to the hydrogenosome. 172 genes were annotated as putative transcription factors based on the following criteria:
(1) GO ID 3700 ("DNA-binding transcription factor activity") or PF00096 ("Zinc finger, C2H2 type") and ( 2) KOG annotation (KOG defline or KOG class) containing "transcription" or "Zn-finger". Velvet regulatory proteins are transcription factors shown to coordinate development and secondary metabolism in filamentous fungi, where they can form complexes with each other and the methyltransferase LaeA (Chen et al., 2024;Hou et al., 2024). By BLASTp comparison (evalue <1E-5), 10 N. californiae velvet proteins were annotated based on homology to A. nidulans velvet proteins (VosA, VeA, VelB, VelC), and no N. californiae LaeA was annotated, based on homology to A. nidulans
LaeA. Lastly, we supply results from OrthoFinder (Emms and Kelly, 2019) in Supplementary Material 3, to describe sequence-similar genes within the N. californiae genome, as well as gene conservation across three other anaerobic gut fungal species (Anaeromyces robustus, Caecomyces churrovis, Piromyces finnis), enabling predictive, homology-based targeting of genes of interest in these species.
We define significant differential expression of individual genes by the following three requirements: (1) q < 0.05, where q denotes adjusted p-values from DESeq2, (2) |log2 fold-change| > 1 for comparing DESeq2-normalized expression levels, and (3) average TPM count greater than 1 for the upregulated condition (Supplementary Fig. 2 and Supplementary Material 3). While there is no universally accepted TPM cutoff for defining gene expression, we include a cutoff here to omit consideration of extremely low abundance transcripts from the differentially expressed gene sets. Negative log2 fold-change refers to higher DESeq2-normalized counts in zoospore-enriched samples, and positive log2 fold-change refers to higher DESeq2-normalized counts in mats.
We define significant differential expression of gene annotations by (1) gene count greater than 10 and (2) p < 0.05, where p denotes pvalues from Fisher's Exact Tests, described below (Supplementary Material 4). For any given KOG, GO, KEGG, InterPro, or CAZyme annotation, a Python script determined total gene count, number of genes upregulated in zoospore-enriched samples, and number of genes upregulated in mats (Supplementary Material 4). These data were used to run statistical Fisher's Exact Tests, summarized in Supplementary Fig. 3, and define significant differential expression of gene annotations (gene count greater than 10 and p < 0.05). The gene count cutoff was implemented because Fisher's Exact Test applied to small sample sizes can yield misleading statistics and is unlikely to be representative of the population of genes under consideration. Supplementary Table 2 lists broad KOG classes that are differentially regulated between mat and zoospore-enriched samples by Fisher's Exact Test statistics, although to glean biological meaning behind differential regulation, it is advised to inspect specific functions, rather than broad categories.
To generate expression profiles for different life stages of anaerobic fungi, we developed a method to separate the two most extreme life stages and extracted high quality RNA from each. The workflow we designed to generate these enrichments targeted the zoospores and fungal mats (Fig. 2). Calkins et al. summarizes previously used collection methods for anaerobic gut fungal zoospores (Calkins et al., 2016), and we implemented culture filtering (typical Miracloth pore size of 22-25 μm, N. californiae zoospore size 10 μm) for its simplicity and suitability for nucleic acid extractions (Tsai and Calza, 1992). This method is conducive to processing large culture volumes (1 L<) to obtain a sizable, zoospore-enriched cell pellet for RNA extraction using liquid nitrogen grinding for cell lysis. We selected a culture length of two days, because we expected RNA quality to decrease with culture age as noted for bulk gut fungal samples (Brown et al., 2023). Notably, after performing RNA-Seq with harvested samples, one zoospore-enriched pellet replicate (HHCCZ) was removed from analysis due to suspected mat contamination (see Materials and Methods) (Supplementary Fig. 1). For overall method validation, zoospore yields were approximated using a hemacytometer to count and classify cells (see Materials and Methods). Two cell types were discernible in cell pellet samples: zoospores and germlings. The concentrations per cell type were 713 ± 64 zoospores/mL and 164 ± 33 germlings/mL (n = 3), and therefore the expected content of zoospore-enriched pellets is ~80 % zoospores. The yield of flowthrough cells was relatively low compared to previous methods but sufficient for RNA extraction and sequencing. RNA extracted from the harvested samples possessed sufficient RNA quality for RNA-Seq, with most RNA quality numbers greater than 8 (Supplementary Table 1).
Thousands of genes were differentially expressed between enriched fungal zoospores and mats (Fig. 3). Analyzing the differentially expressed general and specific annotations informs biological hypotheses. We analyzed detected transcripts that mapped to the genome, totaling 19,968 N. californiae genes with unique amino acid sequences, with most of these genes possessing at least one annotation (see Materials and Methods). Per sample group, we defined expressed genes as having average transcripts per million (TPM) counts >1, and we defined genes that were not expressed as having average TPM ≤ 1 (Supplementary Fig. 2). 12,170 genes (60.9 %) were expressed in mats, with 3828 of those genes not expressed in zoospore-enriched samples. 8978 genes (45.0 %) were expressed in zoospore-enriched samples, with 636 of those genes not expressed in mats. Additionally, 2693 (13.5 %) genes were upregulated in mats, and 3386 (17.0 %) genes were upregulated in zoospore-enriched samples. Taken altogether, mats expressed a larger range of genes, but zoospore-enriched samples specifically upregulated more genes. Of the 4668 (23.4 %) unannotated genes, which lack a GO, KOG, InterPro, or KEGG annotation, 639 genes were expressed only in mats, and 155 genes were expressed only in zoospore-enriched samples. These subsets of unannotated genes likely play life stage-specific roles and are interesting candidates for targeted gene studies.
In this study, we note that zoospore-enriched samples upregulated "glutathione metabolism" (KEGG pathway, p = 1.3E-9) for the oxidative stress response, possibly due to the nature of the harvesting method which briefly exposed these samples to oxygen (Supplementary Fig. 4).
L.V. Butkovich et al. Fungal Genetics and Biology 176 103958 Glutathione (Glu-Cys-Gly in reduced form) is a major antioxidant molecule expected to help protect anaerobic gut fungi from oxidative damage. When zoospore-enriched pellets and fungal mats were harvested, the cells were taken from an anaerobic environment and exposed to an aerobic environment. Relative to the fungal mats that were treated immediately with RNAlater, zoospore pellet samples were more exposed to ambient oxygen prior to RNAlater treatment. As a result, it is possible that oxygen exposure affected the differential regulation between the zoospore and mat sample groups. For future use of this method, we recommend immediate RNAlater treatment of the zoospore-enriched samples (5-10× RNAlater volume relative to the culture flowthrough) to address the effect of oxygen exposure. Transcriptomic profiling of anaerobic gut fungal zoospores enables comparison to studied zoospores of other phyla. Historically, all zoosporic fungi were placed into one early-diverging fungal phylum (Barr, 2001), but now they are divided into three phyla: Neocallimastigomycota, Blastocladiomycota (blastoclads), and Chytridiomycota (chytrids). Anaerobic gut fungi differ in multiple ways from the other zoosporic fungi. While anaerobic gut fungal zoospores are anaerobic, are host symbionts, possess hydrogenosomes, and potentially possess a cell wall based on calcofluor white staining (Lillington et al., 2021), chytrid and blastoclad zoospores are aerobic, are saprobes or parasites, possess mitochondria, and lack a cell wall (Gleason and Lilje, 2009;Gruninger et al., 2014;Li et al., 1993;Money, 2016). Previous studies showed chytrid zoospores of Blastocladiella emersonii are metabolically active (not dormant), and do not produce new DNA, RNA, or proteins until after germination, with dormant ribosomes containing maternally-derived RNA (Barstow and Lovett, 1974;Laundon et al., 2022;Silva et al., 1987). Anaerobic gut fungal zoospores are similarly metabolically active and do not produce DNA, but comparable experiments for transcriptional and translational activity have not been performed for anaerobic gut fungal zoospores. Therefore, the zoospore RNA sequenced in this study was potentially maternally-derived. Fig. 3. Overall statistics compare differential regulation of N. californiae genes between zoospores and mats. A) This plot depicts the number of unique genes with various annotation types and levels of expression in zoospore and mat samples (see Materials and Methods). *"None" refers to the number of genes that do not possess a KOG, GO, IPR, or KEGG annotation. The "Levels of Expression" describes how many genes are transcribed with average transcripts-per-million (TPM) counts >1 in a given grouping. **"Not detected" refers to the number of genes that had raw count = 0 in all mat and zoospore samples. B) More genes were specifically upregulated in zoospores. The volcano plot shows the spread of differentially expressed genes between mat and zoospore samples, with each point representing one gene. Significantly regulated genes satisfy the following cutoffs: q < 0.05, |log2 fold-change| > 1, and average TPM > 1 for the upregulated condition. C) This plot depicts differential regulation of genes between zoospores and mats. "Transcribed" refers to genes with average TPM > 1 across replicates in either the mat or zoospore group. Plots in (A) and (C) were generated with PRISM.
L.V. Butkovich et al. Fungal Genetics and Biology 176 103958 In comparing the zoospore and mature life stages, we expect differential gene expression related to known differences in physiology. Zoospores do not carry out mature functions but are well-prepared to colonize new substrate. Firstly, to describe genes implicated in flagellar motility for zoospores, we considered annotations for cell motility. None of the genes in the KOG class for cell motilitysuch as myosins, light chain dynein, and dynactinwere upregulated in zoospore, and some were upregulated in mats. These genes are either flagella-associated and upregulated in developing zoospores within sporangia, or they have other cellular functions, such as cell division and intracellular transport. Secondly, zoosporogenesis enables reproduction in anaerobic gut fungi, but the genes involved in its onset and progression are not well characterized. For the future study of zoosporogenesis, such genes are likely in the gene set upregulated in mats and not expressed in zoosporeenriched samples. Thirdly, mats upregulated "actin cytoskeleton organization and biogenesis" (GO name, p = 1.5E-3) and "dyneins, heavy chain" (KOG defline, 3.9E-5), with putative functions in the growth, maintenance, and intracellular transport of rhizoid networks. Mats also upregulated the "proteasome endopeptidase complex" (KEGG definition, EC number 3.4.25.1, p = 7.7E-18) for degrading intracellular proteins, while zoospore-enriched samples upregulated ribosomerelevant genes (GO name, p = 3.3E-27). Laundon et al. previously showed that chytrid zoospores of Rhizoclosmatium globosum also upregulate ribosome-relevant genes relative to more mature cells (germlings) (Laundon et al., 2022). These results suggest that mature cells prioritize quality control of existing proteins, while young cells prioritize production of new proteins.
Anaerobic gut fungi possess a variety of regulatory mechanisms to control gene expression, resulting in the observations discussed in this study. One such regulatory mechanism is transcription factors, sequence-specific DNA-binding proteins present in all organisms that control gene expression, via activation or repression, and play a role in cell development. The system of transcription factors in anaerobic gut fungi is not well described, and improved understanding would enable better laboratory manipulation of these organisms. Here, we putatively annotated 172 N. californiae genes as transcription factors (see Materials and Methods), and transcripts for 160 of these were detected in any of the samples. Transcription factors were upregulated in zoosporeenriched samples (p = 4.3E-10), with 60 upregulated in zoosporeenriched samples and 14 upregulated in mats (Supplementary Fig. 5). We highlight the differentially regulated transcription factors as targets for further study. Overall, the zoospore-enriched samples appear to exhibit active and dynamic transcriptional regulation via a range of transcription factors.
Of any sequenced fungi, anaerobic gut fungi possess the most genes for carbohydrate-active enzymes (CAZymes), which contribute to the breakdown of lignocellulosic plant biomass to release sugars and metabolically synthesize fatty acids (Gilmore et al., 2020;Hooker et al., 2019;Solomon et al., 2016a). CAZymes either freely diffuse in the extracellular environment or co-localize to cellulosomes (Haitjema et al., 2017), which are multienzyme complexes that were first described in anaerobic bacteria and function to enhance degradation of recalcitrant plant matter (Artzi et al., 2017). Regulation of CAZymes and cellulosome components has not been transcriptionally studied in separated life stages of anaerobic gut fungi. Previous imaging probe studies by Lillington et al. revealed that anaerobic gut fungal zoospores (Piromyces finnis and N. californiae) displayed cellulosome-associated proteinswith exoglucanase glycoside hydrolase family 48 (GH48) CAZyme domains and dockerin structural domainswhen grown on different substrates, including glucose. In contrast, sporangia only displayed cellulosome-associated proteins on more difficult substrates (Lillington et al., 2021). Additionally, previous transcriptomic work showed CAZymes are broadly catabolite repressed in glucosesupplemented cultures of anaerobic gut fungi with mixed life stages present (Henske et al., 2018a). Even so, previous studies also lead us to expect some detection of CAZyme transcripts from gut fungal cultures with only soluble sugar substrate, likely due to high basal expression levels of these enzymes (Wilken St. et al., 2021;Williams and Orpin, 1987).
Here, we describe differential expression of CAZymes between N. californiae mat and zoospore-enriched samples cultured with glucose as the substrate. Broad CAZyme types were differentially expressed between mat and zoospore-enriched samples (Table 1). Glycoside hydrolases, carbohydrate esterases, and polysaccharide lyases are three catabolic CAZyme types responsible for breakdown of plant polysaccharides, and all three CAZyme types were significantly upregulated in zoospore-enriched samples (Fig. 4A). This finding suggests that zoospores prime for introduction to diverse food sources, regardless of the substrate complexity in the immediate environment or consumed by the parent sporangia. Carbohydrate-binding modules and glycosyltransferases were also significantly upregulated in mats. A subset of CAZymes target endogenous chitin, rather than exogenous substrate, to modulate cell walls and develop rhizoids. For example, carbohydratebinding module family 18 (CBM18), a subset of the carbohydratebinding modules upregulated in mats (Table 1), possesses chitinbinding function (Liu and Stajich, 2015) and was specifically upregulated in mats (p = 8.2E-6) (Supplementary Fig. 6 A). Such matupregulated CAZymes can inform gene targets for future experiments to better understand rhizoid development in mature sporangia.
The catabolic CAZyme types (carbohydrate esterases, glycoside hydrolases, and polysaccharide lyases) were upregulated in zoosporeenriched samples, and carbohydrate-binding modules and glycosyltransferases were upregulated in mats. *Some CAZyme types have pvalues less than 0.05 in both the test for upregulation in zoosporeenriched samples and the test for upregulation in mats. This indicates a subset of genes preferentially expressed in the zoospore-enriched samples and a separate subset are preferentially expressed in the mat samples. In these cases, we make a direct comparison of the p-values to define significant differential expression of an annotation between the sample groups.
In addition to considering regulation of broad CAZyme types, we considered specific CAZyme families to better describe metabolic strategies across life stages. One observation is that glycoside hydrolase families 1 and 3 (GH1/GH3), which contain beta-glucosidases that hydrolyze cellobiose to glucose during cellulose breakdown (Drula et al., 2022), were upregulated in zoospore-enriched samples (GH1: p = 6.2E-6, GH3: p = 1.5E-5) (Supplementary Fig. 6B -C). This observation supports the understanding that zoospores seek out new substrate to colonize and are primed to encounter cellobiose, a key product in plant biomass breakdown in the rumen. A second observation is that the glycosyltransferase family 8 (GT8), which contains glycogenin that stores glucose in the polysaccharide glycogen (Wilson et al., 2010;Zeitz et al., 2019), was likely upregulated in mats, with 6/10 GT8 genes upregulated in mats (Supplementary Fig. 6D). Previous electron microscopy showed that glycogen occupies a large portion of the cytoplasm of Neocallimastix frontalis zoospores (Munn et al., 1981). Here, the transcriptomic analysis supports the hypothesis that mature sporangia assist in preparing the glycogen reservoirs in developing zoospores prior to their release. Additionally, overlaying gene regulation and the metabolic mapping of chitin modulation suggests potential upregulation of glycogen utilization in mats (Supplementary Fig. 7).
Cellulosomes play a major role in how anaerobic gut fungi metabolize plant biomass. Here, we find that cellulosome-associated proteins did not generally differentially express between the mat and zoosporeenriched samples. To approximate cellulosome activity, we considered differential expression of known cellulosome components (Fig. 4B-E): (1) proteins with scaffoldin and dockerin structural domains (Raghothama et al., 2001), and (2) GH48s, which are exoglucanase L.V. Butkovich et al. Fungal Genetics and Biology 176 103958 CAZymes highly expressed in bacterial cellulosomes (You et al., 2023).
Cellulosomes can tether a variety of CAZymes, but in this analysis, we focus on the GH48 CAZyme. Cohesin domains are also known structural components of cellulosomes, but cohesin sequences have not been determined in anaerobic gut fungi. Here, 37 genes with scaffoldin domains were upregulated in mats (p = 7E-4) (Fig. 4C), and 418 genes with dockerin domains were not differentially expressed between mat and zoospore-enriched samples (Fig. 4D). Overall, the GH48 CAZyme annotation was not upregulated in mats (p = 0.75) or zoospore-enriched samples (p = 0.13), but some individual GH48 genes were highly expressed in either group (Fig. 4E). Overall, the results support the understanding that cellulosomes play a role in young and mature cells and that zoospores produce cellulosomes independent of growth substrate complexity. Since previous results showed simple sugars repress cellulosome production in mature cells (Lillington et al., 2021), the expression of cellulosome-associated proteins and upregulation of scaffoldincontaining proteins in mats suggests catabolite repression of cellulosome production occurs post-transcriptionally.
To utilize the breakdown products from plant biomass, anaerobic gut fungi also possess an array of carbohydrate transporters. Some of these are denoted as Sugars Will Eventually be Exported Transporters (SWEETs) (Podolsky et al., 2021;Seppälä et al., 2016), a superfamily of sugar transporters that are abundant in plant genomes and have a wide range of affinities to mono-and disaccharide sugars. Here, 5/9 SWEETs D-E) Overall, dockerin-containing proteins and GH48s, potential cellulosome-associated CAZymes, were not differentially regulated between mat and zoospore-enriched samples. We highlight specific proteins with these annotations that were differentially expressed. One gene (proteinID 703,870), which encodes a protein with two dockerin domains and a GH48 domain (GH48-DOC2-DOC2), was the fourth most expressed gene in zoospores (by average TPM count) and was upregulated in zoospores (q = 2.0E-4). By these statistics, proteinID 703,870 appears to be an important cellulosome-associated GH48 protein, especially early in the anaerobic gut fungal life cycle. For each volcano plot, differentially regulated genes with relatively high average TPM count are displayed by proteinID name. For genes upregulated in zoospores, data point radius scales with average TPM count in zoospores, and for genes upregulated in mats or not differentially regulated, data point radius scales with average TPM count in mats. GH = Glycoside Hydrolase. Figure was created with BioRender.com.
L.V. Butkovich et al. Fungal Genetics and Biology 176 103958 were upregulated in zoospore-enriched samples, but all of the SWEETs were lowly expressed in both sample groups (average TPM count less than 4) (Supplementary Table 3). Interestingly, other transporter annotations from MycoCosm were significantly upregulated in zoosporeenriched samples compared to the mats. These annotations include "predicted transporter (major facilitator superfamily)" (KOG defline, p = 3.6E-18), "carbohydrate transport" (GO name, p = 3.3E-16), and "sugar transporter" (InterPro description, p = 8.0E-20). Upon closer inspection, the genes with these three annotations largely overlap with each other and with the "L-arabinose isomerase" annotation, which was also strongly upregulated in zoospore-enriched samples (KEGG definition, p = 5.3E-17) (Supplementary Fig. 8). These genes are also wellconserved in anaerobic gut fungi, with 56/89 "L-arabinose isomerase" genes sharing a set of homologs in three other anaerobic gut fungal species (see Materials and Methods). L-arabinose isomerase is not a transporter, raising questions about why these genes are annotated with both transporter and L-arabinose isomerase functions. Transporter mechanisms are not yet well understood in anaerobic gut fungi, likely contributing to the confounding annotations for this set of genes, which is clearly upregulated in the zoospore-enriched samples. We hypothesize that these genes upregulated in the zoospore-enriched samples are putative carbohydrate transporters and help zoospores seek out substrate to colonize, but these gene annotations require further validation.
Primary metabolism in anaerobic gut fungi consists of glycolysis, an incomplete citric acid cycle, and processes in the poorly characterized fungal hydrogenosome (Boxma et al., 2004;Marvin-Sikkema et al., 1994;Wilken St. et al., 2021). Anaerobic gut fungi generate most of their energy through glycolysis via anaerobic fermentation, and they possess hydrogenosomes, instead of mitochondria, for some energy production with hydrogen as a notable by-product. Previously, hydrogenase activity was detected in cell-free extracts of zoospores and vegetative growth of anaerobic gut fungi, and hydrogenosomes were detected with electron microscopy, suggesting that hydrogenosomes are present across the life stages (Yarlett et al., 1986). However, the transcriptional regulation of hydrogenosome-associated genes across life stages has not been previously described. The metabolic mapping of hydrogenosomes remains incomplete, but approximate representations have been described (Leggieri et al., 2022;Wilken St. et al., 2021). We performed differential gene expression analysis on primary metabolic enzymes between mat and zoospore-enriched samples to inform hypotheses for preferred Fig. 5. Some primary metabolic pathways were potentially differentially regulated between anaerobic fungal mat and zoospore-enriched samples, based on the regulation of relevant enzymes. This schematic depicts central metabolism (glycolysis, hydrogenosome, citric acid cycle) in N. californiae. Some genes in the pathway for xylose consumption were upregulated in mats (a-c), and 10/10 cytosol-localized PFLs were upregulated in zoospores (k). The following enzyme key describes the Enzyme Commission (EC) number, if available, and differential expression statistics for enzymes involved in primary metabolism. Expression statistics are reported for each enzyme as (number of N. californiae genes with annotations: number of annotated genes upregulated in mats: number of annotated genes upregulated in zoospores). The enzyme annotations presented here lack Fisher's Exact Statistics, due to the small gene count sizes per annotation (see Materials and Methods). Enzymes [with EC numbers if available]: a) xylose isomerase [5.3.1.5], (7:5:1). No genes were annotated with EC number 5.3.1.5, so the InterPro description for xylose isomerase was used instead. Two of the xylose isomerase genes were highly expressed in the mats (average TPM > 1000) with log2 fold-change >3. b) xylulokinase [2.7.1.17], (2:2:0). c) glycoaldehydetransferase [2.2.1.1], (3:1:0). d) triosephosphate isomerase [5.3.1.1], (2:0:0). e) glyceraldehyde 3-phosphate dehydrogenase [1.2.1.12], (7:1:1). Most of these genes were highly expressed in both matsand zoospore-enriched samples (average TPM > 100). f) phosophoglycerate kinase [2.7.2.3], (4:1:0). g) glyceraldehyde 3-phosphate dehydrogenase (NADP+) [1.2.1.9], (1:1:0). h) phosphopyruvate hydratase [4.2.1.11], (2:0:0). i) phosphoenolpyruvate synthase [2.7.9.2], (8:0:3). L.V. Butkovich et al. Fungal Genetics and Biology 176 103958 metabolic pathways of the life stages. Some primary metabolic pathways had differentially regulated genes between the mat and zoospore-enriched samples (Fig. 5). Anaerobic gut fungi typically possess enzymes and transporters to consume xylose, a pentose sugar and main constituent of the hemicellulose xylan in plant biomass (Henske et al., 2018b;Wilken St. et al., 2021). Enzymes for xylose processing, xylose isomerases and xylulokinases, had more genes upregulated in mats (Fig. 5). These observations indicate that mature life stages may be better equipped to process xylan from plant biomass breakdown in the rumen. Some steps for fumarate production also had genes upregulated in mats (Fig. 5). Of the enzymes involved in primary metabolism, pyruvate formate lyases (PFLs) were the most differentially regulated and were upregulated in the zoospore-enriched samples (InterPro description, p = 2.9E-13). PFLs convert pyruvate to acetyl-coA and formate in the cytoplasm and hydrogenosome. Considering only PFLs that were localized to the cytoplasm (see Materials and Methods), we find ten cytosol-localized PFL genes that were all upregulated in the zoospore-enriched samples (Fig. 5). Overall, the transcriptional regulation of primary metabolic pathways can inform the life stage-dependent metabolic strategies in anaerobic gut fungi.
Regarding the hydrogenosome, 14 genes with detected transcripts were putatively annotated and localized to the hydrogenosome (see Materials and Methods) (Supplementary Material 3). Pyruvate formate lyase, acetyl coA-hydrolase, and succinyl coA synthetases can be involved in the ATP synthesis pathway in gut fungal hydrogenosomes, therefore their expression may correlate with hydrogenosome activity (Wilken St. et al., 2021). Similar in expression behavior to the cytosollocalized PFLs, 9/10 hydrogenosome-localized PFLs were upregulated in the zoospore-enriched samples, and 1/10 was upregulated in mats. 2/ 2 acetyl coA-hydrolases were upregulated in the zoospore-enriched samples, but neither encoded protein was predicted to localize to the hydrogenosome. Succinyl coA (SucCoA) synthetases are the putative payoff enzyme in our current understanding of the gut fungal hydrogenosome. Two genes are putatively annotated "SucCoA A" and three are putatively "SucCoA B," but none of these were differentially expressed between mat and zoospore-enriched samples or localized to the hydrogenosome. Overall, the lack of differential expression of the putative payoff enzymes suggests that hydrogenosome activity is not differentially regulated between the mat and zoospore-enriched samples.
Recent work described anaerobic gut fungi as under-explored for secondary metabolites, compounds that often possess bioactivities, like antibiotic activity, and improve survivability of an organism under specific circumstances (Swift et al., 2021). To date, no nonribosomal peptide or polyketide secondary metabolites have been functionally characterized from Neocallimastigomycota, and their native roles are unknown. In higher order fungi, some secondary metabolites control or coincide with fungal development, so we expect some secondary metabolites in anaerobic gut fungi play life stage-specific roles (Gerke and Braus, 2014). For example, the UV-protectant pigment melanin is required for some spore survival and is biosynthesized by either a polyketide pathway or a tyrosinase pathway (Calvo et al., 2002;Cordero and Casadevall, 2017;Eisenman and Casadevall, 2012). In the N. californiae genome, the JGI MycoCosm predicts 39 biosynthetic gene clusters for nonribosomal peptides or polyketides (see Materials and Methods) (Grigoriev et al., 2014;Swift et al., 2021). A biosynthetic gene cluster, as depicted in Fig. 6A, includes a core gene that encodes an enzyme to produce the first biosynthetic intermediate for a secondary metabolite. Typically, fungi co-regulate genes in a biosynthetic gene cluster via complex regulatory systems (Zhgun, 2023). In this study, some putative secondary metabolite core genes were differentially expressed between the mat and zoospore-enriched samples (Fig. 6B), although some predicted accessory genes did not co-differentially express with their corresponding core genes (Supplementary Material 3), possibly due to incorrect annotation of some accessory genes. In filamentous fungi, complexes of velvet regulatory proteins and the methyltransferase LaeA coordinate secondary metabolism and development, and in most known cases, these complexes positively regulate target secondary metabolite genes (Chen et al., 2024;Gerke and Braus, 2014;Hou et al., 2024). Of 10 N. californiae genes annotated as putative velvet proteins (see Materials and Methods), 3 genes (proteinIDs 514,349,514,348,223,018) were upregulated in the zoospore-enriched samples, tentatively suggesting their regulation of the putative secondary metabolite genes upregulated in these samples (proteinIDs 701,295,407,062,704,767). No LaeA homologs were identified in N. californiae (see Materials and Methods), and it remains unknown whether velvet proteins in anaerobic gut fungi form a complex with LaeA-like methyltransferases, as in filamentous fungi. Overall, the differential expression of some putative genes for secondary metabolites and velvet regulatory proteins provides evidence of potential life stage-specific roles.
The transcriptional and metabolic differences between zoospores and mature life stages provide insight into the developmental and cell biology of anaerobic gut fungi. In this study, we developed a method to harvest zoospores for RNA extraction and RNA-Seq. We transcriptomically compared N. californiae fungal mat and zoospore-enriched samples to inform the fundamental roles of the young and mature life stages and hypothesize on their strategies for survival. Despite being grown on glucose, zoospore-enriched samples upregulated catabolic CAZymes and putative carbohydrate transporters, supplying evidence that zoospores prime to encounter more complex plant matter substrates. Additionally, differential regulation of putative transcription factors and secondary metabolites across life stages is essential for deciphering their roles. Fig. 6. Select predicted core biosynthetic genes for secondary metabolites were differentially expressed between fungal zoospores and mats. A) A generic biosynthetic gene cluster is depicted and consists of at least one putative core biosynthetic gene and potentially one or more putative accessory gene that encodes for supporting functions, such as tailoring enzymes, membrane transporters, and mechanisms for self-resistance. B) Multiple secondary metabolite core genes were differentially expressed between mat and zoospore-enriched samples, although overall the core genes across the 39 predicted BGCs were lowly expressed, with 19 having no or low transcripts detected (average TPM ≤ 1) in both mat and zoospore samples. Three core genes (proteinIDs 502,166,28,570,and 454,096) were upregulated in mats and lowly expressed in zoospore-enriched samples (average TPM ≤ 1). Differentially regulated genes are displayed by proteinID name. PKS = polyketide synthase; NRPS = non-ribosomal peptide synthetase. Figure was created with BioRender.com.
L.V. Butkovich et al. Fungal Genetics and Biology 176 103958 Direct follow-up studies to this work could transcriptomically validate and expand upon the results in this study, analyzing zoospore-enriched pellets immediately treated with RNAlater and individually quantified for cell type content. Future studies can consider other anaerobic gut fungal strains, like the non-rhizoidal genera Caecomyces (Gold et al., 1988;Henske et al., 2017) and Cyllamyces (Ozkose et al., 2001), and other substrates, especially plant biomass substrates of different recalcitrance levels. Additionally, innovative methods to precisely isolate intermediate life stages, namely aerotolerant cysts, germlings, and thalli, would contribute to deeper understanding of this unique life cycle. Overall, the transcriptomic profiles generated in this study present a valuable resource for generating hypotheses and guiding future research into the unique biology and biotechnological potential of anaerobic gut fungi.