New Insights into the Ecology and Physiology of Methanomassiliicoccales from Terrestrial and Aquatic Environments

Members of the archaeal order Methanomassiliicoccales are methanogens mainly associated with animal digestive tracts. However, environmental members remain poorly characterized as no representatives not associated with a host have been cultivated so far. In this study, metabarcoding screening combined with quantitative PCR analyses on a collection of diverse non-host-associated environmental samples revealed that Methanomassiliicoccales were very scarce in most terrestrial and aquatic ecosystems. Relative abundance of Methanomassiliicoccales and substrates/products of methanogenesis were monitored during incubation of environmental slurries. A sediment slurry enriched in Methanomassiliicoccales was obtained from a freshwater sample. It allowed the reconstruction of a high-quality metagenome-assembled genome (MAG) corresponding to a new candidate species, for which we propose the name of Candidatus ‘Methanomassiliicoccus armoricus MXMAG1’. Comparison of the annotated genome of MXMAG1 with the published genomes and MAGs from Methanomassiliicoccales belonging to the 2 known clades (‘free-living’/non-host-associated environmental clade and ‘host-associated’/digestive clade) allowed us to explore the putative physiological traits of Candidatus ‘M. armoricus MXMAG1’. As expected, Ca. ‘Methanomassiliicoccus armoricus MXMAG1’ had the genetic potential to produce methane by reduction of methyl compounds and dihydrogen oxidation. This MAG encodes for several putative physiological and stress response adaptations, including biosynthesis of trehalose (osmotic and temperature regulations), agmatine production (pH regulation), and arsenic detoxication, by reduction and excretion of arsenite, a mechanism that was only present in the ‘free-living’ clade. An analysis of co-occurrence networks carried out on environmental samples and slurries also showed that Methanomassiliicoccales detected in terrestrial and aquatic ecosystems were strongly associated with acetate and dihydrogen producing bacteria commonly found in digestive habitats and which have been reported to form syntrophic relationships with methanogens.


. Total DNA extraction
In order to search for natural samples containing Methanomassiliicoccales, an initial large-scale molecular screening was carried out on a wide range of anoxic samples (86 in total).
DNA extractions were carried out in triplicate with a protocol adapted to each type of matrix. Two types of DNA extraction procedures have been implemented depending on the nature of the matrices: (i) Nucleic acids were extracted from sediments, mud and peatland soils by combining a mechanical lysis step to a commercial DNA extraction kit for soil type matrices. The first extraction step consisted of mechanical cell lysis of 0.5 g of sample by bead-beating (5500 rpm, 30 s), using the Precellys 24 ® tissue homogenizer (Bertin instruments, France). Then, the FastDNA™ SPIN Kit for Soil (MP Biomedicals, Thermo Fisher Scientific) was used following manufacturer's instructions. In order to optimize the precipitation of small amounts of DNA and to improve extraction yields, 10 µl of linear acrylamide (5 mg ml -1 , Ambion, Thermo Fisher Scientific) was added to the sample. After a chemical lysis step, DNA was purified on a silica column and eluted in 50 µl of DES buffer. (ii) Nucleic acids were extracted from liquid samples (brines, lake waters, geothermal hot spring fluids) using a phenol/chloroform/isoamyl alcohol (25/24/1; pH 8.0; Sigma-Aldrich) (PCI) DNA extraction protocol. Fifteen ml of sample were first centrifuged (20 min, 13 000 rpm, 4 °C) to pellet the cells. Cell pellets were then resuspended in 1 ml TE-Na-1X lysis buffer, with 10 µl of linear acrylamide (5 mg ml -1 ). Cell lysis and deproteinization were performed by addition of 100 µl sarkosyl (10% w/v), 100 µl sodium dodecyl sulfate (SDS; 10% w/v) and 20 µl proteinase K (20 mg ml -1 , Thermo Fisher Scientific, Illkirch, France), followed by a 3h incubation at 55 °C in a water bath. After incubation, 1 volume of PCI was added to 1 volume of lysis solution, solutions were mixed by gentle inversion and then centrifuged (15 min, 13 000 rpm, 4 °C) to extract the nucleic acids. After centrifugation, the upper aqueous phase was collected and extracted with chloroform. After centrifugation (15 min, 13 000 rpm, 4 °C), the upper aqueous phase was collected and nucleic acid precipitation was carried out by addition of 40 µl sodium acetate (3 M, pH 5.2) and 0.8 volume of isopropanol. In order to favour nucleic acids precipitation, each tube was placed 30 min at -80 °C. After thawing at room temperature, the solutions were centrifuged (20 min, 13 000 rmp, 4 °C) to collect DNA pellets and these were washed with 70% (v/v) ethanol. Finally, DNA were resuspended in 50 µl TE-1X buffer.

Text S2. Origin and detailed description of environmental samples
The eighty-six samples collected from anoxic zones of various ecosystems (marine sediments, deep-sea hydrothermal vents, marine pockmarks, submarine mud volcano, deep-sea hypersaline anoxic brine, peatland soils, lakes, hot springs, river sediments) from 40 worldwide locations, for the pre-screening stage are described in Table A1.
The 22 samples used for an in-depth study were the following ones: (i) five samples, collected in January 2017, were originating from the ombrotrophic Mougau peatland area (Brittany, France). They consisted of four peatland soil samples (MOUG1-MOUG4) collected in January 2017, in the peat surface layer (0-10 cm), and one sediment sample from a freshwater stream (MOUG5) collected in the upper sediment (5-10 cm); (ii) four samples were collected from the center of the deep-sea methane-emitting Håkon Mosby mud volcano, located at 1280 m water depth, in the Barents Sea (Arctic Ocean), during the 2009 ARK-XXIV/2 (PS74) oceanographic cruise. Sedimentary horizons selected after the pre-screening stage were the 0-1 cm (HM1), 1-6 cm (HM2) and 6-11 cm (HM3) sediments from the PS74/168 core, and the 1-6 cm (HM4) surface layer sediments from the PS74/189 core; (iii) three anoxic samples were collected in December 2018, in the water column of the meromictic crater Lake Pavin (Auvergne, France) at 60, 70 and 80 m water depth corresponding, respectively, to the mixolimnion-chemocline interface (PAV60), the chemocline-subchemocline interface (PAV70) and the subchemocline water (PAV80) [4]; (iv) one geothermal hot spring sediment sample was collected, in November 2013, from the Xiada reservoir of the Xiamen botanical garden (China) (XIA); (v) one geothermal spring sediment sample was taken from an ephemeral hot spring located in the tidal sway zone on the beach called Feu de Joie, at the Kerguelen islands (Kerguelen archipelago, France, Indian Ocean) (KERG), in December 2016; (vi) two samples were originating from the athalassic (MgCl2-rich) deep-sea hypersaline anoxic brine (DHAB) Kryos, in the Mediterranean Sea. They were respectively collected from layers at 150 g l -1 salts (mainly MgCl2) (KRY150) and 238 g l -1 (mainly MgCl2) (KRY238) [5], in May 2016; (vii) three sediments samples were collected in a pockmark field during the 2015 PAMELA-MOZ4 oceanographic survey in the Mozambique Channel (off the Madagascar island). Those samples consisted in the first 2-4, 4-6 and 6-11 cm layers of the CSF4 sediment cores referenced as MOZ1, MOZ2 and MOZ3, respectively; (viii) one surficial marine sediment sample (15-20 cm) was collected from the Dourduff-en-Mer marine bay (DOUR), in February 2017, on the coast of Brittany (English Channel, France); (ix) one deep-sea sediment sample was collected from the South West Indian Ocean (COMRA) during the April 2015 COMRA oceanographic cruise; (x) finally, a freshwater sediment sample was taken from a stream (0-10 cm) feeding the coastal river Penfeld in February 2017 (PENF) (Brittany, France). Given that the ultimate goal of this study was to grow Methanomassiliicoccales, all the investigated samples were stored at 4 °C.

Text S3. Total DNA extraction for metabarcoding screening
The standardized DNA extraction protocol which was used on the 22 samples selected for an indepth characterization combined mechanical and chemical lysis. It was applied in triplicate to each sample (+ one negative control). In summary, each DNA extraction was performed onto 0.5 g of environmental matrix. A mechanical lysis of the sample was first carried out by bead-beating (5000 rpm, 15 s), using the Precellys 24 ® tissue homogenizer (Bertin instruments, France), in a FastDNA ® Lysing Matrix tube E (FastDNA ® SPIN kit for soil, MP Biomedicals) containing 1 ml of TE-Na-1× lysis buffer (100 mM Tris, 50 mM EDTA, 100 mM NaCl, pH 8.0) in addition to the sample. Then, the lysate supernatants were subjected to a DNA extraction protocol combining chemical and enzymatic lysis procedures, and a classical phenol/chloroform/isoamyl alcohol (PCI: 25/24/1; pH 8.0) extraction of nucleic acids, as described elsewhere [6]. Nucleic acids were then precipitated with 0.8 volumes of icecold isopropanol, in the presence of 400 µl Na-acetate (3M) and 5 µl linear acrylamide (Ambion, Applied Biosystems). Elution was performed in 40 µl EB buffer (Qiagen) before pooling. DNA extractions from enrichment cultures were performed as described above, except that the mechanical lysis stage has been omitted. For these extracts, elution was performed in 100 µl EB buffer.
Nucleic acid solution quality was determined using the NanoDrop™ 8000 (Thermo Scientific) spectrophotometer. Double-strand DNA concentration was measured using the kit Quantifluor™ dsDNA system (Promega), following the manufacturer's instructions.

Section 4.1. FASTQ generation
Amplicons from bulk samples and enrichment cultures were generated by the company Molecular Research-MrDNA (Shallowater, Texas, USA), using the Illumina MiSeq (2 × 300 bp, pairedend reads) technology, with the kits and according to the manufacturer's instructions. FASTQ were generated using the "Combine FASTA and QUAL into FASTQ (Galaxy v. 1.0.1)" [7] tool on the FROGS pipeline. FASTA and Quality Score files were joined according to a mapping file and the ASCII force quality score to create a single FASTQ block for each read. The quality of raw sequence data were then checked using the "FastQC Read Quality reports" (v. 0.11.4) tool (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Section 4.2. OTU generation and comparison
Once data from enrichment cultures were generated by the company and in order to facilitate OTU comparison between prokaryotic diversity in bulk samples and in enrichment cultures, raw FASTQ files from both sequencing runs were merged into a unique file before processing as described in 2.3. Sequencing and sequence analysis (main manuscript). This resulted on the generation of a unique BIOM file comprising representative 16S rRNA OTU abundance, and taxonomy into each samples.

Section 4.3. Consensus sequences workflow
Extraction of potential reagent contaminants and normalization of representative 16S rRNA gene OTUs were performed using the R studio program (v. 1.3.5001), the phyloseq (v. 1.28.0 R), the decontam (v. 1.4.0) and the metagenomSeq (v. 1.26.3) packages. As the metadata between the bulk samples and in the enrichment cultures were different, two distinct sequence treatments were performed from the previously generated single BIOM file. The BIOM files were converted to a .csv file before being implemented in the R program. Our scripts were based on workflows described elsewhere (https://bioconductor.org/help/coursematerials/2017/BioC2017/Day1/Workshops/Microbiome/MicrobiomeWorkflowII.html and https://benjjneb.github.io/decontam/vignettes/decontam_intro.html, respectively). Workflow used to decontaminate and normalize 16S rRNA gene sequences representatives of the OTUs is listed below. Additional information about the script is preceded by a "#".

Section 4.4. Ecological Network analysis
Beta-diversity between bulks samples was assessed by Network analysis using a Bray-Curtis index, as ecological distance. First, the distance matrix was computed with a script Perl, then the percolation threshold was computed (see methods section) and the Gephi was used to visualize it.

Section 4.5. Phylogenetic analyses of Methanomassiliicoccales OTUs
All 16S rRNA gene sequences representative of Methanomassiliicoccales OTUs, generated from metabarcoding and metagenomic analyses, have been aligned and refined manually with other Methanomassiliicoccales sequences from both Methanomassiliicoccales clades published elsewhere [8][9][10][11][12][13][14][15] using MUSCLE (with 16 iterations maximum) and GBlocks (without contiguous non-conserved positions). Sequences affiliated to various methanogen lineages and sequences belonging to the clades RG8 and Ca. 'Thermoplasmatota' were also added to the dataset, and used as outgroup. The BIONJ method [16] was then applied on sequences with the modifications of Jukes and Cantor and using 1000 resampling bootstrap replicates to generate the final dendrogram. The fatty-acid solution contained, per litre: 25 g valeric acid, 25 g isovaleric acid, 25 g 2methylbutanoic acid and 25 g isobutyric acid. Its pH was adjusted to 7.5 with concentrated NaOH. It was then filter-sterilized (pore-size, ∅ 0.2 µm), placed under a N2 gas phase, and stored at 4 °C.
The porcine hemin solution was prepared as follows, for 100 ml: 50 mg porcine hemin (PanReac AppliChem, Barcelona, Spain) was dissolved into 1 ml NaOH before adjusting the volume to 100 ml with distilled water. The solution was then sterilized by filtration (pore-size, ∅ 0.2 µm) before to be stored at 4 °C, under a N2 atmosphere.
The coenzyme-M solution contained, for 100 ml: 0.164 g Na-2-mercaptoethanesulfonate dissolved into 100 ml sterile water before filtering onto a ∅ 0.2 µm membrane. Headspace was filled with pure N2 gas and stored at 4 °C, in a bottle of brown glass, until use.
The K3 vitamin solution was prepared as follows: 5 mg/ml vitamin K3 was dissolved in 95% (v/v) ethanol and then diluted 1‰ th in distilled water (final concentration 0.05 mg/ml). The solution was then filter-sterilized and stored at 4 °C under N2, in an opaque flask, until use.

Text S6. Quantification of substrates and metabolic products
Headspace CH4, H2, N2 and CO2 were measured using a modified INFICON/Micro GC FUSION Gas Analyzer (INFICON, Basel, Switzerland) fitted with a pressure gauge and two conductivity detectors. Separation was performed using two columns: a 10m molecular sieve column at 65 °C and Ar as a carrier gas; and a 12m RT-Q at 60 °C using He as a carrier gas. Gas concentrations were calculated using the method of Mah and colleagues [20] and the results corrected for the decrease in culture medium volume caused by the withdrawal of samples for analysis by ion chromatography.
Samples for cation, anion and methanol analyses were collected in 2 ml tubes and stored at -20 °C. The dilution of each sample was optimized prior to analysis. Sodium, methylated amines (methylamine, dimethylamine, trimethylamine, choline) and ammonium were analyzed using a Dionex ICS-900 Ion Chromatography System (Dionex, Camberley UK) coupled with a CERS 500 4 mm suppressor and a DS5 conductivity detector (40 °C) and fitted with a RFC-10 Reagent-Free Controller™, an ASDV autosampler, and an IonPac CS16 column maintained at 60 °C in a UltiMate™ 3000 Thermostated Column Compartment (Thermo Scientific, Waltham, MA, USA). The gradient program was as follows: 18 mM methanesulfonic acid (MSA) for 40 min, increase from 10.4 mM MSA min -1 to 70 mM (10 min), and then decrease from 52 mM MSA min -1 to 18 mM (9 min). Acetate, nitrate, thiosulfate and sulfate concentrations were quantified by anion chromatography as described elsewhere [21]. Analysis of methanol concentrations was carried out using an Agilent 6890N Gas Chromatograph (GC) (Agilent, Santa Clara, USA) coupled with a Flame Ionization Detector (FID) instrument. 1 ml of a mixture composed of 400 µl of culture supernatant, 450 µl of ultrapure water and 150 µl of internal standard (Butanol-1) were placed into a 10 ml vial for heating desorption. Then, 500 µl of sample from the gaseous space were injected into the GC-FID using a CTC Combipal autosampler Headspace (Agilent, Santa Clara, USA). Compounds were separated using a Cp-Sil 8CB VARIAN column (30 m, 0.25 mm, i.d; 0.25 µm stationary film thickness) (Agilent, Santa Clara, USA) using the following program: 3 min at 40 °C, followed by an increase of 10 °C min -1 up to 70 °C, held for 15 min; injector and detector temperatures, 250 °C; Nitrogen (2.9 ml min -1 ) was used as a carrier gas.        closest mttB sequences in public databases. This reconstruction was performed using the PhyML algorithm. Two mttB reconstructed sequences were affiliated to Methanomassiliicoccales (bold green) while the other 11 mttB sequences where affiliated to other methanoarchaea (bold purple). Accession numbers are preceded by a "@".