Shotgun metagenomic sequencing reveals a high proportion of undescribed taxa

After quality filtering, the four libraries consisted of 46.624.517 s.d. 21.730.638 read pairs, for a total of 372.996.136 reads (Table 2). The four samples differed markedly in taxonomic composition and community structure (Fig. 2, Additional file 3: Figure S1). Relative abundances of the taxa inferred by MetaPhlAn2 show that Bacteria have the highest abundance 78.4% (196 taxa), followed by viruses with 15.6% (39), Archaea with 4.8% (12), Eukarya with 0.8% (2), and viroids with 0.4% (1). A total of 108 out of 250 (43.2%) unclassified taxa were identified, mainly at the species (84/108, 77.8%) and genus (23/108, 21.3%) levels. In only one case an unclassified order was found, belonging to the class Opitutae. Consistently with 16S profiling, Proteobacteria dominated the water microbiome in all sampling sites, although in hydra and spring, the dominance was less pronounced (68.9% and 72.2% versus 84.4% and 84.5% in storage tank and bathtub, respectively).

 

Open in a separate window

Fig. 2

Taxonomic tree of taxa detected in Comano Terme water by shotgun metagenomic sequencing. The external rings represent microbiome composition at the four sampling sites, with maximum color intensity corresponding to a relative abundance > 1%. Branches are color-coded according to the domain of life: yellow, Archaea; red, Bacteria; green, Eukaryota; purple, viroids; blue, viruses. Shapes of the leaves of branches are coded as follows: small dots, described species; stars, unclassified species; circles, unclassified genera; diamonds, unclassified orders. Taxa are abbreviated as follows: C, Caulobacterales; S, Sphingomonadales; P, Pseudomonadales; Rs, Rhodospirillales; Rb, Rhodobacterales; X, Xanthomonadales; My, Myxococcales; Rc, Rhodocyclales; Mc, Methylocystaceae; Br, Bradyrhizobiaceae; and Mx, Moraxellaceae

The core microbiome included the orders Sphingomonadales, Rhizobiales, and Caulobacterales, and the families Bradyrhizobiaceae and Moraxellaceae. The majority of core taxa, as for instance Actinomycetales and Burkholderiales, were found at a different relative abundance in the four sampling points. Viruses represent the second most abundant source of DNA and are rather homogeneously distributed between all samples. The family Methylocystaceae (alpha-Proteobacteria) was among the major components of the hydra microbiome (39.5%), whereas in spring, storage tank, and bathtub its abundance was much lower (0.4%, 0.03%, and 0.1% respectively) (Additional file 3: Figure S1a).

Isolate genomes and MAGs cover a wide range of taxa and abundance classes

To select the isolate genomes to be sequenced, the culture collection was pre-screened by amplifying and sequencing the 16S rRNA gene from single colonies and by clustering the isolates sharing an identity above 99%. For each cluster, a representative isolate was chosen. Whole genomes of 77 out of 181 isolates were sequenced. After assembly and quality checks, five genomes were excluded. A summary of genome assembly statistics is presented in Additional file 4: Table S3 and Additional file 5: Table S4. The metagenomics assembly resulted in 29.850 s.d. 13.101 contigs, summing up to > 500 Mbp (Table 1). After concatenation and redundancy removal with CD-HIT, the number of contigs was 114.739 with a total length of 255.640.704 bp. The initial automatic binning process resulted in 63 bins (MAGs), but many of them were either contaminated (probably resulting from a mixture of multiple, highly similar genomes) or shorter than 0.5 Mbp. After manual curation and quality filtering, 30 bins (accounting for 58.141.697 bp) passed the thresholds set by MIMAG. Ten out of 30 MAGs had > 90% completion and < 5% contamination, while 12 out of 30 MAGs had > 18 tRNAs, and in 2 MAGs all 3 ribosomal RNA genes were found (Additional file 4: Table S3). According to the single copy gene classifier implemented in Anvi’o, 21 bins belonged to a bacterial phylum, 3 bins were classified as Archaea, and 6 bins were assigned to the CPR super-phylum [12].

The genomes of the bacterial isolates included 4 phyla and 6 classes. Thirty-nine (54.1%) genomes resulted to be unclassified at different levels: 19 at the genus level, 16 at the family level, and 4 at the order level. While MAGs covered 4 phyla and 6 classes, for 13 and 7 of them, it was possible to infer their taxonomy at the order and at the species level, respectively (Additional file 5: Table S4). Isolates and MAGs were widely distributed in the microbial tree of life (Fig. 3).

 

Open in a separate window

Fig. 3

Phylogenetic position of bacterial isolates and MAGs on the bacterial tree of life. Orange dots, isolate genomes; yellow dots, MAGs. Branches that do not contain any isolate genomes or MAGs were collapsed for clarity. Collapsed branches are reported with an abbreviation on the external ring as Arch, Archaea; Ac, Actinomycetaceae; Al, Alteromonadales; Bc, Bacteroidales; Cb, Coriobacteriales; Ep, Erysipelotrichia; Lb, Lactobacillales; Sp, Spirochaetes; Se, Selenomonadales; and Te, Tenericutes. The external ring is colored according to the phylum

While all isolates could be taxonomically identified by PhyloPhlAn2, reaching the number of 100 aligned universal marker threshold for taxonomic placement, 9 MAGs failed to pass this threshold, probably as a consequence of their small genome size. Fragment recruitment analysis, i.e., the percentage of metagenomic reads mapping on a given MAG, revealed that 24 MAGs were reconstructed from hydra sample, including all Archaea and CPR. These accounted for 58.82% of the reads, while the two MAGs reconstructed mainly from the spring sample accounted for 3.38% and the four MAGs reconstructed mainly from the bathtub sample accounted for 9.53% of the reads, respectively.

ANI values calculated between each genome isolated (or assembled) from the spring water and their closest relative in the set of publicly available genomes [19], as inferred by PhyloPhlAn2, showed that 10 genomes belonged to known bacterial species (> 95% ANI over an alignment length > 1Mbp) (Fig. 4), 36 genomes (including 1 MAG) had ANI values ranging 85–92.5% over an alignment length of > 1 Mbp, and could be assigned as undescribed species belonging to known genera. The remaining 55 genomes had low ANI values (< 85%) and/or alignment length (< 1 Mbp). Few isolate genomes and most MAGs belonged to the latter group, with an alignment length ranging from few hundreds to few thousands base pairs and ANI values comprised between 80% and 85% (Fig. 4). In four cases, a MAG clustered within the same clade of an isolate with a relatively narrow phylogenetic distance. For those pairs, ANI values suggest that they may belong to the same family. Specifically, Bin_15_1 showed an ANI value of 72.5% with the isolate L1A1; Bin_9_1 had ANI value of 71.3% with 3R27C2-B; Bin_43, 70.8% with the isolate L1B and Bin_16_1 has an ANI of 68% with isolate 2PA. No MAGs were found to have a species or genus level relationship (ANI value > 95% or > 85%, respectively) with any of the isolates (Fig. 3).

 

Fig. 4

Average nucleotide identity (ANI) values (x-axis) plotted against log₁₀-transformed alignment length (y-axis). Each point indicates the best ANI hit between a given isolate genome (dots) or MAG (triangles) and its closest relative according to the database implemented in PhyloPhlAn2. Dots are color-coded according to the phylum. The threshold used for species and genus delineation is 95% and 85%, respectively

Strain isolation and identification

Water samples of 100 ml were concentrated by filtration using 0.22 μm Pall Supor filters (Pall Corporation, Port Washington, NY, USA). Filters were then placed on standard Petri dishes poured with four different growth media. In order to be able to culture the highest number of bacterial species, three common non-selective media (nutrient agar, tryptic soy agar, plate count agar) and one medium for enumeration of bacteria from potable water (Reasoner’s 2A agar) [33–35] were used for bacterial isolation. Growth media were prepared using filtered water from Comano Terme spring in order to maintain the native chemical composition of this environment.

To maximize the number of cultivable isolates, a large-scale culturing approach was set up following the protocol of Connon and Giovannoni [36]. Based on the microbial load evaluation, 2 μl of Comano Terme water (1–10 cells) were inoculated into 2 ml of growth medium in 96-well plates with square walls and round bottom; edges of the plate were filled with water to prevent edge effect. Four different media were used for liquid cultures, two media for enumeration of bacteria in water (Reasoner’s 2A agar and King’s B [37]), one non selective medium (nutrient broth), and one low-nutrient medium (M9 minimal medium). Each medium was prepared in five versions: standard, two pH variants (pH 5.5 and pH 8), one variant with addition or subtraction of carbon source (glucose, soluble starch, or both), and one amino acid variant (with or without yeast extract). Plates were replicated and incubated at 27 °C and 37 °C in both aerobic and anaerobic conditions. Anaerobic cultures were grown in a Ruskinn Concept 400 Anaerobic Workstation (Ruskinn Technology Ltd., Bridgend, UK) and maintained during incubation using GENBAG, Anaerobic Atmosphere Generator (bioMérieux, Marcy-l’Étoile, France). Evaluation of growth was performed for each well measuring O.D. 600 values with a TECAN Infinite m200 PRO microplate reader (Tecan, Männedorf, Switzerland). Wells displaying growth were then tested for purity performing two rounds of streaking on R2A plates. After 7 days of incubation at 27 °C, glycerol stocks were prepared. The 16S rRNA gene of each isolate was amplified by PCR for taxonomic identification. DNA was extracted using ReliaPrep gDNA extraction kit (Promega, Milan, Italy) and amplified using standard primers 8F-1492R [38]. Amplicons were sequenced using Sanger sequencing service by Eurofins (Eurofins Scientific, Luxembourg). Sequences with length > 400 bp and an average Phred score of > 25 were used for taxonomic identification of isolates using RDP classifier [39]. Assigned taxonomies were accepted only at estimated confidence > 95%.

16S amplicon sequencing

Two samples of 4 l were collected from each sampling point (3 replicates) in February 2016 and in August 2017. Total DNA was extracted using RapidWater DNA extraction kit (MoBio, Carlsbad CA, USA) with minor modifications: at step 5 of the protocol, the PowerWater Beat Tube were heated at 65 °C for 10 min and mechanical cell lysis was extended to 10 min for all samples. All other steps were performed following the manufacturer’s instructions. DNA extraction of 2017 samples was performed on triplicates and extraction yield remained consistent for all samples (data not shown). The V4 hypervariable region of the 16S rRNA gene was amplified by PCR with the 5PRIME HotMasterMix (Quanta BIO, Beverly, MA, USA). The length of the amplicons was 253 base pairs. Negative controls were included during sampling and main wet-lab steps. PCR blanks, DNA extraction blanks and DNA extractions from unused filters were prepared. Amplicons concentration, size range, and purity were measured using Agilent high sensitivity (HS)DNA kit on the Bioanalyzer 2100 (Agilent Technologies Italia S.p.A, Milano, Italy). Based on the molarity estimated using Bioanalyzer, each PCR product was diluted before pooling. The final pool was purified using the Agencourt AMPure XP DNA purification kit, following manufacturer’s instructions. Amplicons were sequenced on an Illumina MiSeq platform with 2 × 300 bp paired-end protocol.

The fastq files were quality checked using FastQC [40], and initial sequence analysis was performed with QIIME1 [41]; after demultiplexing and primer removal, reads were joined using default parameters. Joined reads were then filtered by average Phred quality score > 25 and a minimum length of 250 nucleotides; reads longer than 255 nucleotides were trimmed. Operational taxonomic units (OTUs) were built at 97% sequence identity with uclust [42]. Taxonomic assignment was performed with RDP classifier [39], using the Greengenes [43] reference dataset (release gg 13 5). Diversity analyses were performed using QIIME on a rarefied dataset to prevent statistical artifacts due to different sequencing depths. Alpha diversity was measured with Shannon metrics. Beta diversity measures were measured on unrarefied dataset calculating weighted UniFrac metrics and displayed through 2D PCoA [44].

Shotgun sequencing

Samples of 4 l were filtered using 0.22 μm Pall Supor filters (Pall Corporation, Port Washington, NY, USA). DNA extraction from filter was carried out using the RapidWater DNA extraction kit (MoBio, Carlsbad, CA, USA). Sequencing libraries were prepared using Nextera XT DNA Library Preparation Kit (Illumina, San Diego, CA, USA) following the manufacturer instructions, and quality checked using Caliper LabChip GX. Sequencing was performed on an Illumina HiSeq 2500 platform (Illumina, San Diego, CA, USA) with 100 bp paired-ends and an insert size of 250 bp. Raw sequencing reads were quality checked using fastq-mcf [45], and reads derived from human contamination were removed after mapping on the human genome (version hg19) using bowtie2 [46]. Adaptors and residual synthetic constructs were removed first mapping the reads on the genome sequence of the phage Phy174, and then using trim galore [47].

Microbiome profiling based on a database of taxon-specific marker genes was performed using MetaPhlAn2 [48]. The output of MetaPhlAn2 analyses was plotted using GraPhlAn [49].

Metagenome-assembled genomes (MAGs)

Each of the four samples was assembled separately using megahit [50], with the preset–meta-sensitive and setting the minimum contig length to 1000 bp. The assembled files were then concatenated, and redundant contigs were removed using CD-HIT [51], with 100% coverage and 100% id. The reads derived from the sequencing libraries were mapped on the resulting co-assembly. Downstream analysis was performed using Anvi’o v2.4.0 [52] and included the following steps: (i) taxonomic and functional annotation of the contigs; (ii) profiling the mapping files of single libraries on the co-assembly, to assess the coverage of each contig from each library; (iii) binning contigs into putative MAGs, based on the tetramer composition and differential coverage using CONCOCT [53], as implemented in the Anvi’o pipeline. A further quality check of MAGs relied on the presence of a set of single-copy genes [54]; the proportion of single-copy genes present in each MAG is considered an indicator of genome completedness, while the presence of multiple copies of the same single-copy gene (redundancy) is an indicator of possible contamination. Bins of contigs were manually curated to minimize the contamination by removing contigs (or groups of contigs) that showed markedly different coverage profile and/or tetrameric signature. The annotation of both MAGs and isolate genomes was performed using Prokka version 1.12 [55]. MAGs were classified into high and medium quality, following the guidelines for minimum information about metagenome-assembled genome (MIMAG) [56]. In brief, high-quality genomes had > 90% completeness, N50 > 10 Kb, < 5% contamination, < 500 contigs, encoded all three ribosomal RNA genes and at least 18 tRNAs; medium-quality genomes had > 50% completeness and < 10% redundancy. In addition, a score was calculated as S = completeness percentage − 5 × (redundancy percentage), as in Parks et al. [57]. For subsequent analyses, only high- and medium-quality genomes with an S score > 50 were considered. A dedicated pipeline for the assignment of bins to the candidate phyla radiation (CPR) super-phylum was setup, and the corresponding score was calculated using the CPR-specific single copy genes.

Genome sequencing and assembly

16S rRNA sequences obtained from the isolated strains were clustered using CD-HIT with a threshold of 99%. One representative isolate was selected from each cluster for genome sequencing. DNA extraction and preparation of sequencing libraries was performed as described in the shotgun sequencing section.

The quality of sequencing data was checked with MultiQC [58] and successively processed with Trim Galore [36] for quality and adapter trimming, with a quality threshold set at Phred 25. The filtered reads were de novo assembled using SPAdes v. 3.11.1 [59], using default settings and specifying kmer sizes -k 23,43,63,83,103,113. Quality of the assembled reads was verified using Quast [60]. Taxonomic placement was performed using PhyloPhlAn2 [61]; using the standard phylophlan database with 400 universal marker genes, a 2.0. Supermatrix_aa.cf. configuration, the option “—diversity” set to “high” and the “—fast” option on. The genome of the closest relative to each isolate was used to calculate the average nucleotide identity (ANI), using the software pyani [62].

We would like to thank Federica Armanini and Francesco Beghini for technical support, and Veronica De Sanctis and Roberto Bertorelli (NGS Facility at the Centre for Integrative Biology and LaBSSAH, University of Trento) for NGS sequencing.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.