The Genome of the Anaerobic Fungus Orpinomyces sp. Strain C1A Reveals the Unique Evolutionary History of a Remarkable Plant Biomass Degrader

Culturing, DNA sequencing, and genome assembly. (i) Culturing.

Strain C1A was isolated from the feces of an Angus steer on a cellobiose-switchgrass medium using previously described protocols (16). Orpinomyces strain C1A was grown in an anaerobic, rumen fluid-free basal medium that was reduced by cysteine-sulfide and dispensed under a stream of 100% CO₂ as previously described (17). Cellobiose (3.75 g/liter) was used as the substrate. C1A cultures were scaled up for nucleic acid extraction in 1-liter batches prepared in 2-liter Schott bottles equipped with the stoppered top of a Balch tube to maintain strict anaerobic conditions during fungal growth. Culturing was conducted using the techniques described by Bryant and modified by Balch and Wolfe (18, 19). After autoclaving, the Schott bottles were cooled to room temperature and the gas phase was replaced by vacuuming and repressurization with 100% CO₂ (19). The medium was then amended with penicillin, streptomycin, and chloramphenicol from an anaerobic stock solution in order to provide final concentrations of 50 μg/ml, 20 μg/ml, and 50 μg/ml of each antibiotic, respectively. The medium was then prewarmed at 39°C for approximately 3 to 4 h and inoculated with 50 ml of an actively growing culture of Orpinomyces strain C1A. The cultures were incubated at 39°C for approximately 3 to 4 days, and the fungal cells were harvested during late log phase by centrifugation at 10,000 rpm for 30 min.

(ii) DNA extraction and sequencing.

(a) Illumina sequencing. High-molecular-weight genomic DNA was extracted using a modified cetyltrimethylammonium bromide (CTAB) method for isolation of nucleic acids in anaerobic fungi (20). Four micrograms of high-molecular-weight DNA was used to generate libraries for Illumina Sequence by Synthesis (Illumina-SBS) genome sequencing (21) using the standard Illumina TruSeq DNA protocol (http://genome.med.harvard.edu/documents/illumina/TruSeq_DNA_SamplePrep_Guide_15005180_A.pdf). Post-adaptor ligation size-selected fragments used for flow-cell cluster generation had a mean size of 293 bp as reported by the Agilent 3200 Bioanalyzer. Illumina sequencing was conducted using the services of a commercial provider (Ambrys Genetics, Aliso Viejo, CA, USA) on a HiSeq 2000 sequencing platform using 100-bp paired-end chemistry. Illumina sequencing yielded 29.2 Gb in 146,385,792 quality-filtered paired-end reads (see Table S1 in the supplemental material).

(b) PacBio sequencing. DNA used for Pacific Biosciences (PacBio) single-molecule real-time (SMRT) sequencing (22) was isolated using the Epicentre plant DNA extraction kit (Epicentre Corp., Madison, WI, USA) according to the manufacturer's specification. Selected inserts of 5- to 10-kb read size were prepared from 10 μg of extracted high-molecular-weight DNA by ligation to the SMRTbell sequencing adaptor. SMRT sequencing was conducted using the services of a commercial provider (Expression Analysis, Durham, NC, USA) on a PacBio RS sequencing platform using the second-generation C2 sequencing chemistry with eight zero-mode wavelength (ZMW) SMRT cells. PacBio sequencing yielded 984.8 Mb of quality-filtered data in 463,832 raw long reads (see Table S1 in the supplemental material), with an average read length of 2,124 bp. An additional 26.9 Mb of DNA sequence data was harbored in 16,949 reads (average read length of 1,586 bp) that reached a circular consensus (CCS) during sequencing and was used for long read error correction.

(iii) Genome assembly.

All computational assemblies were conducted using the SGI UV 1000 cache coherent nonuniform memory architecture (cc-NUMA) high-performance computing system Blacklight. Blacklight is an Extreme Science and Engineering Discovery Environment (XSEDE) community-shared computational resource dedicated for high-memory-footprint jobs such as de novo assembly. Blacklight is housed at the Pittsburgh Supercomputing Center (http://www.psc.edu/).

We initially attempted to utilize Illumina paired-end sequencing as the sole mechanism for C1A genome sequencing. Illumina quality-filtered reads were assembled with Velvet 1.1.07 (23), using a kmer value of 63 and a minimum coverage cutoff of 7. The resulting assembly was highly fragmented (see Table S2 in the supplemental material), with an extremely large number of contigs in the final assembly (82,325 contigs), a large proportion of the final assembly (32.4%) harbored in extremely short contigs (300 to 900 bp), and a low N₅₀ (1,666 bp).

Therefore, we sought to improve the assembly by using a hybrid SMRT-Illumina strategy that leverages short-read high-accuracy data formed from Illumina sequencing to correct errors encountered in long reads produced by SMRT sequencing (24). This hybrid approach has two steps: (i) SMRT read correction, where insertion/deletion errors present in SMRT read outputs are removed to produce corrected reads with sufficient accuracy and quality scores, and (ii) de novo assembly of the corrected SMRT reads, either independently or in conjugation with Illumina reads.

SMRT reads were corrected using the PacBioToCA package (http://sourceforge.net/apps/mediawiki/wgs-assembler/index.php?title=PacBioToCA) in wgs 7.0 (25, 26). The high-fidelity read set used for the correction was produced by using a combination of circular consensus SMRT reads, Illumina paired-end reads with sufficient overlap to merge into a single extended accurate read using fast length adjustment of short read (FLASH) (27), and Illumina paired reads without sufficient overlap for extension. Error correction resulted in a total of 570.1 Mbp, in 394,300 long corrected SMRT (C-SMRT) reads with an average phred quality score of 58.5 (28). These C-SMRT reads, which had an N₅₀ of 1,686 bp and ranged between 500 bp and 10,932 bp in length, were subsequently used for de novo assembly. All sequence data included in the final assembly had an average quality score of 59.7. The final assembly was a marked improvement compared to the Illumina-only assembly, as evident from the improved N₅₀/N₉₀ values, the increase in the number of genes with PASA transcript alignment (see below), and the increase in the average length of gene models (see Table S2 in the supplemental material). More importantly, the long C-SMRT reads allowed for the identification of a large number of introns previously undetected using Illumina assembly; the extremely low GC content (8.1%) and the large number of microsatellites within these introns probably hindered their detection and assembly from short Illumina paired-end reads.

We used the core eukaryotic genes (CEGMA) to test the completion of the final assembly (29). Due to the unique nature of C1A, e.g., absence of gluconeogenesis, anaerobic fermentative mode of metabolism, and 22/454 genes being not expected to be present, we identified 408 out of 432 genes within the final assembly, suggesting a sequence completion of ≈94.4%.

RNA sequencing and gene calling. (i) RNA sequencing, assembly, and quantitative analysis.

RNA for RNA-seq analysis (30) was isolated from a log-phase strain C1A subculture grown and propagated on rumen fluid-free basal medium with cellobiose or cellulose using the Masterpure yeast RNA purification kit (Epicentre Corp., Madison, WI, USA). RNA sequencing libraries were generated using the Illumina TruSeq RNA sample protocol. Illumina sequencing was conducted using the services of commercial providers (Ambrys Genetics, Aliso Viejo, CA, USA, using a HiSeq 2000 sequencing platform using 100-bp paired-end chemistry for cellobiose treatment, and Centrillion Biosciences, Palo Alto, CA, using the Illumina Miseq platform for cellulose treatment). Transcripts generated on cellobiose were used in gene calling efforts as described below. The number of reads and Gbp produced per sequencing run is provided in Table S1 in the supplemental material. All quality-filtered reads were assembled into transcript candidates using the de novo transcriptome assembly program Trinity (31). The Trinity de novo assembly was executed using the standard non-strand-specific library settings with the addition of the “–jaccard_clip” option, which minimizes the formation of fusion transcripts, which are often present in fungal genomes, by checking for logical spatial orientation of paired reads on the assembled transcript. Transcripts with a base pair length greater than 300 bp were considered valid candidates for downstream analysis. Expression levels for the transcriptome assembly were calculated by mapping all pair-ended RNA-seq reads with bowtie (31) using the Trinity de novo assembly as the reference index. Quantitative expression estimates of alignment values were calculated in transcripts per million values with the RNA-Seq by Expectation-Maximization (RSEM) package (32). The cellobiose C1A transcriptome was utilized for gene calling through transcript alignment. However, while genomic data were used for genomic analysis in this study, in rare cases within pathway analysis (see Results below), transcriptome data were consulted in case a single gene or a few genes within a specific pathway were inexplicably missing from the genome.

(ii) Gene calling.

Contigs larger than 1,000 bp produced by the assembly were used as input contigs for gene model generation and downstream analysis. Gene calling was conducted using a combination of ab initio gene model prediction using GlimmerHMM (33) and Augustus (34) and transcript alignments from cellobiose treatment using PASA (Program to Assemble Spliced Assemblies) (35). Training parameters for the ab initio programs were generated using the de novo-assembled transcripts aligned with the genome assembly using GMAP (36). Additional gene hint parameters available for the Augustus program using the unassembled RNA-seq read data were generated using the software's recommended protocol (http://bioinf.uni-greifswald.de/augustus/binaries/readme.rnaseq.html).

Ab initio gene-calling algorithms produced 60,595 gene models, which were combined with 14,009 PASA high-quality transcript assemblies and 38,647 Trinity transcripts for genome GMAP alignments. The final single consensus gene model for each theoretical locus was produced by information-based source-weighted integration using EvidenceModeler (EVM) (35). A final total of 16,347 consolidated consensus gene models was generated by EVM (see Table S3 in the supplemental material).

To examine the impact of assembly fragmentation on the overall number of genes identified in the C1A genome, we extracted 300 bp from both ends of every contig in the C1A final genome assembly. To identify sequence redundancy, we queried these ends against NT and NR GenBank databases using BLASTN and BLASTX, respectively. The results show a unique first hit for each contig end, suggesting that the large number of contigs in the final assembly did not result in any false duplication of predicted genes. A similar result was obtained when performing the same analysis against the COG database.

Annotation. (i) Annotation strategy overview.

Annotation of gene models and gene transcripts was achieved using a combination of command line bioinformatics programs, manual curation, and automated online annotation suites. Closest-relative homologs were assigned to each consensus gene and de novo transcript model using the BLAST+ (37) module BLASTP for the gene models and BLASTX for putative transcripts against the NR database. Identified homologs with an E value of e⁻⁴ or less were considered sufficient for evaluation of functional activity assignment. PFAM domains were identified using the hmmscan module of the HMMER suite (38) against the PFAM 26.0 database of conserved protein families (39). All domains identified by hmmscan having a full-sequence E value of e⁻⁴ or less were assigned to the gene models for functional annotation. The integrated microbial genomes system (IMG) was used for automated gene calling and annotation of both the genome gene models and the transcripts. The resulting BLAST and IMG results were manually curated and used for analysis of various cellular processes within the genome, as well as for confirmation of the presence/absence of key metabolic genes using reciprocal BLASTP against the genome. All gene calling and annotation computational work were conducted using the Oklahoma State University high-performance computing clusters Pistol Pete and Spur.

(ii) CAZyme identification and analysis.

Identification of carbohydrate active enzyme (CAZyme) genes in the C1A genome, as well as in multiple genomes that were used for comparative analysis, was achieved via PFAM domain identification and analysis as described previously (40). The CAZyme database classification system (41) was used to classify glycoside hydrolase (GH), carboxyl esterase (CE), pectate lyase (PL), and carbohydrate-binding modules (CBMs). All models harboring a dockerin domain with an e⁻⁴ significance value or less were considered eligible for further analysis and classification. Potential secreted peptides and transmembrane proteins candidates were identified using SingalP 4.0 (42) and TMHMM (43).

Additionally, more-stringent controls were conducted to guard against any possible in silico gene number inflation within the CAZyme data set due to the fragmented assembly. To guard against in silico inflation of gene numbers due to improper assembly, we identified genes with nearly complete (>97%) amino acid identity. Each cluster of genes with nearly identical amino acid sequences was aligned, and genes with less than 80% alignment with the parent model were identified as possible artifacts and removed from the assembly. To guard against possible in silico inflation of gene numbers due to gene fragmentation between two contigs, we manually examined each CAZyme gene model to identify genes with an incomplete PFAM CAZyme domain. Genes with incomplete PFAM domains were removed from the assembly. Collectively, these additional quality control approaches resulted in the removal of 76 GH, 5 CE, and 5 PL genes.

(iii) Repeat identification.

Repetitive DNA sequences (DNA repeats) are defined as sequences present in more than a defined number of copies and that have no apparent biological function (44, 45). DNA repeats can be classified into simple sequence repeats (SSRs) and complex repeats. SSRs (also known as tandem repeats) are classified, depending on the length of the repeated unit, into microsatellites (basic unit length ranging from 1 to 6 bp) (44, 45) and minisatellites (basic unit length ranging from 15 to >150 bp repeated 2 to 100 times) (45). Longer tandem repeats constitute satellites (centromeres) and telomeric repeats. Complex repeats, on the other hand, result from transposable elements (TEs). TEs are further classified into class I retrotransposons (including long terminal repeats [LTRs] and non-LTRs) and class II DNA transposons.

Microsatellite (SSRs) in the C1A genome were identified using PHOBOS (C. Mayer, Ruhr-Universität Bochum, Bochum, Germany) with the minPerfection 100 flag to detect only perfect repeats. SSRs identified have the following minimum number of repeats: mononucleotide with at least 10 repeats; dinucleotides with at least 6 repeats; and tri-, tetra-, penta-, and hexanucleotides with at least 5 repeats.

Complex repeats were identified in the C1A genome as previously described (46) using a combination of RepeatScout (47), RepeatMasker (A. F. A. Smit, R. Hubley, and P. Green, 2003; http://www.repeatmasker.org), LTR_FINDER (48), and BLASTx against RepBase (http://www.girinst.org/repbase/index.html). Briefly, a consensus repeat library was created using the default parameters of RepeatScout. This library was filtered by removing short sequences (<100 bp) and those repeats with significant hits to Uniprot proteins (except repeats with significant hits to transposable elements). The filtered consensus library was then compared to the RepBase database using BLASTx for manual annotation and classification. LTRs were identified in the genome using LTR finder. Similar to other TEs, candidate LTRs were compared to the RepBase for classification. Finally, all candidate repeats classified by RepBase were used to mask the genome in RepeatMasker to identify the number of occurrences and the percent genome coverage of each TE class.

Gene duplication was identified by running local BLASTP using C1A proteins as both the subject and the query. Only the first and second hits were examined. The second hit is a protein with similarity to the query protein present anywhere in the genome. Percent similarity cutoffs of 40% or more were used.

(iv) Noncoding RNA identification.

rRNAs were identified using local BLASTN search with sequences corresponding to published 5.8S (5.8S fragment in GenBank accession number AJ864475.1), 18S (GenBank accession number AY546684.1; Spizellomyces punctatus), 28S (28S fragment in GenBank accession number AJ864475.1), ITS1 (GenBank accession number AF170191.1), and ITS2 (GenBank accession number JN943062.1) as the database. tRNAs were identified using tRNAscan-1.4 (49).

(v) Identification of proteases, protease inhibitors, and transporters.

All proteins were compared to the MEROPS database using BLASTP to identify potential proteases and protease inhibitors. Membrane transporters in the C1A genome were identified by BLASTP comparison against the transporter classification database (TCDB) sequences available at http://www.tcdb.org/seqfile/tcdb using the GBLAST2 program (http://www.tcdb.org/labsoftware.php).

(vi) Identification of hydrogenosomal proteins.

We bioinformatically predicted proteins potentially imported to the hydrogenosomal matrix in strain C1A using a combination of motif search and Mitoprot v1.0 (50). First, C1A proteins were examined for the presence of an N-terminal mitochondrial targeting sequence corresponding to a previously predicted motif similar to ML(S|T|A|C|G|R){0,1}X{0,19}RXF(I|L|F|S|A|G|Q), ML(S|T|A|C|G|R){0,1}X{0,19}R(F|N|E|S|G) (I|L|F|S|A|G|Q), MTLX{0,19}RXF(I|L|F|S|A|G|Q), MTLX{0,19}R(F|N|E|S|G) (I|L|F|S|A|G|Q), MSLX{0,19}RXF(I|L|F|S|A|G|Q), or MSLX{0,19}R(F|N|E|S|G) (I|L|F|S|A|G|Q), where X is any amino acid except tryptophan. Numbers between braces refer to the previous residue repeat number, and parentheses indicate that any of the residues enclosed is possible at that position (7). Mitochondrial import probabilities of proteins harboring this N-terminal motif were then predicted using Mitoprot v1.0 (50), where an arbitrary probability of 0.6 was used as the cutoff. Using these criteria, we identified 21 potential intrahydrogenosomal proteins. Further, Mitoprot was also used to predict the mitochondrial import probabilities of proteins with similarity to known mitochondrial matrix proteins that did not have the above mitochondrial-targeting motif. An additional 25 potential hydrogenosomal proteins were identified using these criteria.

Comparative analysis of the C1A genome to basal fungi, Dikarya, and Opisthokonta genomes.

We used local BLASTP comparison of C1A proteins against all Mycota proteins, as well as against Dikarya proteins. A Mycota BLASTP database was created by downloading proteins of all sequenced fungal genomes available from GenBank and IMG (total of 116 fungal genomes). Of those, 4 belonged to early-branching lineages (Allomyces macrogynus, Batrachochytrium dendrobatidis, Spizellomyces punctatus, and Rhizopus oryzae). The remaining fungal genomes constitute the Dikarya BLASTP database. The number and identity of BLASTP first hits of C1A proteins against both databases at different E value cutoffs (e⁻⁵, e⁻¹⁰, e⁻¹⁵, e⁻²⁰, e⁻²⁵, e⁻³⁰, and e⁻³⁵) were then used to specify C1A proteins that are general fungal proteins (present in both early-branching and Dikarya fungi), early-branching fungus-specific proteins (present only in early-branching fungi but not in Dikarya fungi), and C1A-specific proteins (specific to C1A and absent from other fungal genomes). The last (C1A-specific proteins) were compared to the NR database excluding Mycota using the flag –negative_gilist. Proteins with no hits in the NR database were considered C1A hypothetical proteins. Functional annotation of various C1A-specific and early-branching fungus-specific proteins identified was conducted using the PANTHER classification outline (51).

Lignocellulolytic capabilities of strain C1A. (i) Plant materials and pretreatment.

Samples of mature Kanlow switchgrass (Panicum virgatum var. Kanlow), mature Sorghum bicolor, and mature energy cane (Saccharum officianarum var. Ho02) were obtained from Oklahoma State University experimental plots in Stillwater, OK. Dried alfalfa was obtained from a local farm and ranch supplier. Samples of Bermuda grass (Cynodon dactylon) were obtained from residential lawn clippings in Guthrie, OK. Samples of corn stover from Zea mays were obtained from the Industrial Agricultural Products Center at the University of Nebraska in Lincoln. Untreated wood samples, including cedar (Juniperus sp.), oak (Quercus sp.), and pine (Pinus sp.), were obtained from a local lumberyard in Stillwater, OK. Cottonwood (Populus deltoides) and willow (Salix babylonica) wood samples were harvested from live trees growing in the Stillwater area. All samples were dried at 45°C overnight, milled, and sieved to a final particle size of 2 mm as previously described (52).

Sodium hydroxide (NaOH) treatments were conducted by heating 4 g of dried plant material in 40 ml of a 1% NaOH solution inside a sealed serum bottle at 50°C for 12 h (53). Acid treatment was conducted by heating 4 g of dried plant material in 40 ml of 0.5% H₂SO₄ inside a sealed serum bottle for 1 h (54). Hydrothermolysis-treated switchgrass was prepared by mixing 60 g of switchgrass with distilled water to achieve a 10% dry matter mixture (52). This mixture was placed inside a 1-liter benchtop pressure reactor (Parr Series 4520; Parr Instrument Company, Moline, IL, USA) that was heated to 200°C and agitated at 500 rpm (52). The switchgrass-water mixture was held at 200°C for 10 min and then cooled in an ice bath (52). All of the treated switchgrass samples were recovered from pretreatment incubations by filtration. The sodium hydroxide- and acid-treated switchgrass was washed with deionized water as previously described (53, 54). All of the pretreated switchgrass samples were dried at 45°C for approximately 48 h before they were used in the experiments described below.

(ii) Growth of strain C1A on plant materials.

Experiments to evaluate the growth of strain C1A on different treated and pretreated plant materials were conducted under strict anaerobic conditions in 160-ml serum bottles. All experiments were conducted in triplicate, and unless otherwise specified, 0.5 g of plant material was used as the substrate. Experiments were conducted in a previously described rumen fluid-free basal medium (17). The medium was prepared under strict anaerobic conditions using 100% CO₂ and the techniques of Bryant (18), as modified by Balch and Wolfe (19). Once the basal medium was prepared, it was autoclaved for 20 min at 121°C and 15 lb/in² of pressure and then cooled. Each serum bottle was then amended with the appropriate type of plant biomass inside an anaerobic chamber (Coy Laboratory Products, Grass Lake, MI). After the serum bottles were amended with plant materials, they were removed from the glove bag and the headspace was repressurized with 15 lb/in² of 100% CO₂ (19). Five milliliters of an actively growing culture of strain C1A (approximately 2.6 mg of fungal biomass) was used as an inoculum and added to 45 ml medium in 160-ml serum bottles. In all experiments, serum bottles were incubated at 39°C in a nonshaking incubator. Substrate-unamended controls were included in all experiments to account for any product carryover from the inoculum. Triplicate bottles were sacrificed at different time intervals to quantify substrate loss and product formation.

(iii) Analytical methods.

Fatty acids and ethanol in supernatant fractions were quantified using a high-pressure liquid chromatograph (HPLC) with a refractive index detector (1100 series; Agilent, Santa Clara, CA, USA) and an Aminex HPX-87H column (Bio-Rad, Sunnyvale, CA, USA), which was heated to 60°C. The mobile phase was 0.01 N H₂SO₄, with a flow rate of 0.6 ml per minute. Sugars in supernatant fractions were also quantified using an HPLC with a refractive index detector (1100 series; Agilent, Santa Clara, CA, USA). The HPLC was equipped with an Aminex HPX-87P column (Bio-Rad, Sunnyvale, CA, USA), which was heated to 85°C. Distilled water was used as the mobile phase at a flow rate of 0.6 ml per minute.

The amount of plant material consumed in serum bottles was calculated by subtracting the time final from the time zero dry weights of each plant material. Since the time final pellets contained a mixture of plant and fungal biomass, the amount of fungal biomass at time final was indirectly quantified using formate concentrations as previously described (12). The amounts of cellulose, xylan, hemicellulose, and lignin in the different plant substrates were determined using the standard NREL procedures (55). The procedure included the addition of 3 ml of 72% sulfuric acid to each sample and incubation at 30 ± 3°C for 1 h, with stirring every 5 to 10 min. The samples were then diluted with 84 ml of deionized water, capped, and autoclaved for 1 h to 121°C. The cooled solution was filtered, and this filtrate was used to determine carbohydrate content and soluble lignin. The remaining solids were washed and dried to constant weight at 105°C to determine acid-insoluble residue (AIR) and then converted to ash at 575°C for 24 h (55). Analyses of resulting carbohydrates within the filtrate were done by HPLC with refractive index detection (RID) (Agilent 1100 series; Santa Clara, CA) on an Aminex HPX-87P column at 85°C with a mobile phase of deionized water pumped at 0.6 ml/min for 30 min (55). Twenty microliters of each sample was analyzed for cellobiose, glucose, xylose, galactose, arabinose, and mannose. Contributions of structural constituents to the total biomass composition were determined using the NREL summative mass closure procedure (56). The acid-soluble lignin (ASL) content was determined using a UV spectrophotometer set at a wavelength of 205 nm, as has been previously used to determine ASL in switchgrass (57). As recommended in the NREL procedure, ASL in corn stover was measured at 320 nm, whereas a 240-nm wavelength was used for the remaining biomass types (55).

Nucleotide sequence accession numbers.

The final genome assembly is available in the IMG genome database with accession number 2518645524 and in the GenBank database with accession number PRJNA200719. The final transcriptome is available in the IMG database with accession number 2510461071.

Isolation and general genomic features.

Strain C1A was isolated from the feces of an Angus steer on a cellobiose-switchgrass medium using previously described protocols (16). The isolate displayed polycentric growth and effectively colonized switchgrass. Phylogenetic analysis using the nuclear ribosomal internal transcribed spacer II (ITS-II) region supported the placement of strain C1A as a member of the genus Orpinomyces, and phylogenetic analysis using a concatenated set of 42 housekeeping genes supported the basal, early-diverging position of the Neocallimastigomycota (see Fig. S1 in the supplemental material).

We sequenced the C1A genome using a combination of paired-end short-read Illumina technology (≈290× coverage) and single-molecule real-time (SMRT) Pacific Biosciences technology (≈10× coverage). The transcriptomes of strain C1A grown on cellobiose and cellulose were also sequenced using Illumina technology (see Tables S1 and S2 in the supplemental material). The C1A genome displayed several interesting features (Fig. 1; see also Fig. S2 and S3 and Tables S4 to S8). It had the lowest GC content (17.0%) compared to available genomes of all free-living microorganisms sequenced to date (Fig. 1A). This value is lower than those observed in the notoriously AT-rich Dictyostelium spp. and Plasmodium sp. within the microeukaryotes and is surpassed only by a few proteobacterial obligate endosymbionts, e.g., “Candidatus Zinderia insecticola” (13.5%) and “Candidatus Carsonella ruddii” (16.6%). The GC content was higher in protein-coding genes (26.8%) than in noncoding regions (14.8% in intergenic regions and 8.1% in introns) but still resulted in a marked codon usage skew (see Table S6). The C1A genome was also characterized by a relatively large proportion of noncoding intergenic regions (73.1%) (Fig. 1B). Further, noncoding regions displayed massive proliferation of simple sequence repeats (SSRs) in the C1A genome. The 249,194 SSRs constituted 4.9% of the entire genome, as well as 5.8% and 6.0% of the intergenic region and the introns, respectively (Fig. 1C). These values vastly surpass the number of SSRs that were observed in previously analyzed fungal genomes by at least 1 order of magnitude (46, 58). Homopolymeric A or T monorepeats represented the majority of observed repeats (68.6% and 60.6% of total SSR numbers and length, respectively), with 3,589 identified cases of ≥50-bp stretches of A or T within the assembly (Fig. 1D).

Comparative gene content with basal and Dikarya fungi.

Gene calling resulted in the identification of 16,347 protein-coding genes, a number surpassed by only a few fungal genomes (see Text S1 and Table S5 in the supplemental material). This large number could partly be attributed to gene duplication (Fig. 2A), since 3,252 gene pairs share >90% sequence similarity.

In addition, comparative genomic analysis indicated that only 48.4% of C1A genes have at least one ortholog in all examined Dikarya genomes (n = 116), that 9.5% of C1A genes have at least one ortholog within examined early-branching fungal genomes (n = 4) but not in Dikarya genomes, and that 42.2% (6,886) of C1A genes are unique and have not been previously encountered within the Mycota (Fig. 2B; see also Table S9 in the supplemental material), These unique C1A genes were either genes for hypothetical proteins (n = 5,666); genes with nonfungal, eukaryotic orthologs (n = 578); or genes with bacterial orthologs (n = 642). Eukaryotic, nonfungal C1A genes were mostly encoding cellular processes, e.g., receptors and nucleic acid-binding proteins, highlighting the distinct early-branching fungal position of the Neocallimastigomycota, while C1A genes with prokaryotic orthologs were mainly involved in metabolic processes, e.g., hydrolases, transporters, transferases, and phosphatases, highlighting the potential role of horizontal gene transfer in shaping C1A metabolic capabilities.

Genomic analysis and comparative genomics reveal multiple differences between Neocallimastigomycota and Dikarya.

Detailed analysis of metabolic capabilities, cellular processes, and structural features in the C1A genome is presented as supplemental material (see Text S1 and Fig. S4 in the supplemental material). Briefly, analysis of genes involved in information processing (replication, transcription, and translation), as well as cytoskeletal structure and intracellular trafficking mechanisms, revealed all salient features associated with such processes in eukaryotic cells.

More importantly, comparative genomic analysis identified multiple cellular processes in which either only the C1A genome or all early-branching fungal genomes possess features that appear to be absent from Dikarya genomes but are mostly associated with nonfungal Opisthokonta (choanoflagellates such as Monosiga brevicollis and Capsaspora owczarzaki, an independent unicellular Opisthokonta lineage, and Metazoa) and higher non-Opisthokonta eukaryotes (Table 1; see Text S1 and Tables S10 and S11 in the supplemental material). Five different examples are highlighted. (i) One gene in the C1A genome encodes a member of the metalloprotease site-2-protease (S2P) family, and seven different genes encode various components of the γ-secretase complex, including aspartyl protease presenelin. Both of these types of intramembrane proteases are represented in Metazoa (mostly Chordata, Nematoda, and Arthropoda), with a few representatives in plants, and have no representation in Dikarya (Table 1). (ii) The genome contains two fucosyltransferase genes that mediate fucosylation, a posttranslational modification process that is typically observed in Chordata, Arthropoda, and Viridiplantae but not in the Dikarya (Table 1). (iii) The C1A genome possesses a near-complete focal adhesion (FA) machinery (Table 1). FAs are large multiprotein intracellular assemblies that mediate cell anchorage and mechanical adhesion to the extracellular matrix. They also act as a signaling milieu where signaling proteins are concentrated at sites of integrin binding and connect the cell's cytoskeleton to the extracellular matrix. FA appears to be absent from filamentous fungi and more common in other eukaryotes (Amoebozoa, Metazoa [sponges, placozoans, and cnidarians], and Holozoa) (59). (iv) The C1A genome possesses a complete axoneme and intraflagellar trafficking machinery proteins. The axoneme acts as a scaffold for other protein complexes, including motor proteins (e.g., kinesin and dynein) essential for intraflagellar transport of proteins. Ciliated and flagellated eukaryotic cells are known to possess an axoneme, as do Neocallimastigomycota and other early-branching fungi that produce flagellated zoospores. This feature is absent in other Dikarya fungi that produce nonflagellated spores. (v) Finally, the C1A genome encodes various extracellular protease inhibitors, some of which (serpins) have not previously been encountered in the Dikarya (e.g., serpins of family I4 are present mostly within eukaryotic metazoan phyla [Arthropoda, Chordata, and Nematoda] and have also been identified in Bacteria and Archaea but have not previously been encountered in fungi). Several identified serpins have dockerin domains, confirming their cellulosomal destination and their potential role in combating plant proteases, as previously suggested (60).

Hydrogenosomal structure and function.

Anaerobic fungi lack mitochondria but possess a double-membrane hydrogenosome whose main function is ATP production via substrate-level phosphorylation and hydrogen production (6, 8, 9). The C1A genome encodes a near-complete hydrogenosomal protein import system with components of the TOM outer membrane transport system (4 out of 7 genes), the SAM sorting and assembly complex (4 out of 5 genes) for protein insertion in the outer membrane, the MIA intermembrane space import and assembly complex (2 out of 3 genes), small TIMs (2 out of 4 genes), the TIM22 complex for protein insertion in the inner membrane (6 out of 6 genes), and the inner membrane transport system and associated motor (TIM23 complex, 10 out of 11 genes) (Fig. 3). In comparison, the hydrogenosomal import machinery of Trichomonas vaginalis has been reduced to a few outer membrane proteins (Tom40, Sam50, Hmp35, and Hmp36), a few inner membrane proteins (Tim17/22/23, Tim44, PAM16, and PAM18), and one highly modified intermembrane small TIM (61). Further, examining the phylogenetic affiliation of mitochondrial import proteins in strain C1A clearly demonstrates their fungal origin, since their closest relatives are consistently those from fungal mitochondria.

Using two different bioinformatic criteria, we identified 46 intrahydrogenosomal proteins in the C1A genome. Candidate proteins for import into the hydrogenosomal matrix included several hypothetical proteins, Fe-S cluster assembly and maturation proteins, peptidases, and intraluminal chaperones and cochaperones, as well as pyruvate metabolism and energy production enzymes (see Tables S12 to S14 in the supplemental material).

Hydrogenosomes are the sites of multiple metabolic processes for pyruvate metabolism, ATP production via substrate-level phosphorylation, and regeneration of reduced electron carriers, e.g., NADPH and NADH. The C1A genome carries the genes required for mixed-acid fermentation, the predominant pathway for pyruvate metabolism in anaerobic fungi, as previously suggested (62–64) (Fig. 3; see also Table S13 in the supplemental material). Genomic reconstruction suggests that pyruvate produced from sugar catabolism in the cytosol could be either metabolized cytosolically or imported and metabolized in the hydrogenosome. In the cytosol, pyruvate could be either converted to acetyl coenzyme A (acetyl-CoA) and formate via cytosolic pyruvate formate lyase (PFL) (a cytosolic PFL-activating enzyme is also encoded by the genome), converted to d-lactate via cytosolic d-lactate dehydrogenase, or used to produce tricarboxylic acid (TCA) intermediates required for anaplerotic reactions via an incomplete cytosolic TCA cycle. In the cytosol, acetyl-CoA produced could be converted to ethanol via aldehyde dehydrogenase/alcohol dehydrogenase. In the hydrogenosome, pyruvate could be imported from the cytosol or could be produced from malate via the action of the hydrogenosomal malic enzyme (with the production of CO₂). Hydrogenosomal pyruvate could be then metabolized to acetyl-CoA and formate by a hydrogenosomal PFL. Acetyl-CoA produced in the hydrogenosome could subsequently be converted to acetate via the combined action of hydrogenosomal acetate:succinate CoA transferase/succinyl-CoA synthase to produce ATP via substrate-level phosphorylation. The genome encodes a hydrogenosomal acetyl-CoA hydrolase. A similar enzyme in T. vaginalis was shown to possess an acetate:succinate CoA transferase activity. The Orpinomyces acetyl-CoA hydrolase homolog is most likely performing a similar transferase activity that, coupled to the succinyl-CoA synthase activity, could convert acetyl-CoA to acetate.

In addition, the hydrogenosomal components contain elements for NADH recycling coupled to H₂ production (Fe-only hydrogenase large subunit and NADH dehydrogenase [complex I] subunits E and F). The 2 subunits of NADH dehydrogenase most probably function to reoxidize NADH produced in the lumen (e.g., during fatty acid degradation) and transfer electrons to Fe-only hydrogenase. Since PFL mediates pyruvate metabolism without the production of reduced equivalents, H₂ production via the hydrogenase enzyme is thought to be minor and to be required only to cope with the NADH produced from other intraluminal reactions, e.g., 3-hydroxyacyl-CoA dehydrogenase of fatty acid metabolism (65). Elements of hydrogenosomal NADPH recycling are also present in the genome. NADPH produced from NADP-dependent reactions, e.g., malic enzyme, could possibly be used by the NADP-requiring fatty acid synthesis reactions, e.g., 3-oxoacyl-[acyl-carrier protein] reductase, or recycled by NADPH:quinone reductase to a quinone, where electrons could then transfer to succinate dehydrogenase to reduce fumarate to succinate. Finally, the genome also encodes subunits α, β, γ, and δ of F_oF₁-type ATP synthase that is thought to pump protons to the cytosol, keeping the luminal pH slightly alkaline. ATP synthase is likely functioning in conjunction with an ADP/ATP carrier.

Lignocellulolytic repertoire of strain C1A.

Prior research efforts have identified multiple genes involved in plant biomass degradation in several Neocallimastigomycota isolates (5, 14, 15, 66–88). To provide an overall view of the plant biomass degradation machinery of an anaerobic rumen fungus, we analyzed the lignocellulolytic machinery in the C1A genome. Such an analysis revealed an extremely rich repertoire that consisted of 357 glycoside hydrolase (GH) genes, 24 polysaccharide lyases (PLs), and 92 carbohydrate esterases (CEs) (Fig. 4A; see also Tables S18 to S20 in the supplemental material). Principal-coordinate analysis demonstrated the unique position of the GH catalytic machinery, compared to multiple fungal and bacterial genomes of distinct habitats, ecological roles, phylogenetic affiliation, and oxygen preferences (Fig. 4B). For instance, compared to aerobic fungal biomass degraders of industrial and ecological relevance such as Trichoderma reesei, Postia placenta, Aspergillus oryzae, and Myceliophthora thermophila, the C1A genome shows an expansion of cellulolytic families GH6, GH9, GH45, and GH48 and hemicellulolytic families GH10, GH11, and GH43, as well as the reduction or absence of families GH7, GH16, GH18, GH28, and GH61.

Detailed phylogenetic analysis (see Fig. S5 in the supplemental material) suggests that the GH machinery in strain C1A has evolved from an ancestor with relatively limited cellulolytic capability to a robust cellulolytic and hemicellulolytic organism through the acquisition of genes from multiple bacterial lineages, many of which are known to be prevalent in the bovine rumen. Overall, 247 (69.2%) of GH genes were most closely related to bacterial orthologs, and 141 (39.5%) of GH genes were most closely related to bacterial orthologs from lineages that are prevalent in the bovine rumen. Such lineages include the families Lachnospiraceae, Clostridiaceae, Eubacteriaceae, and Ruminococcaceae within the order Clostridiales, the family Streptococcaceae within the order Bacillales, and the family Prevotellaceae within the order Bacteroidales, as well as the phylum Fibrobacteres. Cellulose degradation machinery in strain C1A consists of GH5, GH8, GH9, and GH45 endoglucanases and GH6 and GH48 cellobiohydrolases. GH8 is an exclusively prokaryotic gene family (41), and phylogenetic analysis of GH5 and GH9 endoglucanases indicates their close affiliation with endoglucanases from multiple bacterial sources, including the ruminal genera Clostridium, Ruminococcus, and Eubacterium. On the other hand, strain C1A also possesses the distinctively eukaryotic fungus-affiliated GH45 endoglucanases that have rarely been observed in bacterial genomes. Similar to endoglucanases, a dual prokaryotic/eukaryotic origin of strain C1A cellobiohydrolases was observed. Strain C1A possesses multiple GH48 reducing-end cellobiohydrolases, hallmarks of cellulosomal cellobiohydrolases, as well as multiple GH6 non-reducing-end cellobiohydrolases that are common in fungal genomes but rarely observed in anaerobic cellulolytic bacteria.

Unlike cellulose metabolism, hemicellulose degradation machinery in strain C1A appears to be entirely of prokaryotic origin. The C1A genome contains all genes required for the degradation of xylans (glucuronoarabinoxylans and arabinoxylans), mannans (galactoglucomannans and glucomannans), and mixed β-(1,3-1,4) glucans. Strain C1A appears to be highly adapted to the degradation of xylans, the prevalent hemicelluloses in grasses (order Poales) (89). This is evident by the identification of 109 different xylanases, xylosidases, arabinofuranosidases, and α-glucuronosidases belonging to families GH10, GH11, GH39, GH43, and GH67, in addition to multiple glucuronoarabinoxylan- and arabinoxylan-debranching enzymes (acetylxylan esterases, ferulic acid esterases, and polysaccharide deacetylases). Phylogenetic analysis of GH10 xylanases suggests their close affiliation with multiple bacterial lineages, including the ruminal genera Butyrivibrio, Clostridium, and Eubacterium. Phylogenetic analysis of GH11 xylanases suggests that they have been solely acquired from Fibrobacter species, important constituents of rumen microbiota. A similar bacterial origin was also observed for GH39 and GH43 xylosidase/arabinofuranosidases as well as GH67 α-glucuronosidases, with potential bacterial donors being members of the genera Clostridium, Ruminococcus, Butyrivibrio, Cellulosilyticum, Eubacterium, and Prevotella. Finally, GH26 mannosidases and GH16 β-(1,3-1,4)-glucanase, mediating the breakdown of mannans and mixed glucans, also had similar bacterial origin, with several sequences affiliated with ruminal lineages, e.g., Acetovibrio, Fibrobacter, and Streptococcus.

Anaerobic fungi produce cellulosomes: extracellular structures that harbor multiple extracellular enzymes bound to scaffoldins (2). Cellulosome-bound genes in anaerobic fungi usually harbor a fungal dockerin domain (FDD) that is similar in structure to carbohydrate-binding module family 10 (CBM10) (5). We identified a total of 220 genes (see Table S21 in the supplemental material) with FDDs, 108 of which contained dual glycoside hydrolase-fungal dockerin domains (GH-FDDs). GH-FDD genes identified suggest that cellulosome-bound enzymes play a role in the degradation of cellulose and hemicellulose but not chitin, starch, or pectin. Within the remaining FDD-containing genes, we identified multiple putative activities that could either aid in biomass degradation (e.g., polysaccharide deacetylases, tannase, lipases, swollenin, and expansin module proteins) or act as cellulosomal preservation and defense mechanisms (e.g., protease inhibitors [serpins]), as well as multiple conserved hypothetical and hypothetical proteins.

Carbohydrate-binding modules (CBMs) are noncatalytic domains that are often encountered in lignocellulolytic enzymes and promote the association of the enzyme with the substrate. A total of 103 genes harboring carbohydrate-binding module (CBM) domains belonging to 6 different CBM families were identified (see Table S22 in the supplemental material). The majority (75.7%) of CBMs were members of the exclusively fungal CBM1 domain. Twenty-six genes with CBM domains were associated with GH enzymes, and 7 were associated with PL enzymes (see Table S22). Within GH-CBM dual domain genes, CBM1 domains were associated with several GH10 and GH11 xylanases, CBM18 domains were associated with GH18 chitinases, and CBM48 domains were associated with GH13 amylases. No CBM domains were identified in GH genes putatively involved in cellulose metabolism in the C1A genome. No CBM2 or CBM3 domains, the prevalent CBMs in bacterial plant biomass-degradation genes and in rumen anaerobic cellulosomal bacteria, respectively, were identified in the C1A genome (see Table S22).

Comparative transcriptomic analysis of strain C1A was conducted on cellobiose-grown versus microcrystalline cellulose-grown cultures (Fig. 5; see also Fig. S6 and Tables S23 and S24 in the supplemental material). A total of 172 GH genes were expressed under both conditions, while 39 and 4 GH genes were identified only in cellobiose-grown and cellulose-grown cultures, respectively. In cellulose-grown cultures, transcripts belonging to GH5 cellulases, as well as GH9 and GH48 cellobiohydrolases, were drastically upregulated compared to cellobiose-grown cultures. GH8 and GH45 cellulases were only slightly upregulated, and their overall transcriptional levels were relatively low (Fig. 5A and B). GH1 and GH3 β-glucosidases, essential for substrate degradation under both conditions, were either not significantly changed or only slightly upregulated in cellulose-grown cultures (Fig. 5A and B).

Analysis of expression profiles of all GH genes identified under both conditions revealed that while several cellulase and cellobiohydrolase genes were clearly upregulated in cellulose-grown cultures, the majority of such genes were not significantly (>2-fold) affected by the growth condition (Fig. 5C), and few were even significantly downregulated. Interestingly, few of the genes upregulated in cellulose-grown cultures belong to GH families associated with the degradation of plant polymers other than cellulose, e.g., GH10 and GH11 xylanases, GH18 chitinases, and GH26 mannosidases (see Fig. S6 in the supplemental material).

Strain C1A is an effective, versatile biomass degrader.

Strain C1A effectively metabolized a variety of sugars and polysaccharides, including crystalline cellulose and xylan (see Fig. S7 in the supplemental material). More importantly, strain C1A grew readily on untreated, as well as mild acid-, mild alkali-, and hydrothermolysis-treated switchgrass, with the concurrent utilization of cellulose and hemicellulose fractions, but not lignin (Fig. 6; also see Tables S25 and S26 and Fig. S8 and S9). Dry weight losses of substrate ranged from 18.6% (28.7% of nonlignin fraction) in untreated switchgrass to 40.8% (53.9% of nonlignin fraction) in NaOH-treated switchgrass. Further, adjustments to the inoculum/substrate ratios resulted in an increase in the amount of switchgrass metabolized up to 42.8% and 58.4% of the dry weight of untreated and NaOH-treated switchgrass, respectively. Strain C1A performed extremely well on NaOH-treated switchgrass, since this method of pretreatment retains the majority of the hemicellulose content (90, 91), which is degradable by strain C1A. Strain C1A also grew well on hydrothermolysis-treated switchgrass, presumably due to the fact that the removal of hemicellulose resulted in a greater accessibility to cellulose fibers. End product analysis indicated that lactate, acetate, and formate are the main end products of plant biomass degradation. Only minor amounts of ethanol were produced, ranging between 0.045 and 0.096 mg ethanol/mg biomass (Table S26).

In addition to switchgrass, we tested the capability of strain C1A to utilize several other types of energy crops (e.g., alfalfa, sorghum, and energy cane), agricultural residues (e.g., corn stover), and grasses (e.g., Bermuda grass). We chose these specific plant materials due to the variations in the percentages of cellulose, hemicellulose, and lignin in these plants. The results (Fig. 6E; see also Tables S25 and S26 in the supplemental material) demonstrate the versatility of strain C1A, since it was able to metabolize all different types of examined plant biomass. Within both untreated and NaOH-treated experiments, strain C1A was most effective in metabolism of corn stover, with 40.6% and 62.3% dry weight loss, 51.0% and 75.8% loss in cellulose fraction, and 43.0% and 74.3% loss in hemicellulose fractions in untreated and NaOH-treated corn stover, respectively.