Metabolic versatility of small archaea Micrarchaeota and Parvarchaeota

20 Pages·2017·2.84 MB·English

TheISMEJournal(2018)12:756–775
https://doi.org/10.1038/s41396-017-0002-z
ARTICLE
Metabolic versatility of small archaea Micrarchaeota and
Parvarchaeota
Lin-Xing Chen 1 ● Celia Méndez-García 2,3 ● Nina Dombrowski 4 ● Luis E. Servín-Garcidueñas5 ●
Emiley A. Eloe-Fadrosh6 ● Bao-Zhu Fang1 ● Zhen-Hao Luo1 ● Sha Tan1 ● Xiao-Yang Zhi7 ● Zheng-Shuang Hua1 ●
Esperanza Martinez-Romero8 ● Tanja Woyke6 ● Li-Nan Huang1 ● Jesús Sánchez2 ● Ana Isabel Peláez2 ●
Manuel Ferrer9 ● Brett J. Baker4 ● Wen-Sheng Shu10
Received:5April2017/Revised:26August2017/Accepted:9October2017/Publishedonline:8December2017
©TheAuthor(s)2018.Thisarticleispublishedwithopenaccess
Abstract
Smallacidophilicarchaea belongingtoMicrarchaeota andParvarchaeota phyla areknowntophysically interact with some
Thermoplasmatales members in nature. However, due to a lack of cultivation and limited genomes on hand, their
biodiversity, metabolisms, and physiologies remain largely unresolved. Here, we obtained 39 genomes from acid mine
1234567890 dwroarinlda.ge16(ASMrDR)NaAndgheontespbraisnegdenavniarloynsmesenrtesvaeraoleudndththaet
ElectronicsupplementarymaterialTheonlineversionofthisarticle Parvarchaeota were only detected in AMD and hot spring
(https://doi.org/10.1038/s41396-017-0002-z)containssupplementary
habitats, while Micrarchaeota were also detected in others
material,whichisavailabletoauthorizedusers.
including soil, peat, hypersaline mat, and freshwater,
* BrettJ.Baker suggestingaconsiderablehigherdiversityandbroaderthan
[email protected] expected habitat distribution for this phylum. Despite their
* Wen-ShengShu small genomes (0.64–1.08Mb), these archaea may con-
[email protected]
tributetocarbonandnitrogencyclingbydegradingmultiple
saccharides and proteins, and produce ATP via aerobic
1 StateKeyLaboratoryofBiocontrol,GuangdongKeyLaboratory respiration and fermentation. Additionally, we identified
ofPlantResources,CollegeofEcologyandEvolution,SunYat-
several syntenic genes with homology to those involved in
SenUniversity,Guangzhou510275,People’sRepublicofChina
iron oxidation in six Parvarchaeota genomes, suggesting
2 DepartamentodeBiologíaFuncional-IUBA,Universidadde
their potential role in iron cycling. However, both phyla
Oviedo,Oviedo,Spain
lackbiosyntheticpathwaysforaminoacidsandnucleotides,
3 CarlR.WoeseInstituteforGenomicBiology,Universityof
suggesting that they likely scavenge these biomolecules
IllinoisatUrbana-Champaign,Urbana,USA
from the environment and/or other community members.
4 DepartmentofMarineScience,UniversityofTexasAustin, Moreover,low-oxygenenrichmentsinlaboratoryconfirmed
MarineScienceInstitute,PortAransas,TX78373,USA
our speculationthat both phyla are microaerobic/anaerobic,
5 LaboratoryofMicrobiomics,NationalSchoolofHigherStudies based on several specific genes identified in them.
Morelia,NationalUniversityofMexico,Morelia,Michoacan
58190,Mexico Furthermore, phylogenetic analyses provide insights into
the close evolutionary history of energy related functional-
6 DepartmentofEnergyJointGenomeInstitute,WalnutCreek,CA
94598,USA ities between both phyla with Thermoplasmatales. These
results expand our understanding of these elusive archaea
7 YunnanInstituteofMicrobiology,YunnanUniversity,
Kunming650091,People’sRepublicofChina byrevealingtheirinvolvementincarbon,nitrogen,andiron
cycling, and suggest their potential interactions with
8 DepartmentofEcologicalGenomics,CenterforGenomic
Sciences,NationalUniversityofMexico,Cuernavaca, Thermoplasmatales on genomic scale.
Morelos62210,Mexico
9 InstituteofCatalysis,ConsejoSuperiordeInvestigaciones
Científicas(CSIC),Madrid,Spain
10 SchoolofLifeSciences,SouthChinaNormalUniversity,
Guangzhou510631,People’sRepublicofChina
MetabolicversatilityofARMAN 757
Introduction anaerobic, thus enrichments with inoculum from AMD
systems were performed for these elusive archaea (Sup-
Archaea constitute a considerable portion of microbial plementary Fig. 1).
diversity,andplaysignificantrolesinmanybiogeochemical
cyclesonEarth[1].However,comparedtobacteriatheyare
much less understood and fewer genomes have been Materials and methods
sequenced [2]. The recently delineated superphylum
DPANN includes several phyla of archaea with small cell ARMAN genomes and related metagenomes in
and genome sizes and limited metabolic capabilities [3–6]. public database
To date, 48 DPANN draft genomes are available (Supple-
mentary Table 1; see references therein) and only two In NCBI database, there are two Micrarchaeota and two
symbiotic Nanoarchaeota co-cultures have been obtained Parvarchaeota genomes reported from Iron Mountain:
[5, 7]. ARMAN-1 (NCBI accession number, PRJNA349044),
TwoDPANNphyla,MicrarchaeotaandParvarchaeota, ARMAN-2 (PRJNA38565), ARMAN-4 (PRJNA38567),
referred to as Archaeal Richmond Mine Acidophilic and ARMAN-5 (PRJNA38569), which were included for
Nanoorganisms (ARMAN), were first reported in acid analyses in this study.
mine drainage (AMD) biofilms of Iron mountain (Rich- A metagenomic data set sampled from the Fankou mine
mond, CA, USA) and are among the smallest micro- tailings AMD outflow in 2010 (sample abbreviation:
organisms described to date [8, 9]. The AMD biofilms in FK_AMD_2010), which was published previously [15],
Iron Mountain have been comprehensively studied for was included in this study to retrieve ARMAN genomes
microbial ecology and evolution [10]. Four genomes of from newly assembled contigs. The ultra-small cells were
ARMAN have been obtained from this site, including collected by filtering the AMD outflow through 0.22μm
ARMAN-1 and ARMAN-2 from Micrarchaeota, and pore size filters, which were preserved for DNA extraction
ARMAN-4 and ARMAN-5 from Parvarchaeota. The [15].
metabolicfunctionsofARMAN-2,-4and-5intheAMD The available ARMAN genomes were also used as
biofilms have been speculated based on metaproteomic references to search for additional ARMAN sequences in
analyses [11], while ARMAN-1 was published recently publiclyavailable contigs/scaffolds of >3000 metagenomic
onlytoreportitsCRISPR–Cassystem[12].Interestingly, data sets deposited in the IMG/M system (alignment
ARMAN cells were observed having interactions with length≥250bp, similarity≥75%). This analysis retrieved
Thermoplasmatales cells via pili-like structures [11], and oneadditionalmetagenomicdatasetfromObsidianPoolin
this phenomenon was further documented using cryo- the Yellowstone National Park (SRA accession,
genic transmission electron microscope technology [13], SRP099390; IMG genome ID, 3300002966) with
while the ecological significance of such interactions ARMAN-related genomic sequences (sample abbreviation:
remains unclear. Moreover, ARMAN-specific PCR pri- YNP_OP). Obsidian Pool (N44.36°, W110.26°) was
mersandmetagenomicshaverevealedtheiroccurrencein reported as a mesothermal hot spring (56°C; Supplemen-
many other AMD-related environments [14–17], indi- tary Table 2) [18].
cating wide distributions of related microorganisms in
nature.
Study sites, sampling, DNA extraction, and
Despite these investigations described above, we know
sequencing
littleabouttheirbiodiversity,environmentaldistribution(in
other acidic and non-acidic environments), physiologies,
Metagenomic data sets from two AMD and two hot spring
androlesinbiogeochemicalcycling.Toaddressthesegaps,
environments were collected and analyzed in this study
we assembled and binned 39 new genomes from metage-
(SupplementaryFig.2andSupplementaryTable2);thesite
nomic datasets obtained from two AMD and three hot
descriptions and experimental procedures (e.g., sampling,
spring related environments around the world, and the
DNA extraction and sequencing) are detailed below.
environmental distribution of Micrarchaeota and Parvarch-
aeota taxa were evaluated by analyzing 16S rRNA gene (1) Fankou AMD outflow (FK_AMD_2014). Fankou
sequencesfromthosenewgenomesandNCBIandIMG/M mine is located in north Guangdong province of
databases. Metabolic potentials of Micrarchaeota and Par- China (N25.05°, E113.67°), and the microbial ecol-
varchaeotawerepredictedbasedonfunctionalannotationof ogy of this site was deeply studied previously, to
theirgenes,torevealtheirmetabolicfunctionsandpotential reveal the microbial diversity and composition at
roles in nature. Additionally, genomic information likely spatial and temporal scales of the mine tailings [19],
suggested that ARMAN spp. are microaerobic and/or the shifts of microbial composition and function
758 L-XChenetal.
across the tailings acidification processes [20], and sequenced on an Illumina HiSeq2000 platform with
also the metatranscriptomic activities in the related a PE 100bp kit.
AMD systems [14, 15]. Microorganisms in AMD (4) Tengchong geothermal area (TC_Endo). The Teng-
outflow were collected in September, 2014, using chong geothermal area is located in the south of
filters with different pore sizes (i.e., 1.6, 0.45, and Yunnan province of China (N24.95°, E98.44°),
0.22μm). A total of 20L AMD outflow was representing the largest geothermal area in China.
successively filtered through the abovementioned Endolithic communities were detected in the surface
three filters, then the 0.22μm filters with ultra-small ofarocknearahotspringfumarole(referredtoasthe
cells on them were packaged in 50mL sterile tubes samplingsiteof“Drty”in[22]).Endolithiccommunity
and preserved in liquid nitrogen when transporting to sampleswerecollectedwithasterileironspoon(~0.5
the lab, and stored at −80°C until DNA extraction cm from the outside) in December, 2014, and stored
within24h.CommunitygenomicDNAwasextracted inaniceboxwhentransportingtothelabwithin24h,
fromthefiltersasdescribedpreviously[14].Thetotal and stored at 4°C until DNA extraction within 24h.
genomic DNA was sequenced using Illumina paired- DNA extraction methods were the same as used for
end (PE) 125 and 250bp kits on a HiSeq2500 AMD outflow filters, except that no liquid nitrogen
platform with an insert length of ~400bp. It should was used for breaking the cells. The total genomic
bepointedoutthat,thissamplingapproach(alsoused DNA was sequenced using Illumina paired-end (PE)
for FK_AMD_2010) may miss some smaller mem- 125and250bpkitsonaHiSeq2500platformwithan
bersofthepopulation,forsomemicrobialcommunity insert length of ~400bp.
members could pass through 0.22μm pore size filters (5) Los Azufres National Park (Me_Mat). The Los
[8, 21]. Azufres National Park is located in the west of
(2) Fankou AMD sediment (FK_Sedi). A core was Mexico, 65km east of the city of Morelia (N19.47°,
sampled from the Fankou AMD pond using a W100.39°). Microbial layer samples were collected
sediment-sampling collector in March, 2015. The duringApril2013attheLosAzufresgeothermalfield
sediment core was divided into five layers based on in a rocky area full of cracks from which hot
their different colors, which generally reflect their fumaroles emerge, just in front of the muddy lagoon
oxidativestatus.Eachofthefivelayerswerecollected known as Los Azufres Natural Spa. DNA was
in 50mL sterile tubes, stored in an icebox and extracted using the MOBIO Mega Prep Soil DNA
transported to the lab within 4h after sampling, and extraction kit, and sequenced on an Illumina Miseq
storedat4°CuntilDNAextractionwithin24h.DNA platform using a PE 250bp kit.
of the AMD sediment was extracted using the
FastDNA SPIN kit (MP Biomedicals, Irvine, CA,
USA) according the manufacturer’s instructions. The Metagenome assembly and genome binning
extracted DNA was individually amplified using
multiple displacement amplification (MDA) with the All the raw data of the metagenomic data sets were filtered
GenomiPhi DNA amplification kit (v3) (Amersham to remove low quality bases/reads using Sickle software
Biosciences, Piscataway, NJ), and purified for library (version1.33;[23]),withtheparameterssetto“-q=15and
construction and sequenced on an Illumina −l=50”(SupplementaryTable3).Afterwards,oneormore
HiSeq2500 platform with PE 150bp kits. quality datasets of each sampling site were assembled (or
(3) Los Rueldos (LR_AMD). This sampling site is co-assembled for the five AMD sediments, and for the two
located in the Morgao stream, 2km northeast of the Los Rueldos AMD streamers) using SPADES (version
town of Mieres and 20km southeast of Oviedo, the 3.7.1; [24]) with a proper kmer set (Supplementary
capital city of Asturias, NW Spain (N43.27°, Table 3), under the “--meta” model. To calculate scaffold
W5.77°). Previous analyses have revealed the exis- coverage, all high-quality reads were mapped to the
tence of ARMAN taxa in the AMD of this site [16]. assembled scaffolds using BBMap, with the parameters set
Twomacroscopicstreamersamples(oneoxicandone as “minid=0.97, local=t”. The coverage of each scaffold
anoxic) thriving along an AMD outflow were was determined based on the mapping results using Meta-
collected, and filtered with 0.45μm filters. The BAT [25]. For scaffolds assembled from co-assembled
genomicDNAwasextractedfromthefiltrationwithin quality datasets, each metagenomic dataset was separately
2h as previously described [16], and amplified using mapped to the scaffolds, and the coverage of each dataset
the GenomiPhi DNA amplification kit (v3) (Amer- was calculated. The scaffolds were binned using MetaBAT
sham Biosciences, Piscataway, NJ). The DNA with default parameters, considering both tetranucleotide
samples were used for library construction, and frequencies and coverage of the scaffolds (multiple
MetabolicversatilityofARMAN 759
coverages for co-assembled AMD sediments and Los set as “minid=0.97, local=t”. The mapped quality
Rueldos AMD streamers). reads of each genome bin were assembled using
SPADES (version 3.7.1) with the k-mer set the same
Identification of ARMAN associated genomic bins as during the metagenomic assembly, under the
“--careful” model.
The bins obtained from MetaBAT were evaluated for (4) removal of contaminations from the reassembled
taxonomic assignment, genome completeness and potential genome bins. This step was performed as stated
contamination, using the “lineage_wf” pipeline in CheckM above in step (2). The optimized genome bins were
[26].TodeterminewhichofthebinsbelongtotheARMAN visualizedusingESOM(SupplementaryFig.3).After
phylum, all bins that were assigned to the lineage of genome optimization, a set of 54 conserved archaeal
“Archaea” by CheckM were further assessed using the fol- single-copy genes (SCGs) [3] were retrieved to
lowing procedures: estimate genome completeness (Supplementary
Table4).Additionally,thegenomebinswereassessed
(1) comparing their predicted genes (amino acid
using CheckM (Packs et al. 2015) to estimate the
sequences, which were received from the individual
degree of contamination (Supplementary Table 5).
scaffoldsusingProdigalincludedinCheckM)against
the NCBI-nr database using the DIAMOND software
[27] under the BLASTp model. Genome bins analyses
(2) comparing the scaffolds against the SILVA 16S/18S
dataset using BLASTn. All genomic bins with >50% The optimized genome bins were submitted to the IMG/
of genes having a besthit to ARMAN spp. from Iron MER system for gene calling and annotation (Supplemen-
Mountain (ARMAN-1, −2, −4, and −5; see above), tary Table 5). Protein-coding genes were also assigned to
and/or with a 16S rRNA gene sequence closely functionalorthologsoftheKEGGOrthologydatabaseusing
related to ARMAN spp. (sequence identify >75%), KAAS [29] (parameters: GHOSTX, Genes dataset, SBH
were considered as ARMAN genome bins for further assignment method). Carbohydrate-active enzymes were
evaluation. identified using the carbohydrate-active enZYmes (CAZy)
databaseonthedbCANwebserver [30].Asubset ofCAZy
genes, which were shown to be involved in specific car-
Genome bins optimization
bohydrate degradation pathways [31], was selected for
further analysis. Peptidases wereidentifiedusingMEROPS
TheARMANgenomebinswerefurtheroptimizedtoobtain
[32]. All the genome bins were also uploaded to RAST
high-quality genomes as follows:
(Rapid Annotation using Subsystem Technology) [33, 34]
(1) combination of genome bins. The automatic binning for gene prediction and annotation, using the Network-
tool MetaBAT may separate a “true” genome bin into Based SEED API using the command svr_submi-
two ormore smaller, separate bins. Tocombinethese t_RAST_job, selecting the domain Archaea and the RAST
smaller ARMAN genome bins all scaffolds were gene caller method. Metabolic pathways were constructed
visualized using ESOM [28] and bins that clustered for Micrarchaeota and Parvarchaeota based on gene anno-
together and shared a similar coverage range were tation results. The structural prediction of proteins was
combined and evaluated by CheckM. performed by the online tool PRED-SIGNAL with default
(2) removal of potential contaminations. Protein-coding parameters [35].
genes from each ARMAN genome bin (also those
combined ones) were determined using Prodigal KEGG Orthology based cluster analyses of ARMAN
included in CheckM (see above) and the acid amino genomes
sequences were compared against the NCBI-nr
database using DIAMOND (“blastp” model; [27]). To have a global comparison of metabolisms of the
Each scaffold was evaluated for contamination based ARMAN genomes, the occurrence of KEGG Orthology
onthetaxonomicalassignmentofitsgenes;andthose (KO) in each Micrarchaeota and Parvarchaeota genome
scaffolds with abnormal coverage and/or GC content were summarized based on KAAS annotation results (see
were deleted or moved to another ARMAN bin. above) (Supplementary Table 6), then the data was trans-
(3) re-assembly of genome bins. Quality metagenomic ferredto(0,1)data(foragivenKOinagivengenome,ifthe
reads were mapped to the scaffolds of each genome count number is 1 or more, then set as 1, if the count
bin obtained from the former two processing steps number is 0, then set as 0). Hierarchical clustering analysis
described above, using BBMap with the parameters of all the Micrarchaeota and Parvarchaeota genomes was
760 L-XChenetal.
respectively performed based on the (0,1) data, with cor- tree topology was uploaded to iTOL v3 [43] for visualiza-
relation distance and complete-linkage method using the R tion and formatting.
package of “pheatmap” (version, 0.7.7).
Relative abundance and community composition
16S rRNA gene sequences retrieval analyses
ARMAN-related 16S rRNA gene sequences were retrieved The single-copy ribosomal protein S3 (rpS3) was used as a
from NCBI GenBank (Supplementary Fig. 1). The 16S phylogeneticmarkerfor taxonomicassignmentandrelative
rRNA gene sequences of ARMAN-1 and ARMAN-2 abundance analyses as previously reported [3]. For each
(Micrarchaeota) and ARMAN-4 (Parvarchaeota), were metagenomic assembly result, all rpS3 sequences were
used as query to search for related 16S rRNA gene identified by CheckM [26], and the coverage was deter-
sequences in GenBank using BLASTn. Hits with an mined based on the coverage from the corresponding
alignment coverage ≥50% and sequence similarity ≥75% scaffold (see above). The relative abundance of each taxa
(phylum level threshold as determined using 16S rRNA was calculated by dividing the coverage of its rpS3 by the
genesequences;[36])weredownloadedfrom theGenBank total coverage of all rpS3 proteins in the community.
datasets. Additionally, 21 Micrarchaeota 16S rRNA gene
sequences were obtained from the IMG/M system, which Micrarchaeota and Parvarchaeota enrichment
were also included in this study [37]. For phylogeny ana- experiment
lyses of these retrieved 16S rRNA gene sequences with
those from the obtained genome bins (see below for An enrichment attempt of ARMAN cells was performed
methods), an initial tree was constructed to remove those under microaerobic/anaerobic conditions in laboratory, to
sequences not related to ARMAN, as determined by the illustrate they could grow with low O [11, 17]. We have a
2
phylogenetic topology and sequence similarity between stimulated AMD systemlocated inthecampusofSun Yat-
eachother(≥75%).Finally,weobtained176Micrarchaeota sen University, which is a 1×1×1m cement-lined pit
and 114 Parvarchaeota related 16S rRNA gene sequences with approximately 5 t natural pyrite particles (from Yunfu
with a length ≥800bp (Supplementary Table 7). sulfur/iron mine of Guangdong, China), and purified water
was injected into the pit monthly to accelerate the pyrite
Phylogenetic analyses acidification processes, and the microbial communities
dynamics have been reported previously [44].
Atotalof16ribosomal proteinsequences(i.e.,L2,L3, L4, In August of 2016, the AMD outflow samples collected
L5, L6, L14, L15, L16, L18, L22, L24, S3, S8, S10, S17, from thestimulated systemwere filtered throughpore sizes
and S19) from the ARMAN genomes and NCBI Archaea of 0.45μm filters, and the filtrate was collected into sterile
genomes (including selected genomes from all known tubes and stored at 4°C for enrichment within 24h (the
archaeal phyla) were individually aligned using Muscle filtrate will be referred as “inoculum” hereafter). Subse-
[38], and filtered by TrimAL [39] to remove columns quently,50mlinoculumwascentrifuged(4°C,12,000×g,
comprised of more than 95% gaps, then the 16 filtered 30min) and the supernatant was discarded, and genomic
ribosomal protein sequences were concatenated for each DNA was extracted using the FastDNA SPIN kit (MP
genome.Alltheconcatenatedsequencesrepresentedatleast Biomedicals, Irvine, CA, USA) according the manu-
50% of the expected alignment columns (29 out of 33 facturer’s instructions, and PCR analyses using ARMAN-
genomes spanned >80% of the expected alignment col- specificprimerset(ARM979F/ARM1365R;[11])indicated
umns). The concatenated sequences were used to recon- the occurrence of ARMAN cells in the inoculum (data not
structaphylogenetictreeusingRAxMLversion8.1.24[40] shown).
withtheparameterssetas“-fa-mPROTGAMMALG”.16S The culture medium of the novel and extremely acid-
rRNA gene sequences from all ARMAN genomes (see ophilic, cell-wall-deficient archaeon Cuniculiplasma divul-
above; deleted insertion with a length ≥10bp) and other gatum (primarily known as G-plasma of
archaea [2] were aligned using the SINA alignment algo- Thermoplasmatales),whichwasreportedtoco-culturewith
rithm [41] through the SILVA web interface [42], the filterable archaea with reduced genome size [45], was used
alignment was filtered by TrimAL to remove columns forARMANenrichment.Theculturemediumwasmodified
comprised of more than 95% gaps and a phylogenetic tree by adjusting the pH value to 3.0 with H SO , and addition
2 4
wasconstructedusingRAxMLversion8.1.24[40]withthe of1gNa S O each100ml,thensterilizedat115°Cfor25
2 2 3
parameters set as “-f a -m GTRGAMMA -n boot -c 25 -p min.Toeach100mlofmodifiedandsterilizedmedium,1g
12345 -x 12345”. The resulting newick file with the best of sterilized nutriment was added, that is, (1) yeast extract,
(2) beef extract and (3) peptone, additionally, a control
MetabolicversatilityofARMAN 761
Table1 Generalfeaturesofthe
Newlyreconstructed Previouslypublisheda
ARMANgenomesthatnewly
reconstructedandpreviously Micrarchaeota Parvarchaeota ARMAN- ARMAN- ARMAN- ARMAN-
publishedARMANgenomes 1 2b 4 5
No.ofgenomes 27 12 1 1 1 1
No.ofscaffolds 3–59(22) 7–78(32) 66 1 44 73
Totallength(Mb) 0.66–1.08 0.64–0.96 0.87 1.01 0.79 0.92
(0.85) (0.80)
GCcontent(%) 28.3–54.2 32.9–44.5 45.6 48 34.4 34.2
(44.5) (37.4)
No.oftRNAgenes 27–48(37) 30–48(42) 36 38 43 45
No.ofprotein-coding 743–1181 701–1095 993 1069 933 1076
genes (912) (897)
Genelength(bp) 793–891(843) 767–857(820) 793 859 767 770
Overlappedgenes(%) 7–15(10) 17–20(19) 9 10 18 20
Codingdensity(%) 87.3–92.7 90.4–93.4 89.4 91.3 90.7 90.8
(90.4) (92.2)
Estimatedcompletenessc 78–100(94) 78–98(92) 89 100 80 85
Contamination(%)d 0–2.8(0.7) 0–2.9(1.2) 0 0.9 3.3 5.6
aThe previously published ARMAN genomes were analyzed (annotation and others) as for newly
reconstructedones
bTheclosedgenomeofARMAN-2isreportedinthisstudy
cFor54archaealsingle-copygenes(SCGs)previouslyusedforgenomecompletenessevaluationasin[3].
FourMicrarchaeotagenomesweredetectedwithallthe54SCGs(SupplementaryTable4),whilenotclosed
dThecontaminationofeachgenomewereevaluatedwithCheckM[26]usingasetof149archaealSCGs.
Thoseinbracketsareaveragenumbers.SeeSupplementaryTable5fordetails
withoutnutrimentadditionwasincludedinthisexperiment. ng/μl),1μlwasdiluted×10times,and1.2μlofthediluted
Afterwards, 1mL of inoculum (see above) was added into DNA(~10ng)wasamplifiedusingMDAintriplicates.The
100mL of each of the four mediums (in 250mL culture three amplified DNA samples and also the raw inoculum
bottle).Theculturebottlesweresealedandlocatedinadark DNAsample(non-MDA)weresequencedformetagenomic
and non-shaking incubator for growth (37°C). analyses(qualitycontrol,co-assembly,coveragecalculation
After 14 days, cells from all four mediums (~100ml) as described above). The relative abundance of taxa in the
were individually collected by centrifugation (4°C, non-MDA sample and three MDA samples were evaluated
10,000×g,15min),genomicDNAwasextractedfromthe based on the rpS3 sequences as described above.
precipitates [14]. PCR using ARMAN-specific primer set
was performed to detect the ARMAN occurrence, the Phylogenetic analyses of interested genes
resultsshowedallexceptingthenegativecontrol,contained
ARMANcells(SupplementaryFig.4).Duetotherelatively Phylogenetic analyses were performed individually for
lowbiomassofthesecommunities,thethreeextractedDNA specificgenesinvolvedintheTCAcycleandtherespiratory
samples with ARMAN cells were amplified by MDA chain.Forexample,foragivengene,theproteinsequences
(GenomiPhi DNA amplification kit (v3); Amersham Bios- in all available Micrarchaeota genomes (33 in total) and
ciences, Piscataway, NJ), and sequenced for metagenomic Parvarchaeota genomes (16 in total) were identified and
analyses. De novo assembly (co-assembly) and genome extracted, then compared against the NCBI-nr database for
binning was performed, and the community composition theirclosedhomologiesusingBLASP(e-valuethreshold=
and relative abundance of taxa were analyzed based on 1e−5). The 20 top hits were extracted and double checked
ribosomal protein S3 sequences as described above. for their functional annotation, those ones with correct
PreviousstudiesshowedthatMDAtreatmentcouldskew function were retained, and aligned with those from the
therelativeabundanceoftaxaandcompromisequantitative abovementioned genomes using Muscle with default para-
analysisofmetagenomes[46],wethuscomparedtheresults metres [38], and filtered by TrimAL [39] to remove col-
ofnon-MDAandMDAmetagenomesoftheinoculum(see umns comprised of more than 95% gaps, then a
above).FortheinoculumDNAsample(concentration=89 phylogenetic tree was built using RAxML version 8.1.24
762 L-XChenetal.
[40]withtheparameterssetas“-fa-mPROTGAMMALG” gene length (820bp vs. 843bp), higher coding density
(corresponding gene sequences from E.coli strains were (92.2% vs. 90.4%) and higher frequency of overlapping
used for outgroup). The resulting newick file with the best genes(19%vs.10%)(Table1andSupplementaryTable5).
tree topology was uploaded to iTOL v3 [43] for visualiza-
tion and formatting, all the trees were rooted between the Relative abundance and microbial composition
outgroup and others.
Though25outof39newgenomeswerereconstructedfrom
Accession codes metagenomic data sets obtained by filtration for small cells
(Materials and methods), 17 out of these 25 genomes
All newly constructed Micrarchaeota and Parvarchaeota showed a relative abundance of <1% (Supplementary
genomes reported in this study have been deposited in the Fig. 5). However, one Micrarchaeota genome (FK_Sedi_-
IMG/MER system under the accession numbers listed in bin_12_4) was detected with a relative abundance of
Supplementary Table 5. 5.3–21.3% in AMD sediment layers of FK_Sedi_C1_3,
C1_4 and C1_5), indicating its potential significance in the
corresponding AMD sediments.
Results and discussion Interestingly, while the microbial composition of the
analyzedcommunitieswereremarkablydifferentfromeach
New Micrarchaeota and Parvarchaeota genomes other (Supplementary Fig. 6), the “alphabet plasmas” (e.g.,
reconstructed from metagenomes A-, E-, G-, I-plasma; [47]) belonging to the order Ther-
moplasmatales (phylum, Euryarchaeota; class, Thermo-
A total of 10 unpublished metagenomic datasets were plasmata) were detected in all of them. “Alphabet plasma”
obtainedfromtwoAMD(oneAMDoutflowandfiveAMD are abundant inhabitants in AMD ecosystems [47, 48], our
sediment metagenomes from Fankou mine tailings of analyses indicated they may also favor acidic hot spring
China, and two AMD streamer metagenomes from Los related environments.
Rueldos of Spain) and two hot spring related environments
(one hot spring endolithic metagenome from Tengchong of Phylogeny and biodiversity of Micrarchaeota and
China, one hot spring mat metagenome from Los Azufres Parvarchaeota
National Park of Mexico), and another two data sets were
retrieved from public databases (one AMD outflow meta- Based on a phylogenetic analysis of 16 concatenated ribo-
genome from Fankou mine tailings of China [15]; and one somal protein sequences, both Micrarchaeota and Par-
hot spring sediment metagenome from Yellowstone varchaeota branched in the DPANN superphylum of the
National Park of USA) (Supplementary Figs. 1 and 2). archaeal domain (Fig. 1), Micrarchaeota clustered with
Similar to the Iron mountain where ARMAN were first Diapherotrites, and Parvarchaeota clustered with
reported, all these sampling sites except the Yellowstone Nanoarchaeota, Pacearchaeota, and Woesearchaeota. The
National Park Obsidian pool were characterized by low pH 16S rRNA gene sequences based phylogeny showed a
values (Supplementary Table 2). different pattern as reported in a previous study [3], with
A total of ~200Gbp of raw sequences were obtained Parvarchaeota clustered with Pacearchaeota and Woe-
from all these 12 metagenomes (Materials and methods; searchaeota, and then clustered with Micrarchaeota (Sup-
SupplementaryTable3).Qualitycontrol,denovoassembly plementary Fig. 7). With the availability of sufficient
and subsequent genome binning from these metagenomes genomes obtained in this study, both phylogeny revealed
resultedinatotalof27Micrarchaeotaand12Parvarchaeota thatMicrarchaeotaandParvarchaeotaarenotmonophyletic
genomes (Supplementary Fig. 3). Based on the occurrence but two distinct phyla, and should not be combined as one
of54archaealsingle-copygenes[3],thenewlyconstructed [4].
genomes showed a completeness of 78–100% (93% on Based on similarity thresholds of 16S rRNA gene
average), and a very low degree of genomic contamination sequences for taxonomic level definition (genus, 94.5%;
estimates (0.8% on average) (Supplementary Table 4). family, 86.5%; [36]), the previously published ARMAN-1
These genomes have a genome size approximating and -2 represent two genera in a single family, as well as
0.64–1.08Mb, and on average encoding 912 genes for ARMAN-4and-5(SupplementaryFig.7).Inthisstudy,the
Micrarchaeota, and 897 genes for Parvarchaeota, which newly reconstructed Micrarchaeota genomes represented at
were comparable to those of the four published ARMAN least 12 genera within two families, and the Parvarchaeota
genomes [11, 12] (Table 1 and Supplementary Table 5). genomesrepresentedatleastthreegenerawithinonefamily,
ComparedwiththeMicrarchaeotagenomes(unpairedttest, greatly expanding the phylogenetic and genomic diversity
P<0.05), Parvarchaeota have significantly shorter average of these two lesser known phyla.
MetabolicversatilityofARMAN 763
Fig.1 Phylogeneticanalysesof 50 Euryarchaeota
Micrarchaeotaand
ParvarchaeotagenomesThe 3 Heimdallarchaeota
phylogenetictreewas 1Odinarchaeota
Asgard Hot spring related
constructedbasedon16 2 Lokiarchaeota
YNP hot spring (public)
concatenatedribosomalprotein 4 Thorarchaeota TC endolithic community (this study)
sequencesfromeach 1Korarchaeota Mexico microbial mat (this study)
Micrarchaeotaand
13 Thaumarchaeota AMD related
Parvarchaeotagenomeand
FKAMD outflow 2014 (this study)
referencearchaealgenomes(or 4 Aigarchaeota FKAMD outflow 2010 (public)
scaffoldswiththetarget 2 Geothermarchaeota FKAMD sediment (this study)
ribosomalproteinsequences) TACK IMAMD biofilm (public)
(Materialsandmethods).The 13 Bathyarchaeota LRAMD streamer (this study)
threearchaealsuperphyla 33 Crenarchaeota
TACK,AsgardandDPANN, 5 YNPFFA
andthenumberofgenomes To Bacteria
3 Verstraetearchaeota
includedforeachphylumare
shown.Previouslypublished 5 Altiarchaeales (Euryarchaeota)
genomesareindicatedby
3 Diapherotrites
arrows,andthemetagenomic
datasetsof“YNPhotspring”and YNP_OFPK__bAinM_D2_421010_bin_24
“FKAMDoutflow(2010)”were FK_AMD_2014_bin_90
FK_AMD_2014_bin_107
retrievedfrompublicdatabases.
FK_AMD_2014_bin_113
Bootstrapvaluesarebasedon FK_Sedi_bin_35
100replicates FK_AMD_2014_bin_115
FK_Sedi_bin_28
TC_Endo_bin_11
FK_AMD_2014_bin_52
FK_AMD_2014_bin_40
FK_Sedi_bin_12_4 Micrarchaeota
FK_AMD_2014_bin_97
LR_AMD_bin_1
FK_AMD_2014_bin_20
FK_AMD_2010_bin_7
FK_AMD_2014_bin_88
FK_AMD_2014_bin_31
TC_Endo_bin_32
ARMAN-2
FK_AMD_2010_bin_22 DPANN
FK_AMD_2014_bin_121
FK_AMD_2014_bin_73
FK_AMD_2014_bin_50
TC_Endo_bin_6
FK_AMD_2014_bin_108
FK_AMD_2014_bin_36
Me_Mat_bin_1
ARMAN-1
9 Nanohaloarchaeota
5 Aenigmarchaeota
9 Woesearchaeota
10 Pacearchaeota
4 Nanoarchaeota
FK_AMD_2014_bin_57
FK_AMD_2014_bin_27
TC_Endo_bin_24
LR_AMD_bin_3
ARMAN-5
TC_Endo_bin_23
ARMAN-4 Parvarchaeota
TC_Endo_bin_47
Tree scale: 0.5 FK_AMD_2014_bin_17
FK_AMD_2010_bin_6
Bootstrap ≥ 70% and < 90% FK_AMD_2014_bin_2
FK_AMD_2014_bin_10
Bootstrap ≥ 90% FK_AMD_2010_bin_5
FK_Sedi_bin_55_70
Environmental distribution of Micrarchaeota and environments, including hypersaline mats, soils, peat,
Parvarchaeota freshwater lakes, and underground water, indicating their
potential highly adaptive capacities and phylogenetic
Comparison of Micrarchaeota 16S rRNA gene sequences diversity of at least two classes (Supplementary Fig. 8). In
from around the world revealed their considerable biodi- contrast for Parvarchaeota, all the retrieved 16S sequences
versity, especially within AMD and hot springs. Micrarch- from NCBI GenBank were from AMD-related environ-
aeota were also found in a variety of non-acidic ments, and were affiliated within the same family as found
764 L-XChenetal.
withintheTengchongacidicendolithiccommunityreported showed the opposite trend (Supplementary Fig. 5 and
in this study (Materials and methods). The limited dis- Supplementary Table 2).
tribution of the known Parvarchaeota may be due to their
lowerphylogeneticdiversity,indicatingmoreParvarchaeota Comparative genomic analysis based on KEGG
clades remain to be discovered, while their low abundance Orthology
in nature may hinder this process.
For those Micrarchaeota and Parvarchaeota spp. with To reveal if there are any clade specific metabolisms for
genomic information obtained in this study, some of them bothphyla,ahierarchicalclusteringanalysiswasperformed
showed clear environmental preference. Micrarchaeota basedontheoccurrenceofKEGGOrthology(KO)patterns
genus 12 and Parvarchaeota genus 1 and 2 members in each genome (Fig. 2) (Materials and methods). Com-
(Supplementary Fig. 7) were all obtained from environ- paredtothephylogenyanalysesresultsshowninFig.1(left
ments with higher temperature (Supplementary Table 2). panels of Figs. 2a, b), different clades in both phyla gen-
And the three Micrarchaeota members from Fankou AMD erallyhavetheiruniquegenecontents,indicatingvariations
sediment tended to have higher relative abundance in the in metabolic potentials between clades. For Micrarchaeota,
lower layers, which were characterized by a higher Fe2+ the non-oxidative pentose phosphate pathway, aspartate-
concentration and lower Fe2+/Fe3+ ratio, while the only semialdehyde dehydrogenase, saccharopine dehydrogenase
Parvarchaeota member from these layers (FK_Sedi_55_70) and zinc transporter genes were found only in group 1,
a
Phylogeny cluster KOs cluster
YNP_OP_bin_241 FK_AMD_2014_bin_107
FK_AMD_2010_bin_24 FK_AMD_2010_bin_24
FK_AMD_2014_bin_90 FK_AMD_2014_bin_113
FK_AMD_2014_bin_107 Group 1 FK_AMD_2014_bin_90
FK_AMD_2014_bin_113 FK_Sedi_bin_28
FK_Sedi_bin_35 YNP_OP_bin_241
FK_AMD_2014_bin_115 FK_AMD_2014_bin_88
FK_Sedi_bin_28 FK_AMD_2014_bin_115
TC_Endo_bin_11 FK_Sedi_bin_12_4
a FK_AMD_2014_bin_52 Group 2 FK_AMD_2014_bin_97
ot FK_AMD_2014_bin_40 FK_AMD_2014_bin_73
ae FK_Sedi_bin_12_4 FK_AMD_2014_bin_121
ch FK_AMD_2014_bin_97 Group 3 FK_Sedi_bin_35
Micrar FLFLRRKK____AAAAMMMMDDDD____bb22ii00nn11__4411__bbiinn__2200 Group 4 FK_AMD_201A0R_MbiAn_N2-22
FFKK__AAMMDD__22001100__bbiinn__77 TC_Endo_bin_11
FK_AMD_2014_bin_88 FK_AMD_2014_bin_40
FK_AMD_2014_bin_31 Group 5 FK_AMD_2014_bin_52
TC_Endo_bin_32 TC_Endo_bin_6
ARMAN-2 FK_AMD_2014_bin_108
FK_AMD_2010_bin_22 Me_Mat_bin_1
FK_AMD_2014_bin_121 Group 6 ARMAN-1
FK_AMD_2014_bin_73 FK_AMD_2014_bin_36
FK_AMD_2014_bin_50 FK_AMD_2014_bin_50
TC_Endo_bin_6 FK_AMD_2014_bin_20
FK_AMD_2014_bin_108 FK_AMD_2010_bin_7
FK_AMD_2014_bin_36 Group 7 LR_AMD_bin_1
Me_Mat_bin_1 TC_Endo_bin_32
ARMAN-1 FK_AMD_2014_bin_31
b
Phylogeny cluster KOs cluster
FK_AMD_2014_bin_57 FK_AMD_2014_bin_57
FK_AMD_2014_bin_27 FK_AMD_2014_bin_27
TC_Endo_bin_24 Group 1 TC_Endo_bin_24
LR_AMD_bin_3 TC_Endo_bin_23
a ARMAN-5 ARMAN-4
ot TC_Endo_bin_23 ARMAN-5
e
a ARMAN-4 FK_AMD_2014_bin_2
h
varc FTKC__AEMndDo__2b0in1_44_7bin_17 Group 2 FK_FAKM_SDe_d2i_0b1i4n__b5i5n__7100
ar FK_AMD_2010_bin_16 FK_AMD_2010_bin_5
P
FK_AMD_2014_bin_2 LR_AMD_bin_3
FK_AMD_2014_bin_10 TC_Endo_bin_47
FK_AMD_2010_bin_5 Group 3 FK_AMD_2014_bin_17
FK_Sedi_bin_55_70 FK_AMD_2010_bin_6
Fig. 2 Comparative analyses of gene contents of Micrarchaeota and genome;seeMaterialsandmethods)werecompared,andthegenomes
Parvarchaeotagenomes.ForaMicrarchaeotaandbParvarchaeota,the weremanuallyassignedtoseveralclades/groupsbasedonthecluster
phylogenyclusterpattern(fromFig.1)andthecorrespondingKEGG patterns.Thesameclade/groupinaphylogenyclusterandgenecon-
Orthology clustering pattern (based on occurrence of KOs in each tentsclusterswerelinkedwithasolidline
MetabolicversatilityofARMAN 765
Carbon monoxidedehydrogenase CooS andarginasegenes antibiotics, heat shock, heavy metals, and oxidative stress
occurredonlyingroup3,andammoniumtransportergenes (Fig. 3 and Supplementary Table 8). Interestingly, 15
wereonlydetectedingroup 7(Fig. 2a).ForParvarchaeota, Micrarchaeota and 10 Parvarchaeota genomes were detec-
glutamine synthetase and L-asparaginase genes were ted with genes encoding resistance protein for fosmido-
exclusively found in group 1 and heptose III glucur- mycin,whichisaninhibitorofisoprenoidbiosynthesis[52].
onosyltransferase and inositol transporter genes in group 3 Almost all genomes of both phyla carry a superoxide dis-
(Fig.2b).Itisreasonabletospeculatethatthespecificgene mutaseandanalkylhydroperoxidereductaseallowingfora
contents of different clades may help them to inhabit dis- quick response to oxidative stress. Notably, Micrarchaeota
tinct niches when they coexist with each other. A similar carry a Fe and the Parvarchaeota carry a Mn superoxide
trendhasbeenobservedformembersofThermoplasmatales dismutase, likely indicating their distinct evolutionary his-
AMD archaea, the “alphabet plasmas”, which were also tories of obtaining this function. Interestingly, evidence of
detectedinourmetagenomicsamples.Ithasbeenfoundthat MicrarchaeotaactingasapathogentargetingotherBacteria
theThermoplasmatalesAMDarchaeadifferentiatebysubtle was detected in three genomes that encode for a lytic
genomic differences that allow their co-existence even if murine transglycosylase [4], which could be transported
they share a great number of metabolic capabilities [49]. outside of the cell via the Sec-SRP secretion system
(Fig. 3a).
Metabolic potentials of Micrarchaeota and
Parvarchaeota Amino acid and nucleotide biosynthesis
Inordertoresolvethephysiologicalcapabilitiesencodedin Although Micrarchaeota are able to generate alanine and
these genomes (including the four published ARMAN glutamate respectively from pyruvate and aspartate
genomes), we used a variety of functional gene database (Fig. 3a), these organisms lack biosynthetic pathways for
comparisons for annotation (Materials and methods). other amino acids, and Parvarchaeota lack biosynthetic
Though the metabolisms of Micrarchaeota and Parvarch- pathways for all amino acids. However, both phyla encode
aeota have been previously predicted based on three gen- for a variety of extracellular (e.g., S16, S53), membrane
omes [11]; with a higher number of genomes representing (e.g., M48, M50), and cytoplasmic (e.g., M01, M24) pep-
higher phylogenetic diversity, we now have a better tidases and transporters (e.g., amino acid permease) poten-
understanding of their metabolic versatility. The metabolic tially allowing them to scavenge amino acids from the
potentials of Micrarchaeota and Parvarchaeota based on all environment and thus to contribute to nitrogen cycling
43 genomes are detailed in the following several sections. (Fig. 3 and Supplementary Table 8). Both phyla appear to
lack genes for de novo biosynthesis of nucleotides from
Cell membrane biosynthesis phosphoribosyl pyrophosphate (PRPP) (Supplementary
Fig. 10), albeit some Micrarchaeota might produce PRPP
Isoprenoidsareessentialinalllivingorganisms,andvitalforcell from ribose-5-phosphate (Fig. 3a). However, complete
wallandmembranebiosynthesis.Asdetectedinotherarchaea, pentose phosphate pathways (for ribose-5-phosphate gen-
Micrarchaeotapossessgenesinvolvedinthemevalonatepathway eration) were not detected in these phyla, though some
for biosynthesis of isoprenoid precursors (i.e., isopentenyl Micrarchaeota genomes harbored genes for the non-
diphosphate and dimethylallyl diphosphate; IPP and DMAPP) oxidative pentose phosphate pathway. To date, complete
[50].Additionally,thepresenceofisopentenylphosphatekinase pentosephosphatepathwaysinArchaeahavebeendetected
genes indicates that they also have the alternative pathway for in only two Nanohaloarchaea [53, 54] and two Woe-
isopentenyl diphosphate synthesis (Supplementary Fig. 9). searchaeota [3]. Consequently, both phyla may acquire
Notably,nogenesofthemevalonatepathwayorthemethylery- nucleotides from the environment and/or other community
thritol 4-phosphate pathway [51] for IPP and DMAPP bio- members, for survival and cell proliferation, though free
synthesisweredetectedinParvarchaeota.However,considering nucleotidesmaybeveryunstableinanacidichabitatwhere
theirhighgenomecompleteness(SupplementaryTable4),itis MicrarchaeotaandParvarchaeotaareoftenfound(pHvalue
possiblethatthesepathwaysarepresentandtherelatedgenestoo 0.5–4.0; Supplementary Table 2).
noveltobeidentifiedviasequencesimilaritycomparison.Alter-
natively,theymayobtainisoprenoidsfromtheenvironment. Carbon fixation
Stress reponse Noneofthesixknowncarbonfixationpathways/cycleswas
detected in the Micrarchaeota and Parvarchaeota genomes,
Many genes related to environmental resistance were likely indicating that both phyla are heterotrophic. How-
detected in both phyla, including those for drugs, ever, the ribulose 1,5-bisphosphate carboxylase/oxygenase