Next generation sequencing data of a defined microbial mock community

Generating sequence data of a defined community composed of organisms with complete reference genomes is indispensable for the benchmarking of new genome sequence analysis methods, including assembly and binning tools. Moreover the validation of new sequencing library protocols and platforms to assess critical components such as sequencing errors and biases relies on such datasets. We here report the next generation metagenomic sequence data of a defined mock community (Mock Bacteria ARchaea Community; MBARC-26), composed of 23 bacterial and 3 archaeal strains with finished genomes. These strains span 10 phyla and 14 classes, a range of GC contents, genome sizes, repeat content and encompass a diverse abundance profile. Short read Illumina and long-read PacBio SMRT sequences of this mock community are described. These data represent a valuable resource for the scientific community, enabling extensive benchmarking and comparative evaluation of bioinformatics tools without the need to simulate data. As such, these data can aid in improving our current sequence data analysis toolkit and spur interest in the development of new tools.

Generating sequence data of a defined community composed of organisms with complete reference genomes is indispensable for the benchmarking of new genome sequence analysis methods, including assembly and binning tools. Moreover the validation of new sequencing library protocols and platforms to assess critical components such as sequencing errors and biases relies on such datasets. We here report the next generation metagenomic sequence data of a defined mock community (Mock Bacteria ARchaea Community; MBARC-26), composed of 23 bacterial and 3 archaeal strains with finished genomes. These strains span 10 phyla and 14 classes, a range of GC contents, genome sizes, repeat content and encompass a diverse abundance profile. Short read Illumina and long-read PacBio SMRT sequences of this mock community are described. These data represent a valuable resource for the scientific community, enabling extensive benchmarking and comparative evaluation of bioinformatics tools without the need to simulate data. As such, these data can aid in improving our current sequence data analysis toolkit and spur interest in the development of new tools.

Background & Summary
By definition, benchmark studies aim to provide standards that can be used to evaluate the performance of a process. The field of nucleic acid sequencing and sequence data processing has witnessed immense developments towards optimizing the balance of sequencing cost, precision and overall applicability to real-world questions. This progress has routinely relied on experimental setups of defined nature to critically rate novel approaches. In recent years, mock communities have been assisting in a variety of laboratory and computational test experiments, which resulted in quantitative and qualitative evaluation of corresponding studied methods. For example, mock communities were generated for the comparison of DNA extraction methods 1-3 , for the development of a dual-index sequencing and curation pipeline for Illumina MiSeq generated amplicon sequence data [4][5][6][7][8] , and to evaluate the Ion Torrent sequencing platform for gene-targeted studies 9,10 . Similarly, Pabinger et al. [11][12][13] used a mock community to benchmark MEMOSys, a web-based platform for metabolic models. The jumpstart consortium human microbiome project (HMP) data generation working group established a standardized protocol for ensuring high throughput consistency of 16S rRNA gene amplification and sequencing protocols by implementing a synthetic mock community of 21 known organisms, before finalizing their HMP 16S 454 protocol [14][15][16] . The HMP DNA and sequence data resources have not only enabled comprehensive characterization of the human microbiota, e.g. [17][18][19] , but also the use and development of a variety of advanced analysis tools. For example, chimera screening tools UCHIME and Chimera Slayer 1,3 , the OTU construction pipeline UPARSE 4,6-8 , and fine-tuned workflows for amplicon gene studies 9 used HMP data generated from mock communities. In contrast to the HMP mock, the synthetic community described here, MBARC-26 (Mock Bacteria ARchaea Community), is composed of organisms isolated from heterogeneous soil and aquatic environments as well as derived from human, bovine and frog (Table 1). MBARC-26 consists of 23 bacterial and 3 archaeal strains, belonging to the phyla Acidobacteria, Actinobacteria, Bacteroidetes, Deinococcus-Thermus, Firmicutes, (Alpha-and Gamma-)Proteobacteria, Spirochaetes, Thermotogae, Verrucomicrobia and Euryarchaeota. Genome sizes span 1.8-6.5 Mbp, GC contents vary between 28.4-72.7%, and repeat content ranges from 0-18.3% (Fig. 1, Table 1). All genomes are available as finished genome sequences in GenBank (Table 1). MBARC-26 DNA was shotgun sequenced on Illumina HiSeq 2000 and PacBio RSII sequencing platforms ( Table 2). We provide detailed descriptions of organism characteristics (Table 1), sample processing, including DNA extraction and quantification, sequencing library creation, and sequencing procedures ( Table 2). Data statistics encompass sequencing throughput characteristics (Table 2), community structure according to read mapping to reference genomes and according to molarity (Fig. 2 (Fig. 3b), and GC content analysis ( Supplementary Fig. 3). Due to inherent sequencing technology differences 11,13 , these two datasets are characterized by platform-, run mode-, and chemistry-specific read length, data throughput, GC and amplification bias, and error rate. We point out that our quantitative results are directly correlated to the respective sample preparation and sequencing methods used, as these have been shown to critically affect community representation 14,20 .
To date, several studies already utilized MBARC-26 and took advantage of its purposefully selected characteristics. Availability of complete reference genomes and relative abundance spread of individual constituents enabled determining lower limits of various metagenome library preparation protocols 14 . MBARC-26 was also used to develop a new full-length 16S rRNA gene amplicon sequencing protocol called PhyloTags 17 and allowed for quantitative comparison of amplicon to shotgun sequence data and bias evaluation associated with GC content. Using the MBARC-26 Illumina metagenome dataset and corresponding single-cell sequence data Bremges et al. developed MeCorS, a metagenome-enabled singlecell read correction tool 21 . To further encourage the use of this mock community, we report the release of molarity and shotgun sequence datasets of MBARC-26.
Perpetual community efforts to develop improved DNA sequence analysis software with various applications for shotgun sequence data requires standardized and well-characterized data for benchmark experiments. MBARC-26 was validated according to the specific sample processing tools using a variety of commonly used quality control methods, is accompanied by data statistics, and meant to enable method development and evaluation while enabling reproducibility of research findings.

Methods
These methods are expanded from descriptions in our previous work 17 .
Cultivation and DNA extraction DNA from Escherichia coli, Salmonella bongori, Salmonella enterica, Clostridium perfringens, Clostridium thermocellum and Streptococcus pyogenes was purchased from the American Type-Culture Collection (ATCC, Manassas, VA, USA). DNA from Fervidobacterium pennivorans, Thermobacillus composti and Corynebacterium glutamicum was extracted using phenol-chloroform extraction, as described in (ref. 22 DNA from Hirschia baltica was extracted using the Blood and Cell Culture DNA Maxi Kit (Qiagen, Valencia, CA, USA). DNA from Meiothermus silvanus, Nocardiopsis dassonvillei and Segniliparus rotundus was extracted using the Qiagen Genomic 500 DNA Kit (Qiagen, Hilden, Germany). DNA from Pseudomonas stutzeri was isolated using the Wizard Genomic DNA Purification Kit (Promega Corp., Madison, WI, USA). DNA from Coraliomargarita akajimensis, Halovivax ruber, Natronobacterium gregoryi and Spirochaeta smaragdinae was extracted using the Masterpure Gram-Positive DNA Purification Kit (Epicentre, Madison, WI, USA). All DNA extracts were quantified using the PicoGreen assay and the Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA) ( Supplementary Fig. 1). Each sample was quantified in quadruplicate. Samples were pooled at varying ratios to generate the mock community (Fig. 2, Supplementary Table 1).

Library creation and sequencing
For Illumina library creation, 100 ng of genomic DNA of MBARC-26, brought up to a total of 100 μl using TE, was sheared to 300 bp using the Covaris LE200 (Covaris, Inc., Woburn, MA, USA) and size- selected using AMPure XP beads (Beckman Coulter, Brea, CA, USA): 60 μl of beads were added to 100 μl of sample. The sample was then incubated at room temperature (RT) for 5 min. Beads were pelleted using a magnetic particle concentrator (MPC) (Thermo Fisher Scientific, South San Francisco, CA, USA) until liquid was clear. The supernatant was removed and transferred to a new tube. 30 μl of AMPure XP beads were then added for the second bead size selection. The mixture was pulse vortexed, quickly spun and incubated at RT for 5 min. Beads were pelleted using a magnetic particle concentrator (MPC) (Thermo Fisher Scientific, South San Francisco, CA, USA) until liquid was clear. The supernatant was then discarded without disturbing the beads and 200 μl of freshly prepared 75% ethanol (EtOH) was added, followed by a 30 s incubation to wash the beads. EtOH was discarded before the wash step with EtOH was repeated for a total of two washes. Afterwards, the sample was placed on a thermocycler (Eppendorf, Hamburg, Germany) with the lid open and incubated at 37°C until the beads were dry and residual EtOH had evaporated. The beads were re-suspended in 53 μl of EB buffer (Qiagen, Redwood City, CA, USA), vortexed, quickly spun and incubated at RT for 1 min. Beads were pelleted using an MPC until liquid was clear (Thermo Fisher Scientific, South San Francisco, CA, USA). 50 μl of supernatant was then transferred to a new tube. The DNA fragment size was assessed using the Agilent Bioanalyzer 2100 High Sensitivity Kit (Agilent Technologies, Palo Alto, CA, USA) before proceeding to end repair. The fragments were treated with the Kapa Library Preparation Kit ORIGIN (Kapa Biosystems, Wilmington, MA, USA) for the following steps: For end-repair 26 μl MilliQ water, 9 μl 10X End Repair Buffer, and 5 μl End Repair Enzyme were combined in a 1.5 ml tube. The cocktail was vortexed and quickly spun, then stored on ice. 40 μl of End Repair cocktail was added to the 50 μl DNA sample. The mixture was vortexed and quickly spun, before incubation at 30°C for 30 min in a thermocycler (Eppendorf, Hamburg, Germany). After incubation, 126 μl of AMPure XP beads (Beckman Coulter, Brea, CA, USA) were added to 90 μl of End Repair sample, pulse vortexed, quickly spun, and incubated at RT for 5 min. Beads were pelleted using a MPC until liquid was clear. The supernatant was then discarded without disturbing the beads. The beads were washed twice with 200 μl of freshly prepared 75% EtOH with an incubation time of 30 s. After washing, the sample was incubated at 37°C in a thermocycler with the lid open until residual EtOH had evaporated. For DNA elution, 17.5 μl of EB buffer was added. The sample was vortexed, quickly spun, and incubated at RT for 1 min, before beads were pelleted on a MPC. 15 μl of supernatant was then transferred to a new tube. For A-tailing, 9 μl of MilliQ water, 3 μl of 10X A-Tailing Buffer and 3 μl of A-Tailing Enzyme were combined in this order in a 1.5 ml tube. The cocktail was vortexed and quickly spun. 15 μl of the A-Tailing cocktail was added to the 15 μl sample. The mixture was vortexed and quickly spun. The samples were then incubated in a thermocycler at 30°C for 30 min, followed by 5 min at 70°C.
Adaptor ligation was immediately performed thereafter: 9 μl of 5X Ligation Buffer and 5 μl of ligase were combined in a 1.5 ml tube, vortexed and spun. The mixture was pulse vortexed and quickly spun. 14 μl of adaptor ligation cocktail were added to the 30 μl sample, before 1 μl of 18 μM adaptor was added to the ligation mixture for a final concentration of 400 nM. The mixture was incubated in a thermocycler at 20°C for 15 min.
After adaptor ligation, 5 μl of EB Buffer was added to 45 μl of adaptor-ligated sample. The sample was size-selected and washed twice with 45 μl of AMPure XP beads as described previously. After the first clean-up step, the sample was eluted with 52 μl of EB Buffer and 45 μl of supernatant was transferred to a clean tube. After the second clean-up step, the sample was eluted with 25 μl of EB Buffer. 23 μl of supernatant was transferred to a clean tube. The sample was quality-controlled and quantified using an Agilent Bioanalyzer 2100 High Sensitivity Kit.
The prepared Illumina library was further quantified by using the Kapa Biosystems next-generation sequencing library qPCR kit according to the manufacturer's guidelines (Kapa Biosystems, Wilmington, MA, USA). The amplification products were run on a Roche LightCycler 480 real-time PCR instrument  . The observed discrepancy between molarity and % mapped PacBio and Illumina sequences in T. composti is likely due to contamination as T. composti was previously found to occur as laboratory contaminant in various shotgun metagenome datasets (unpublished data). The smaller discrepancies are expected due to DNA quantification spreads and platform biases. Colors denote phylum association as defined in Fig. 1 Table 2). This resulted in 355,875,608 raw reads. For PacBio library creation, 5 μg of gDNA was sheared using a Covaris LE220 to generate 2 kb fragments (Covaris, Inc., Woburn, MA, USA). The sheared DNA fragments were then prepared according to the SMRTbell template preparation kit guidelines (Pacific Biosciences, Menlo Park, CA, USA). Briefly, DNA fragments were treated with DNA damage repair mix, end-repaired, and 5' phosphorylated. PacBio hairpin adapters were then ligated to the fragments to create SMRTbell template for sequencing. The SMRTbell templates were purified using exonuclease treatments and size-selected using AMPure PB beads (Pacific Biosciences, Menlo Park, CA, USA) ( Table 2).
Sequencing primers were annealed and v. P4 sequencing polymerase was bound to the SMRTbell templates. The prepared SMRTbell template libraries were then sequenced on a Pacific Biosciences RSII sequencer using v. C2 chemistry and 1x120 min sequencing movie run times (Pacific Biosciences, Menlo Park, CA, USA). This resulted in 300,584 raw reads (Table 2).

Sequence QC
Illumina shotgun reads were filtered using BBDuk (filterk = 27, trimk = 27; http://jgi.doe.gov/data-andtools/bb-tools/) to remove Illumina adapters, known Illumina artifacts, phiX, and to quality-trim both ends to Q12. Resulting reads containing more than one 'N', or with quality scores (before trimming) averaging less than 8 over the read, or length under 40 bp after trimming, were discarded. Remaining reads were mapped to a masked version of human HG19, dog, cat, and mouse with BBMap (http://jgi. doe.gov/data-and-tools/bb-tools/), discarding all hits exceeding 93% identity. This resulted in 347,963,988 filtered reads with average insert size of 219 ± 43 bp.
Quality filtering and error correction of PacBio sequences was performed using the RS_ReadsOfInsert protocol v. 2.3.0 in SMRT Portal (minimum subread length: 50 bp; minimum read quality: 75%). This resulted in 53,654 quality-filtered subreads with average read length of 1,041 ± 576 bp.

Mapping, repeat regions, and phylogenetic tree construction
High quality Illumina and PacBio sequences were mapped to their bacterial and archaeal reference genomes using BBMap with parameters bbmap.sh, ambig = toss (Illumina) and mapPacBio.sh, ambig = toss (PacBio), respectively. Numbers of mapped sequences were normalized to the respective whole genome and chromosome lengths of reference organisms (Supplementary Table 1). Unmapped sequences amounted to 2,105 (3.92%) and 3,777 (7.04%) PacBio sequences, when mapped against genome and chromosome references, respectively. In the Illumina dataset, 8,981,844 (2.58%) and 18,088,260 (5.20%) Illumina sequences remained unmapped, when mapped against genome and chromosome reference, respectively. Repeat regions reported here were retrieved from NCBI GenBank 23 on May 16, 2016. They include tandem, inverted, flanking, terminal, direct and dispersed repeat types. For phylogenetic tree construction, full-length 16S rRNA gene sequences were aligned using the SINA aligner 24 including 10 neighbors at 95% minimum identity for classification against the SILVA, RDP, greengenes, LTP, and EMBL databases 25 . The alignment was masked using the SILVA-compatible 1,349 Lane mask 26 . Tree construction was performed using FastTree 27 .

Data Records
Filtered shotgun sequences generated on the Illumina and PacBio platforms are publically available through NCBI (Data Citation 1 and Data Citation 2).

Technical Validation
To assess the quality of genomic DNA received, we used the PicoGreen assay and the Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA). Each sample was quantified in quadruplicate. Samples were pooled at varying ratios to generate the mock community ( Supplementary Fig. 1).