Genomic analysis of methanogenic archaea reveals a shift towards energy conservation

14 Pages·2017·2 MB·English

by  Sean P. Gilmore

Gilmoreetal.BMCGenomics (2017) 18:639
DOI10.1186/s12864-017-4036-4
RESEARCH ARTICLE Open Access
Genomic analysis of methanogenic
archaea reveals a shift towards energy
conservation
Sean P. Gilmore1, John K. Henske1, Jessica A. Sexton1, Kevin V. Solomon1,6, Susanna Seppälä1,2, Justin I Yoo1,
Lauren M. Huyett1, Abe Pressman1, James Z. Cogan3, Veronika Kivenson4, Xuefeng Peng1,4, YerPeng Tan5,
David L. Valentine4 and Michelle A. O’Malley1*
Abstract
Background: Themetabolismofarchaealmethanogensdrivesmethanereleaseintotheenvironmentandiscriticalto
understandingglobalcarboncycling.Methanogenesisoperatesataverylowreducingpotentialcomparedtoother
formsofrespirationandisthereforecriticaltomanyanaerobicenvironments.Harnessingoralteringmethanogen
metabolismhasthepotentialtomitigateglobalwarmingandevenbeutilizedforenergyapplications.
Results:Here,wereportdraftgenomesequencesfortheisolatedmethanogensMethanobacteriumbryantii,
Methanosarcinaspelaei,Methanosphaeracuniculi,andMethanocorpusculumparvum.Theseanaerobic,methane-
producingarchaearepresentadiversesetofisolates,capableofmethylotrophic,acetoclastic,andhydrogenotrophic
methanogenesis.Assemblyandanalysisofthegenomesallowedforsimpleandrapidreconstructionofmetabolismin
thefourmethanogens.ComparisonofthedistributionofClustersofOrthologousGroups(COG)proteinstoasample
ofgenomesfromtheRefSeqdatabaserevealedatrendtowardsenergyconservationingenomecompositionofall
methanogenssequenced.Furtheranalysisofthepredictedmembraneproteinsandtransportersdistinguisheddiffering
energyconservationmethodsutilizedduringmethanogenesis,suchaschemiosmoticcouplinginMsar.spelaeiand
electronbifurcationlinkedtochemiosmoticcouplinginMbac.bryantiiandMsph.cuniculi.
Conclusions:Methanogensoccupyauniqueecologicalniche,actingastheterminalelectronacceptorsinanaerobic
environments,andtheirgenomesdisplayasignificantshifttowardsenergyconservation.Thegenome-enabled
reconstructedmetabolismsreportedherehavesignificancetodiverseanaerobiccommunitiesandhaveledto
proposedsubstrateutilizationnotpreviouslyreportedinisolation,suchasformateandmethanolmetabolismin
Mbac.bryantiiandCO metabolisminMsph.cuniculi.Thenewlyproposedsubstratesestablishanimportant
2
foundationwithwhichtodecipherhowmethanogensbehaveinnativecommunities,asCO andformateare
2
commonelectroncarriersinmicrobialcommunities.
Keywords:Methanogenesis,Archaea,Metabolism,Anaerobes,Energy
*Correspondence:[email protected]
1DepartmentofChemicalEngineering,UniversityofCalifornia,SantaBarbara,
California,USA
Fulllistofauthorinformationisavailableattheendofthearticle
©TheAuthor(s).2017OpenAccessThisarticleisdistributedunderthetermsoftheCreativeCommonsAttribution4.0
InternationalLicense(http://creativecommons.org/licenses/by/4.0/),whichpermitsunrestricteduse,distribution,and
reproductioninanymedium,providedyougiveappropriatecredittotheoriginalauthor(s)andthesource,providealinkto
theCreativeCommonslicense,andindicateifchangesweremade.TheCreativeCommonsPublicDomainDedicationwaiver
(http://creativecommons.org/publicdomain/zero/1.0/)appliestothedatamadeavailableinthisarticle,unlessotherwisestated.
Gilmoreetal.BMCGenomics (2017) 18:639 Page2of14
Background methanogenic species are of great interest, as is eluci-
Methanogenic archaea are key players in anaerobic dating the details of ATP generation via carbon dioxide
communities and as such contribute largely to global reduction. Comparing and contrasting the metabolic
warming and energy production, with interesting po- pathways in methanogenic species of diverse evolutionary
tential roles in human and animal health [1, 2]. These originswillbevitaltoprogressinthisarea.
extraordinary microbes perform a type of anaerobic Currently, methanogens are grouped in three ways: by
respiration known as methanogenesis by reduction or energy requirements, based on phylogeny, or by catalytic
dismutation of carbon dioxide, methyl compounds, or abilities. Metabolic grouping depends on energy source,
acetate to methane or methane and carbon dioxide in which can be i) hydrogen/formate and carbon dioxide,
several ecosystems and consortia [3]. Methanogenesis ii) methyl compounds, or iii) acetate [1]. Phylogenetic
operates at a very low reducing potential compared to groupingbasedon16SrRNA sequencingcomprisesonly
other forms of aerobic and anaerobic respiration [4]. As two classes: Class I and Class II [16]. This classification
such, methanogens can be found in such diverse envi- also reflects metabolic differences, with basal Class I lin-
ronments as the deep ocean, rice paddies, wetlands, eages requiring hydrogen and carbon dioxide, and some-
landfills, and the gastrointestinal tracts of termites, ru- times also formate, and the more recently evolved Class
minants, and humans, where more favorable electron II lineages requiring either methyl compounds or acetate
acceptors like oxygen or inorganic ions are unavailable [17, 18]. Recently, Anderson et al. conducted a new type
[1]. Because they constitute a major source of atmos- of phylogenetic analysis and proposed that Class II be
pheric methane, there has been a drive to increase defined to include only species of the Methanomicro-
understanding of methanogen metabolism and their in- biales order, while the Methanosarcinales would form a
teractions with other organisms [5]. This understanding new class: Class III [18]. Recently, new genomes and
of methanogen metabolism is necessary for optimal de- metagenomes sequenced have identified new clades of
sign of biochemical processes aimed at generating or methanogens outside of the early-established class sys-
removing methane [6, 7]. tems, such as the Methanomassiliicoccales [19] or the
Many studies seek to mitigate negative effects con- newly discovered Methanonatronarchaeia [20].
tributed by methanogens. Livestock, as a component of We sought to better understand metabolic capabilities
agriculture, are the largest anthropogenic source of at- of methanogens and the biochemical pathways behind
mospheric methane [8]. Certain methanogens thrive in them. To this end, we obtained four methanogenic
the rumen compartment of cattle, where they partner species from diverse environments and sequenced their
with cellulose-degrading gut microbes, consuming the genomes. Under the three class system [18], Methano-
hydrogen and carbon dioxide these microbes produce bacterium bryantii and Methanosphaera cuniculi are
as waste [9]. This phenomenon, known as enteric emis- Class I methanogens requiring H /CO and H /metha-
2 2 2
sion, produces approximately 30% of global methane nol, respectively. Methanocorpusculum parvum is a
emissions, a greenhouse gas 25 times more potent than member of Class II and requires H /CO , formate, or
2 2
carbon dioxide [5, 8]. The effect of methanogens on 2-propanol/CO . Finally, Methanosarcina spelaei is a
2
human and animal health is less clear. Although there member of Class III and utilizes H /CO , acetate, or
2 2
are conflicting reports, many studies describe a direct methyl compounds for growth and methanogenesis. By
correlation between the presence of methanogens in sequencing the genomes of each member and recon-
human guts and the occurrence of diseases such as structing their metabolisms, we seek to elucidate the
periodontitis,obesity,andinflammatoryboweldisease[2]. metabolic pathways and mechanisms implemented by
The high energy content of methane has recently each of the three classes and their contributions to
called attention to methanogens as a tool in waste-to- energy, health, and agriculture.
energy technologies [7, 10]. As in the rumen and other
natural anaerobic environments, methanogens can part- Methods
ner with other microbes in specialized anaerobic reac- MethanogencultureandDNAisolation
tors to convert organic biomass into methane [11]. This Methanobacterium bryantii strain M.o.H. (DSM 863),
process, termed anaerobic digestion, results in gener- Methanosarcina spelaei strain MC-15 (DSM 26047),
ation of potentially-useful biogas [11]. Recently, anaer- Methanosphaera cuniculi strain 1R-7 (DSM 4103), and
obic microbes have received attention for their ability to Methanocorpusculum parvum strain XII (DSM 3823)
thrive on crude biomass [12–14]. Methanogens enable wereorderedfromtheLeibnizInstituteDSMZ–German
the growth of other microbes by removing self-limiting Collection of Microorganisms and Cell Culture
waste products, which methanogens assimilate or use as (https://www.dsmz.de/). Cultures were grown in a
electron donors/acceptors inmethanogenesis[3,15].For modified version of M2 medium adapted from Teunis-
this reason, the substrate utilization capabilities of sen et al. [21]. The base medium contained 150 mL of
Gilmoreetal.BMCGenomics (2017) 18:639 Page3of14
Solution A (0.45 g/L KH PO , 0.45 g/L (NH ) SO , Librarypreparationandsequencing
2 4 4 2 4
0.9 g/L NaCl, 0.09 g/L MgSO ·7 H O, 0.09 g/L CaCl ·2 Genomic DNA (gDNA) was prepared for high through-
4 2 2
H O), 150 mL of Solution B (0.45 g/L K HP0 ), 12 g/L put sequencing (HTS) using the TruSeq DNA PCR-
2 2 4
NaHCO , 1 mL of Vitamin Supplement (ATCC MD- Free library prep kit supplied by Illumina, Inc. (San
3
VS) supplemented with 0.9 μg cyanocobalamin and Diego, CA). Briefly, purified gDNA was first fragmented
40 mg CoM, 10 mL trace element solution (2.5 mg/mL usingaCovaris(Woburn,Massachusetts)M220Focused
MnCl ·4H O, 2.5 mg/mL NiCl ·6H O, 2.5 mg/mL Ultrasonicator, followed by end repairs, size selection
2 2 2 2
NaMoO ·2H O, 2.5 mg/mL H BO , 2.0 mg/mL FeS- (~330 bp), end adenylation and paired-end adapters
4 2 3 3
O ·7H O, 0.5 mg/mL CoCl ·6H O, 0.5 mg/mL SeO , ligation using the kit. Prepped libraries were then quan-
4 2 2 2 2
0.5 mg/mL NaVO ·4H O, 0.25 mg/mL ZnCl , 0.25 mg/ tified using Qubit (Life Technologies, Carlsbad, CA)
3 2 2
mL CuCl ·2H O, all dissolved in 0.2 M-HCl), 10 mL and TapeStation (Agilent, Santa Clara, CA), before
2 2
hemin solution (1 mg/L hemin, 0.01% (v/v) ethanol, pooling. HTS was performed with an Illumina Next-
dissolved in 0.05 M NaOH), 1 mL of resazurin solution Seq500 sequencer using a 150 cycle, mid output kit
(0.1% (w/v) resazurin), 1 g L-cysteine-HCl [22]. The (2 × 75 paired-end). Resequencing of Methanosarcina
total volume of the medium was adjusted to 1 L with spelaei was completed using a 75 cycle mid output kit
MilliQ H O, boiled vigorously to drive out the oxygen, (75 bp single-end).
2
cooled under CO , and aliquoted under an 80/20 mix-
2
ture of H /CO before autoclaving. Vitamin solution Genomeassemblyandannotation
2 2
and methanol were filter sterilized and added after Genomes were assembled from 75 bp paired end Illu-
autoclaving. Mbac. bryantii medium was supplemented mina NextSeq reads (Illumina, San Diego, CA) at a
with 2 g/L sodium acetate, 4 g/L sodium formate, 2 g/L minimum of 49× coverage and annotated with the De-
yeast extract, and 4 g/L bacto casitone. Mcor. parvum partment of Energy Systems Biology Knowledgebase
medium was supplemented with 2 g/L sodium acetate, (KBase, http://kbase.us) automated pipeline. Briefly,
4 g/L sodium formate, 2 g/L yeast extract, 4 g/L bacto reads were preprocessed with BayesHammer [23] before
casitone, and 10 μM sodium tungstate. Msar. spelaei being assembled de novo into contigs with the SPAdes
medium was supplemented with 4.5 g/L sodium chloride, and Velvet genome assembly algorithms [24, 25]. Gen-
2g/Lsodiumacetate,2g/Lyeastextract,2g/Lbactocasi- omic features including ORFs, large repeat regions,
tone, and 0.5% (v/v) methanol. Msph. cuniculi medium rRNAs, CRISPRs, and tRNAs were then identified and
was supplemented with 3.6 g/L sodium acetate, 2 g/L annotated with the Rapid Annotations using Subsys-
yeastextract,2g/Lbactocasitone,and1%(v/v)methanol. temsTechnology toolkit (RASTtk) [26]. These gene an-
Allculturesweregrownat39°Cwithoutshaking.Experi- notations were combined with biochemical information
mentstestinggrowthonsubstratesalternativetoH were from the Kyoto Encyclopedia of Genes and Genomes
2
conductedin100%CO headspace. (KEGG) [27] to reconstruct the metabolism of each
2
To isolate genomic DNA, 1 mL of methanogens was methanogen. The assembly for Methanosphaera cuni-
inoculated into 40 mL of media in 60 mL Wheaton culi contained a number of low coverage (~1X) contigs,
serum bottles until stationary phase (OD ~ 0.2–0.5) which were removed from the draft assembly and ana-
600
and then harvested by centrifugation for 30 min at lysis. Scaffolding of assembled contigs was performed
10,000×g at 4 °C. Cell pellets were resuspended in using SSPACE [28].
0.5 mL TE Buffer (10 mM Tris, 1 mM EDTA, pH 8.0).
Sodiumdodecylsulfatewasaddedtoafinalconcentration COGanalysis
of 0.5%, proteinase K (New England BioLabs, Ipswitch, Coding domains were assigned to Clusters of Ortholo-
MA) was added to 100 μg/mL, and RNaseA (MoBio La- gous Groups (COG) classes by using the RPSBLAST
boratories, Carlsbad, CA) was added to 100 μg/mL. The program, version 2.2.26 against the CDD database [29],
mixture was incubated at 37 °C for 1 h. NaCl was added which was downloaded from the NCBI website on April
to0.5M,and0.5mLofphenol:chloroform:isoamylalco- 24th, 2015. The COG annotations were used to generate
hol (25:24:1) was added. The solution was mixed and genome maps using a combination of BLAST Ring
then centrifuged at 13,000×g for 10 min at 4 °C. The Image Generator (BRIG) [30] and CGview [31]. A ran-
aqueousphasewastransferredtoanewtubeand0.6mL dom sample of 100 genomes from the RefSeq database
of isopropyl alcohol was added. The mixture was incu- [32] was downloaded on July 7th, 2016. The random
bated at −20 °C for ~16 h and then centrifuged at genomes were assigned to COG classes using the
13,000×g for 5 min at 4 °C. The pellet was washed with previously described method. To estimate error in the
70% ethanol, centrifuged at 13,000×g for 5 min at 4 °C, sampling, the 100 genomes were sampled with replace-
and finally resuspended in 10 mM Tris buffer pH 8.0 ment 1000 times. Significance was determined through
andstoredat−20°C[17]. Fisher’s exact test, and verified through difference of
Gilmoreetal.BMCGenomics (2017) 18:639 Page4of14
medians analysis [33]. For comparative purposes, the Nucleotidesequenceaccessionnumbers
fraction of genes assigned to a given COG class was The final genome assembly and annotation information
defined as the number of genes assigned to the class is available in the GenBank database with accession
divided by the total number of genes assigned to any numbers PRJNA300714, PRJNA300715, PRJNA300716,
COG class, which accounted for differing genome sizes and PRJNA300717 for Mbac. bryantii, Msar. spelaei,
and unknown gene content. Msph. cuniculi,andMcor.parvum, respectively.
The evolutionary history of the methanogen isolates
was inferred using the Neighbor-Joining method [34]. Results and discussion
A sum of branch length of 1.91983843 was used. A Genomesequencingandassembly
bootstrap test with 1000 replicates was used to com- In order to compare methanogen metabolism and meth-
pute the final tree [35]. The evolutionary distances anogenesis across archaeal genera, four methanogens
were computed using the Maximum Composite Likeli- (Mbac. bryantii, Msar. spelaei, Msph. cuniculi, and
hood method [36] and are in the units of the number Mcor. parvum) were selected for genomic sequencing
of base substitutions per site. The analysis involved 91 and analysis. Collectively, the methanogens chosen rep-
nucleotide sequences. All positions with less than 95% resent a diverse set of isolates from four different genera
site coverage were eliminated. That is, fewer than 5% and all three classes [18], capable of methylotrophic,
alignment gaps, missing data, and ambiguous bases acetoclastic, and hydrogenotrophic methanogenesis, as
were allowed at any position. There were a total of shown in Table 1 [44–48]. Draft genomes for Mbac.
1246 positions in the final dataset. Evolutionary ana- bryantii, Msar. spelaei, Msph. cuniculi, and Mcor. par-
lyses were conducted in MEGA7 [37]. Redundant and vum were assembled and annotated through the Depart-
unnecessary branches were trimmed for a more easily ment of Energy Systems Biology Knowledgebase
viewed tree and formatted using the Interactive Tree platform (KBase), a summary of which can be seen in
of Life [38]. Table 2. The size of each genome was as expected for
methanogens of their respective genera, with Msar.
Additionalannotationandanalyses spelaei the largest (5.1 Mb), Mcor. parvum and Msph.
Possible modes of methanogenesis as described on cuniculi the smallest (1.7 and 1.9 Mb, respectively), and
metacyc.com [39] and in Blaut [40] were investigated Mbac. bryantii in between (3.5 Mb). Mcor. parvum had
for each genome. Using the R programming language, the highest GC content (51%), while the other three
an algorithm was developed to search the functional were more AT-rich (~30% GC content). The GC bias is
annotations provided by the KBase pipeline for genes genome-wide, extending to both coding and intergenic
that encode enzymes involved in methanogenesis. regions.Ananalysisofcodonfrequency(Additionalfile1:
Methanogenicgeneswere identified bysearching for ei- TableS1)reflectedtheGCbias,whereMcor.parvumhad
therEnzymeCommission(EC) numbers[41]orspecific theloweststopcodonfrequencyofTAA(42%),compared
enzyme names and functional descriptions. Hydroge- to 67%, 73%, and 54% for Mbac. bryantii, Msph. cuniculi,
nases were identified by searching the KBase annota- andMsar.spelaei,respectively[49].
tions specifically for Eha, Ehb, Mvh, Frh, Ech,Vho, and
cytochrome b enzymes. Since these enzymes are multi- AnalysisofthedistributionofCOGproteinsrevealsa
subunit complexes, the components of each hydrogen- shiftingenomiccompositiontowardsenergy
ase were identified. BLAST-based comparisons were conservation
made using BLASTP version 2.2.30+ against the gen- The predicted genes were assigned to COG classes
ome indicated, each downloaded from the RefSeq data- using RPSBLAST. A summary of the number of genes
base [32] on August 1st, 2016. assigned to each COG class can be seen in Fig. 1 and
Transmembrane domains were predicted using a Additional file 1: Table S2. The largest classes of COGs
Transmembrane Hidden Markov Model (TMHMM) represented across the methanogens include Energy
[42]. Transporter analysis was conducted by running Production and Conversion [C], Coenzyme Transport
BLASTP version 2.2.30+ against the transporter classi- and Metabolism [H], and Inorganic Ion Transport and
fication database (TCDB) [43], downloaded on January Metabolism [P]. Analysis of the predicted membrane
13th, 2015. A custom python script was used to select proteins further demonstrated the trend towards energy
the blast hit with the lowest e-value, with the strict re- conservation (Additional file 1: Figure S1). The three
quirement that the blast hit matched at least 70% of the largest classes of membrane proteins and transporters
length of the subject and query. Membrane proteins belonged to energy conservation, ion transport, and
and transporters were assigned classes based on TCDB metal or cofactor acquisition.
number (transporters) or by manual curation of KBase In order to quantitativelycomparethe COG categories
annotations (membrane proteins). of the sequenced methanogens, a random sample of 100
Gilmoreetal.BMCGenomics (2017) 18:639 Page5of14
Table1MethanogensCharacterizedinthisStudy.Displaysasummaryofthefourmethanogenssequenced:Mbac.bryantiiisolation
datafrom[46,47],Msar.spelaeiisolationdatafrom[48],Msph.cuniculiisolationdatafrom[45],Mcor.parvumisolationdatafrom[44]
Methanobacteriumbryantii Methanosphaeracuniculi Methanosarcinaspelaei Methanocorpusculumparvum
CurrentClassification Archaea,Euryarchaeota, Archaea,Euryarchaeota, Archaea, Archaea,
Methanobacteria, Methanobacteria, Euryarchaeota, Euryarchaeota,
Methanobacteriales, Methanobacteriales, Methanomicrobia, Methanomicrobia,
Methanobacteriaceae, Methanobacteriaceae, Methanosarcinales, Methanomicrobiales,
Methanobacterium,bryantii Methanosphaera,cuniculi Methanosarcinaceae, Methanocorpusculaceae,
Methanosarcina,spelaei Methanocorpusculum,parvum
MethanogenClass[18] ClassI ClassI ClassIII ClassII
Gramstain Variable Gram-positive Gram-negative Gram-negative
Cellshape Rod Coccus Sarcina-likecoccus Irregularcoccus
Motility Non-motile Non-motile Non-motile Weaklymotilebysingle
flagellum
Sporulation Nonsporulating Nonsporulating Nonsporulating Nonsporulating
Optimaltemperature 37–45°C 35–40°C 33°C 15–45°C
range
Oxygenrequirement Strictlyanaerobic Strictlyanaerobic Strictlyanaerobic Strictlyanaerobic
CarbonAssimilation CO Autotrophy RequiresAcetate CO Autotrophy RequiresAcetateor
2 2
YeastExtract
Energysource H/CO H/methanol H/CO,acetate,methanol, H/CO,formate,
2 2 2 2 2 2 2
monomethylamine, 2-propanol/CO
2
dimethylamine,trimethylamine
Biosafetylevel BSL1 BSL1 BSL1 BSL1
Isolationsource Syntrophiccultureisolated Intestinaltractofa Subsurfacesulfurouslake Anaerobicsourwheydigester
fromsewagesludge rabbit inoculatedwithsewagesludge
Table2GenomeSequencingStatisticsforStrainsinthisStudy
Organism Mbac.bryantii Msph.cuniculi Msar.spelaei Mcor.parvum
Numberofreads 7,552,398 6,883,020 9,312,063 4,790,894
Readlength(bp) 75 75 75 75
Coverage 166 102 137 213
Totallength(Mb) 3.5 1.9 5.1 1.7
Largestscaffold(kb) 1030 188 103 255
Numberofscaffolds(>1000bp) 14 29 293 17
GC% 33 28 39 51
N50(kb) 764 138 37 116
L50 2 6 48 5
Numberofuniquegenespredicted 3526 1658 5913 1844
GenesassignedtoCOGs 2243 1139 2814 1322
Geneswithsignalpeptides 300 104 596 167
Genesassignedtotransporterclassification(TCDB) 416 170 556 247
Genesencodingtransmembranehelices 934 331 1367 376
#predictedgenes(n≥0bp) 3526 1658 5913 1844
#predictedgenes(n≥300bp) 2831 1451 3584 1522
#predictedgenes(n≥1500bp) 335 228 500 203
#predictedgenes(n≥3000bp) 38 41 54 12
Totallength(n≥1000bp) 3,463,789 1,930,335 5,029,712 1,709,622
Gilmoreetal.BMCGenomics (2017) 18:639 Page6of14
Fig.1Annotatedgenomemapshighlightkeymethanogenesisgenes.CirculargenomemapsforMsar.spelaei,Mbac.bryantii,Mcor.parvum,and
Msph.cuniculi.Concentricringsrepresentthefollowingfromoutermosttoinnermost:1)AssembledScaffoldboundaries.2)Genesresponsible
formethanogenesisontheforwardstrand.3)PredictedORFs,coloredbyCOGclass,andpredictedtRNAandrRNAontheforwardstrand.4)
PredictedORFs,coloredbyCOGclass,andpredictedtRNAandrRNAonthereversestrand.5)Genesresponsibleformethanogenesisonthe
reversestrand.6)GCContent.7)GCSkew
genomes from the NCBIRefSeq Database was annotated the enrichment in Energy Production and Conversion
ina similarmethod.As showninAdditional file1:Figure genes, hasbeen noted previously [19, 50–52]. Assuch,a
S2, the relative distribution of COG proteins from the focus on the energy conservation mechanisms of the
RefSeq database achieved stability by 100 genomes sam- methanogensstudiedherewascrucialtoaccuraterecon-
pled. Additional file 1: Figure S3 shows a bootstrap ana- structionsofmetabolism.
lysisofthesampling,whichindicatesastandarderrorof
less than 0.3% for each COG category distribution. As BLAST-basedcomparisontocloselyrelatedspecies
showninFig.2,allmethanogenssequencedweresignifi- reflectssimilaritieswithinmethanogengenera
cantly enriched for genes in Energy Production and The evolutionary relationship of the methanogens was
Conversion (C), Coenzyme Transport and Metabolism determined by placing them into a phylogenetic tree as
(H), and Translation, Ribosomal Structure and Biogen- shown in Fig. 3. Each methanogen clustered closely
esis (J). Furthermore, they had significantly fewer genes with other methanogens from its assigned genus, with
assigned to Carbohydrate Transport and Metabolism the Class I methanogens (Mbac. bryantii and Msph.
(G), Lipid Transport and Metabolism (I), and Secondary cuniculi) distinguished from the Class II methanogen
Metabolites Biosynthesis, Transport, and Catabolism (Mcor. parvum) and the Class III methanogen (Msar.
(Q).Theabnormalgenomecomposition,andinparticular spelaei). Methanosphaera and Methanocorpusculum
Gilmoreetal.BMCGenomics (2017) 18:639 Page7of14
Fig.2DistributionofCOGgenesincategoriessignificantlydifferentthantheRefSeqSample.Thefractionofproteinsineachmethanogen
belongingtoCOGcategoriesshowndiffersignificantlycomparedtotheRefSeqsample.Allfoursequencedmethanogenshavesignificantly
moregenescategorizedasenergyproductionandconversionandcoenzymetransportandmetabolism,whichconfirmstheobservationthat
themethanogenshaveaconfirmedshifttowardsenergyconservationinordertooccupytheirecologicalniche.Significancewasdetermined
throughFisher’sexacttest,andverifiedthroughdifferenceofmediansanalysis[24],*representsp<0.01
Fig.3Sequencedmethanogensclusterwithothermethanogensintheirrespectiveclasses.Thephylogenetictreerepresentstheevolutionary
relationshipofthesequencedmethanogens(highlightedincolorbyClass)comparedtoothercloselyrelatedsequencedarchaea.Asshown,
methanogensinthisstudydisplayedacloseevolutionaryrelationshiptoothermethanogenswithintheirgenera.Withinthistree,the
Methanomicrobiales(representedbyMcor.parvum)arephylogeneticallyclosetotheMethanosarcinales(Msar.spelaei),butmetabolically
moresimilartotheMethanobacteriales(Msph.cuniculiandMbac.bryantii).Bootstrapvaluesareindicatedforeachnode
Gilmoreetal.BMCGenomics (2017) 18:639 Page8of14
are both underrepresented in number of species se- group of proteins annotated in the datasets belonged
quenced to date. As such, Msph. cuniculi and Mcor. to either Type I or Type II restriction-modification
parvum were both compared to the other sequenced system, suggesting slight differences in specificity
species in their genera. A BLAST-based comparison against foreign DNA. Taken together, these observa-
was performed to determine the proteins found in tions suggest that Msph. cuniculi and Mcor. parvum
Msph. cuniculi or Mcor. parvum not found in other se- are metabolically similar to other Methanosphaera and
quenced species. The results of the BLAST-based com- Methanocorpusculum, respectively, and those other se-
parison are shown in Additional files 2 and 3 for quenced isolates represent good templates for meta-
Msph. cuniculi and Mcor. parvum, respectively. The bolic reconstruction.
most striking observation from this comparison is the
number of genes annotated as either “hypothetical pro- Metabolicreconstructionofmethanogenesispathways
tein” or with no annotation. Two hundred fifty three of In order to further investigate proteins involved in en-
the 277 proteins (91.3%) found in Msph. cuniculi, but ergy conservation mechanisms, the specific pathways of
not in Methanosphaera stadtmanae, lacked a descrip- methanogenesis were reconstructed for each sequenced
tive annotation, and the same was true for 230 of the isolate (Fig. 4). Of particular interest were the distinct
293 proteins (78.5%) found in Mcor. parvum, but not energy conservation mechanisms employed by each
Methanocorpusculum bavaricum or Methanocorpuscu- isolate and the predicted substrate utilization routes.
lum labreanum. Of the few proteins containing a de- Together, these two aspects provide the backbone for
scriptive annotation, none of them play a significant metabolic models, which can be used for predicting be-
role in methanogenesis in either isolate. The largest havior in microbial communities.
Fig.4Metabolicreconstructionofmethanogenesisrevealstwomechanismsofenergyconservation.Energyconservationthrough
methanogenesisisdetailedforelectronbifurcation(a)andchemiosmoticcoupling(b).Electronbifurcationconservesenergyin
methanogenesisbytakingtwopairsofelectronsfromtwoseparatehydrogenmoleculesandsplittingthemintoahighenergystate(CO or
2
ferredoxinreduction)andalowenergystate(CoB-S-S-CoMheterodisulfidereduction).GenesforthismechanismwerefoundinMsph.cuniculi,
Mbac.bryantii,andMcor.parvum,althoughthecoupledhydrogenasewasnotidentifiedinMcor.parvum.Msar.spelaeiutilizeschemiosmotic
couplingforenergyconservation,whereNa+orH+transportintothecellislinkedtoH oxidation,andtransportoutofthecellislinkedto
2
methyltransferaseandheterodisulfidereductaseactivity,establishinganetoutwardgradientforATPproduction.Thehydrogenasedepicted
in(a)representstheMvhhydrogenase(Msph.cuniculiandMbac.bryantii),andthehydrogenasein(b)representstheEha(Mbac.bryantiiand
Mcor.parvum),Ehb(Mbac.bryantiiandMsph.cuniculi),orEchhydrogenase(Msar.spelaeiandMcor.parvum).Methylotrophicmethanogenesis
pathwaysaredisplayedin(c).MethanolutilizationpathwaysarefoundinMsph.cuniculi,Msar.spelaei,andMbac.bryantii.Acetateandmethylamine
utilizationpathwaysarefoundonlyinMsar.spelaei.Theacetyl-CoAsynthasecomplexisfoundinMbac.bryantiiandMsar.spelaei,allowing
themtofixCO
2
Gilmoreetal.BMCGenomics (2017) 18:639 Page9of14
Methanobacteriumbryantii CDS.6040–4)alongwithaformatetransporter(CDS.6039),
As a member of the Class I methanogens, Mbac. bryantii suggesting the possibility of growth on formate, although
was expected to utilize a cytoplasmic heterodisulfide re- this phenotype in Mbac. bryantii has not been reported
ductase complex (Hdr, CDS.4429–30), coupled to a non- previously [58]. It also contains the complete set of genes
F reducing hydrogenase (Mvh, CDS.4299–302) [53]. necessaryformethanogenesisfrommethanol/H ,including
420 2
These two enzymes couple with the formylmethanofuran a methyl transferase (CDS.5844), corrinoid protein
dehydrogenase complex (Fwd, CDS.4483–89), which al- (CDS.5843), and corrinoid activating protein (CDS.5842).
lows for electron bifurcation between the Mvh and Hdr, Asulfitereductasegeneisalsopresent(CDS.5502),which
wheretwopairsofelectronsaresplitbetween athermo- likely explains the previously reported observation that it
dynamically unfavorable reaction (CO reduction to for- resistssulfiteinhibition[59].Finally,Mbac.bryantii’sgen-
2
mylmethanofuran or reduction of ferredoxin) and a ome possesses an alcohol dehydrogenase (CDS.3886) and
thermodynamically favorable reaction (heterodisulfide NADP-dependent F reductase (CDS.3694), suggesting
420
reduction) [51]. As shown in Fig. 4 and Table 3, Mbac. the possibility of growth on isopropanol or isobutanol,
bryantii contains these components and likely catalyzes sometimes seen in other Methanobacterium species [46].
methanogenesis through electron bifurcation. Mbac. However, attempts at growth on formate, methanol, and
bryantii has all the genes necessary for its previously re- isopropanol were unsuccessful under conditions tested in
ported phenotypic growth on H /CO and contains this study, suggesting that the genes responsible are in-
2 2,
genes for an acetyl-CoA synthase complex (CDS.3915– active or the proper conditions for gene activation were
3920),allowingittofixCO forcentralmetabolism.The unmet. The presence of the noted genes are important,
2
Fwd complex, which requires either tungsten or molyb- however, when consideringhow Mbac. bryantii might act
denum as a cofactor [54], was found in close proximity in co-culture or in its native environment. For example,
to genes related to molybdenum transport (CDS.4477– both formate and hydrogen are commonly used to
80), molybdopterin synthesis (CDS.4490), and tungsten transfer electrons in microbial consortia [60], so it is
transport (CDS.4475), as seen in other methanogens possible that formate metabolism is triggered by low
[55]. The clustering of genes suggests that they may be partial pressures of H .
2
co-regulated to aid in complex formation. Two sets of
energy-conserving hydrogenases, Eha (CDS.5584–90) Methanosphaeracuniculi
and Ehb (CDS.6156–70), were found in the genome of A genomic analysis of Msph. cuniculi revealed that it is
Methanobacterium bryantii. Eha and Ehb have been as- metabolically similar to the other Methanosphaera
sociatedwithanaplerosis[56]andCO assimilation[57], sequenced, Msph. stadtmanae. Msph. stadtmanae is
2
respectively, during hydrogenotrophic methanogenesis known to have a restrictive metabolism, requiring both
inMethanococcusmaripaludis. methanol and H for growth [61]. Like Msph. stadtma-
2
Mbac. bryantii possesses several genes for metabolic nae, Msph. cuniculi contains the genes encoding the
substrate utilization not previously observed from Mbac. non-F -reducing hydrogenase, Mvh (CDS.1980–2)
420
bryantii in isolation. It contains several copies of the [61]. As with Mbac. bryantii, this hydrogenase likely
formatedehydrogenasegenes(CDS.3844–7,CDS.333–4, couples with the heterodisulfide reductase complex,
Hdr (CDS.2562–4) [51]. The heterodisulfide reduction
Table3HydrogenaseComponentsoftheMethanogensin is likely coupled to oxidation of two H molecules and
2
thisStudy ferredoxin reduction. The ferredoxin is then either used
Species Hydrogenases Components for anabolic processes or oxidized by Ehb (CDS.2129–
Mbac.bryantii Eha B,C,F,G,H,J,M,N,O,P,P3,R,gene2(V) 49) in an energy-conserving mechanism producing H2
Ehb D,E,F,H,I,K,L,M,N,O,Q and transporting ions out of the cell, as proposed re-
Mvh A,D,G
cently by Thauer et al. [5].
Frh A,B,G
The utilization of methanol likely occurs through
Msph.cuniculi Ehb A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,Q
the action of three proteins: a corrinoid protein
Mvh A,D,G
Frh A,B,G (CDS.2507), a methyltransferase corrinoid activation
Msar.spelaei Ech A,B,C,D,E,F protein (CDS.2506), and a methanol methyltransferase
Cytochromeb I,II protein (CDS.2508) [62]. These proteins transfer a me-
Vho A,C,G
thyl group from methanol, to the corrinoid protein,
Frh A,B,G
and finally to methyl-coenzyme M. Beyond methanol
Mcor.parvum Ech A,B,C,D,E,F metabolism, Msph. cuniculi contains all the genes ne-
Mbh A,B,C,D,E,F,G,H,I,J,K,L,M,N
Eha B,C,D,E,F,G,H,J,M,N,O cessary for reduction of CO2 for methanogenesis, as
Ehb Q shown in Fig. 4a. However, previous characterization
Frh A,B,G
has concluded that it requires methanol for growth
Gilmoreetal.BMCGenomics (2017) 18:639 Page10of14
[45]. The lack of growth is most likely explained by andpyrrolysyl-tRNAsynthetase,allofwhichwerefound
the requirement of molybdopterin by the formyl- in the same region of the genome as the monomethyla-
methanofuran dehydrogenase, and the lack of enzymes mine methyltransferase complex (CDS.7548–51) [67].
capable of producing molybdopterin in Msph. cuniculi. Finally, an acetate kinase (CDS.6562) and phosphotran-
However, it is also possible that genes are not functional sacetylase (CDS.6561) were found, explaining growth on
in vivo. Like Mbac. bryantii, it is important to note the acetate[68].
potentiallynovelsubstrateroutes,asMsph.cuniculimay
be capable of utilizing CO and H in a native environ- Methanocorpusculumparvum
2 2
mentorinconditionsnotseeninisolation. As shown in Fig. 4, Mcor. parvum has all the genes ne-
cessary for methanogenesis with H /CO and formate as
2 2
Methanosarcinaspelaei substrates. Mcor. parvum lacks the acetyl-CoA synthase
Insteadofutilizingelectronbifurcationtocoupletheun- (Fig. 4c), uncoupling methanogenesis from CO fixation
2
favorable reduction of CO to the favorable reduction of into acetyl-CoA, which explains its phenotypic require-
2
the CoM-S-S-CoB heterodisulfide, H -utilizing Metha- ment for acetate or yeast extract for carbon assimilation
2
nosarcina couple each step to ion transport, termed [44]. It contains the most diverse set of hydrogenases of
chemiosmotic coupling, shown in Fig. 4b [5]. The favor- the methanogensstudied here, asshowninTable 3.Two
able step of heterodisulfide reduction is coupled to separate gene clusters were found within the genome
transportofH+orNa+outofthe cell,establishingagra- containingtheformylmethanofurandehydrogenase com-
dient [5]. This step is performed by a methanophenazine plex, Fwd. In one case, the Fwd genes (CDS.2339–40)
reducing hydrogenase (Vho, CDS.8117–25) / heterodi- are clustered with an energy conserving hydrogenase,
sulfide reductase (Hdr, CDS.5944–6) complex [5]. Part Eha (CDS.2325–37). Several of the Eha components are
of this gradient is then used by the energy conserving predicted to be membrane associated (Additional file 4)
hydrogenase (Ech, CDS.6169–74) to reduce ferredoxin with a predicted function of coupling Fwd activity to
with H [5]. Finally, the ferredoxin is utilized to reduce transmembrane sodium or proton transport according
2
CO by the formylmethanofuran dehydrogenase com- to the transporter analysis (Additional file 5) [5]. A
2
plex (Fmd, CDS.9330–4) [5]. As shown in Fig. 4b and second cluster contains genes for the Fwd complex
Table 3, Msar. spelaei’s genome contains all of the genes (CDS.3326–9) along with the heterodisulfide reductase
required for this process. Additionally, the membrane (HdrABC, CDS.3322–4), including two HdrA subunits
protein and transporter analysis (Additional files 4 and and one each of the HdrB and HdrC subunits. None of
5) revealed that the Ech and Hdr complexes are pre- the associated subunits contain predicted transmem-
dicted to be membrane-bound and catalytically active brane domains, suggesting that the Hdr complex is cyto-
transporters, which has been experimentally validated solic, as previously determined [5]. Unlike members of
in other isolates [63, 64]. The genome for Msar. spelaei the Class I methanogens or the Methanocellales species
contains the F -reducing hydrogenase (FrhABG, [69], the complementary [NiFe] hydrogenase (MvhAG)
420
CDS.9055–8), used for CO reduction. It also contains was not identified in the Mcor. parvum genome, but the
2
a methanophenazine-dependent F H dehydrogenase MvhD subunit was present. The FrhABG hydrogenase
420 2
complex (Fpo, CDS.7883–95), which reduces methano- (CDS.2912–15) was identified as well. The presence of
phenazine with F H while transporting two protons two separate gene clusters coding for the Fwd complex
420 2
across the membrane. This allows for the interconver- suggest that they may be expressed in response to differ-
sion of reduced cofactors (methanophenazine and F ) ent growth conditions, or play different roles in meth-
420
while conserving energy, which is important for methy- anogenesis. The Eha hydrogenase has been shown to be
lotrophic methanogenesis [65]. critical for replenishing intermediates of methanogenesis
Genome sequencing of Msar. spelaei presents an ex- lost to leaky electron bifurcation, dilution due to growth,
planation for its observed methanogenesis on methanol, or biosynthesis [56]. It was also proposed that a separate
acetate, and methylamines, as shown in Fig. 4c. Methyl- ferredoxin was utilized in order to separate electron
transferases specific to methanol (CDS.8578–80), mono- pools for anabolism and replenishment of intermediates
methylamine(CDS.10223–4),dimethylamine(CDS.10167– [56]. The weak similarity of the two Fwd complexes in
8),andtrimethylamine(CDS.9576,9578)wereallidentified Mcor. parvum (~36–44% amino acid identity of each
within the genome. Each methylamine methyltransferase component), could represent differences in specificity of
includes an amber stop codon that putatively encodes a the ferredoxin used – one specific for ferredoxin re-
pyrrolysineresidue,whichhasbeenshowntobecriticalfor duced by Eha, and one specific for ferredoxin reduced
activityinMethanosarcinaacetivorans[66].Supportingthe by a putative electron bifurcating Hdr complex. Separate
synthesis of this non-canonical amino acid are proline ferredoxin proteins are encoded in each of the Fwd gene
reductase, pyrrolysine synthetase, proline 2-methylase, clusters, further supporting this hypothesis. The exact