extfieldsofstudy: Medicine Biology

text:
Metagenomic Analysis of Flaviviridae in Mosquito Viromes Isolated From Yunnan Province in China Reveals Genes From Dengue and Zika Viruses

More than 6,000 mosquitoes of six species from six sites were collected and tested for their virome using metagenomics sequencing and bioinformatic analysis. The identified viral sequences belonged to more than 50 viral families. The results were verified by PCR of selected viruses in all mosquitoes, followed by phylogenetic analysis. In the present study, we identified the partial dengue virus (DENV), Zika virus (ZIKV), and Japanese encephalitis virus (JEV) sequences in mosquitoes. Metagenomic analysis and the PCR amplification revealed three DENV sequences, one of which encodes a partial envelope protein. Two ZIKV sequences both encoding partial nonstructural protein 3 and one JEV sequence encoding the complete envelope protein were identified. There was variability in the viral titers of the newly isolated virus JEV-China/YN2016-1 of different passage viruses. The newly identified Zika virus gene from ZIKV-China/YN2016-1 was an Asian genotype and shared the highest nucleotide sequence identity (97.1%) with a ZIKV sequence from Thailand isolated in 2004. Phylogenetic analysis of ZIKV-China/YN2016-1 and ZIKV-China/YN2016-2 with known Flavivirus genes indicated that ZIKV has propagated in Yunnan province, China.

INTRODUCTION
Mosquitoes are important hosts and insect vectors for a range of infectious viruses including Japanese encephalitis virus (JEV), Zika virus (ZIKV), dengue virus (DENV), and Getah virus (GETV), which pose a significant threat to human health and represent an economic burden worldwide (Pham et al., 2005;Klungthong et al., 2007). Mosquitoes acquire viruses when consuming blood from hosts undergoing viremia. The viruses then replicate and propagate in the insect host before being introduced into further victims during biting and blood feeding (Ritchie et al., 2007;Motooka et al., 2015). Mosquitoes feed on a wide range of hosts including humans and other vertebrates, invertebrates, and plants (Shi et al., 2014). Yunnan province in China harbors a diverse range of mosquito-borne viruses (Feng et al., 2015); therefore, regional surveillance is imperative. The detection of viruses in mosquitoes is usually performed by using reverse transcription polymerase chain reaction (RT-PCR) and nested PCR approaches (Almeida et al., 2005;Houghton, 2009;Li et al., 2015;Sim et al., 2015). However, compared with Illumina sequencing (Alquezar-Planas et al., 2013;Cholleti et al., 2016;Ergünay et al., 2017), these traditional methods are time-consuming and labor-intensive to detect low-level viromes in mosquito vectors (He et al., 2013;Miesen et al., 2016). Therefore, metagenomic analysis of mosquitoes is likely to be of great value to avoid missing the detection of viruses with high infectivity and pathogenicity as well as the detection of previously unknown viruses.
The present study aimed to build a valid surveillance method to monitor the distribution of Flaviviridae from mosquitoes in Yunnan province, China and to provide useful insights into viral isolation, prevention, and control. Diverse and abundant viromes from mosquitoes isolated in Yunnan province were investigated using metagenomic sequencing and nested PCR. The presence of DENV, JEV, and ZIKV viruses was confirmed, and JEV was isolated. This preliminary exploration of the metagenomes of mosquito-borne viruses lays the foundation for further research on the territorial distribution, diversity, and surveillance of mosquito-borne viruses in China and other countries.

Mosquito Collection
More than 6,000 living or freshly dead female mosquitoes were collected in Yunnan province, China, during August and September 2016 (Figure 1). Sample I comprised Culex tritaeniorhynchus (C. tritaeniorhynchus). Sample II was a mixture of Armigeres obturbans (Ar. obturbans) and Aedes albopictus (Ae. albopictus). Sample III was a collection of Anopheles sinensis (An. sinensis), Uranotaenia hebes (U. hebes), and Armigeres durhami (Ar. durhami) ( Table 1). The mosquito samples were grouped based on species. C. tritaeniorhynchus was the predominant transmitting vector of JEV, and Ar. obturbans and Ae. albopictus were the predominant transmitting vectors of DENV and ZIKV (Byrd et al., 2013). The identified species of mosquitoes were collected separately and stored at −80 • C. This research was approved by the Ethics Committee and the Research Ethics Committee of Yanbian University Medical College. All experiments involving active viruses were performed in a biosafety level 3 laboratory.

Preparation of Mosquito Samples
Briefly, each mosquito sample was mixed and ground using a tissue grinder with SM buffer (50 mM Tris, 10 mM MgSO4, 0.1 M NaCl, pH 7.5) (MeilunBio, Dalian, China). To remove mosquito debris and other materials, the mixed samples were centrifuged at 13,000× g for 20 min, and the supernatants were used to extract the viral nucleic acid.

Extraction of Viral Nucleic Acid and Reverse Transcription
To remove the free nucleic acid and eliminate the contaminating host genomic DNA, 14 U Turbo DNase (Ambion, Austin, TX, USA), 25 U Benzonase Nuclease (Novagen, San Diego, CA, USA), 20 U RNase I (Fermentas, Ontario, Canada), and 10× DNase buffer (Ambion) were added to 127 µL of the supernatants to a final volume of 150 µL, followed by digestion at 37 • C for 1 h. The viral nucleic acid in the obtained products was extracted using a virus nucleic acid isolation kit (Bioer Technology, Hangzhou, China) according to the manufacturer's instructions. The total viral nucleic acids were reverse-transcribed using anchored random primers and Superscript III reverse transcriptase (Invitrogen, Carlsbad, CA, USA). The anchored random primers ( Table 2) were added separately to the viral nucleic acid and incubated at 75 • C for 5 min and then placed on ice for 5 min for denaturation. To obtain the reverse-transcribed product, 40 U of RNase OUT (Invitrogen, Carlsbad, USA), 200 U of SuperScript III reverse transcriptase (Invitrogen), 1 µL of 0.1 M dithiothreitol (DTT) (Invitrogen), 1 µL of 10 mM dNTPs (TaKaRa, Dalian China), 4 µL of 5× first-strand buffer (Invitrogen), and RNasefree H 2 O (TaKaRa) were added to a final volume of 20 µL and incubated at 25 • C for 10 min, followed by 50 • C for 60 min and then 75 • C for 10 min.

Synthesis of Double-Strand cDNA (dscDNA)
Before the synthesis of dscDNA, 1 µL of RNase H (TaKaRa) was added to the obtained products to degrade the free RNA. To synthesize dscDNA, anchored random primers were added and incubated 75 • C for 5 min and then on ice for 5 min for denaturation. Then, 1 µL of Klenow fragment (TaKaRa), 1 µL of 10 mM dNTPs (TaKaRa), 2 µL of 10× Klenow buffer (TaKaRa), and 6 µL of dd H 2 O (TaKaRa) were added and the samples were incubated at 37 • C for 60 min, followed by 75 • C for 10 min. Then, 0.5 µL of exonuclease I (TaKaRa), 1 µL of shrimp alkaline phosphatase (SAP, TaKaRa, Dalian, China), 5 µL of 10× phosphatase buffer (TaKaRa), and 24 µL of DEPC H 2 O (TaKaRa) were added and incubated at 37 • C for 60 min followed by 75 • C for 10 min to eliminate the phosphates and the free single-strand nucleic acid in the dscDNA reaction.

Sequence-Independent Single-Primer Amplification (SISPA) and Purification of PCR Products
To increase the quantity of the viral nucleic acids, SISPA was applied to amplify the dscDNA. The 50 µL reaction comprised 10 µL of dscDNA, 2 µL of barcode primer (

Metaviral Sequencing
The purified PCR products from the three samples were sent to the Wuhan Genome Institute (BGI, Shenzhen, China) for Illumina sequencing. Briefly, to obtain ∼180 bp DNA fragments, the purified PCR products were ultrasonicated, and dATP and Klenow fragment were added to produce 3 ′ -dA overhangs. To establish genomic DNA libraries, DNA fragments were bound to Illumina adaptors and amplified using PCR with adaptor primers. Amplicons were ligated to flow cells to which fluorescently labeled dNTPs were added. DNA sequences were identified by using the sequencing-by-synthesis method (SBS, Illumina). Base calling was monitored by the program GAPipeline (BGI), with default settings. No-calling reads and adaptor sequences were eliminated. The remaining sequences were assembled into contigs using SOAPdenovo software (BGI, Shenzhen, China). Contigs and sequences longer than 100 bp were defined as significant data for further in silico analysis.

Computational Analysis
The contigs and sequences were aligned using the blastx and blastn with the nonredundant and viral reference sequences in the GenBank database. (https://www.ncbi.nlm.nih.gov/ genbank/) BLAST hits with an E-value ≤ 10e −5 were considered significant. After eliminating the bacterial and eukaryotic sequences, we analyzed the virus-like sequences.

Identification of Detected Viruses
Based on the alignment outcomes of the viral contigs and the match position of the viral contigs with the corresponding viruses in GenBank, specific primers were designed and synthesized to identify the detected viruses. Viral nucleic acids were extracted by Bioer Technology (Hangzhou, China) and amplified using the designed primers and a PCR Master Mix (Tiangen, Beijing, China). To avoid false positives in Illumina sequencing, viral identification was performed three times in an independent manner.

Phylogenetic Analysis
The obtained PCR products were sequenced, and the sequences were aligned against sequences of the representative viruses using CLUSTAL W version 2.0. (Accession numbers are shown in the phylogenetic trees). To achieve better accuracy in phylogenetic analysis, neighbor-joining phylogenetic trees were produced using MEGA 7 with 1,000 bootstrap replicates Li, 2015).

Isolation of Viruses
Viral isolation was conducted with the positive ground mosquito supernatants. Briefly, the ground mosquito supernatants were diluted seven-fold with DMEM containing 2% FBS and cultured with BHK-21 cells for 5-7 days. Cultures were examined daily for evidence of a virally induced cytopathic effect (CPE). Cultures without a CPE were blind-passaged three times.
Identification of the Antiviral Effect Using an Indirect Immunofluorescence Assay (IFA) BHK-21 cells were seeded in 12-well plates, and after 24 h of culture, they were rinsed with phosphate-buffered saline (PBS) (Solarbio, Beijing, China) and fixed using methanol (Solarbio) for 10 min at room temperature. The cells were washed three times with PBS and then incubated in blocking buffer [3% bovine serum albumin (BSA) (BOSTER, Beijing, China) in Trisbuffered saline (Solarbio) and Tween20 (Solarbio)] for 1 h at room temperature. The cells were then incubated with an anti-E monoclonal antibody (Abcam, Cambridge, UK) (1:20) in PBS with 3% BSA overnight at 4 • C. After the cells were washed three times in PBS, they were incubated with a fluorescently labeled secondary antibody-fluorescein isothiocyanate-conjugated goat antimouse antibody (ZSGB-Bio, Beijing, China)-in the dark at 37 • C for 1 h. The nuclei were washed with PBS three times and stained with 4,6-diamidino-2-phenylindole (DAPI) (Solarbio). Slides were imaged under a fluorescence microscope (Olympus, Tokyo, Japan).

Observation by Negative-Stain Electron Microscopy
BHK-21 cells with suspected JEV-induced CPE were collected by repeated freezing and thawing three times, followed by centrifugation at 10,000× g for 15 min. The obtained virus suspension was centrifuged at 60,000× g for 4 h and the supernatant was removed gently and discarded. Subsequently, the viral precipitate was resuspended with an isometric mixture of 6.1% (v/v) pH 7.2 glutaraldehyde (HEDEBIO, Beijing, China) fixative and DMEM, 25 µL of which were added to the copper grid. After desiccation, one drop of 3% phosphotungstic acid (JINDU, Shanghai, China) was added for negative staining.
Before observation under an electron microscope (FEI, Hillsboro, USA), the grid was placed in an incubator (SANYO, Osaka, Japan) at 37 • C for desiccation.

Amplification and Phylogenetic Analysis of JEV Complete Genome
The JEV strain was isolated, and then the complete sequence was amplified using RT-PCR. The designed primers are displayed in Table 3. The complete JEV genome was assembled by PCR products using SeqMan version 7.1 (DNAStar, Madison, USA) with the default value. Subsequently, phylogenetic analysis was performed based on the complete sequence of JEV with other representative JEV strains reported in China and neighboring countries.

Detection of the Viral Titer of Different Passages
The viral titers of different passages of JEV (passage 5, passage 10, passage 20, and passage 30) were evaluated in BHK-21 cells in triplicate, with DMEM as a control. Cells were incubated with serial ten-fold dilutions of JEV or DMEM for 120 h at 37 • C. The viral titer TCID 50 (50% tissue culture infective dose) was calculated by using the Reed-Muench method (Krah, 1991), according to the cytopathic effect (CPE) in cells caused by viruses. It was evaluated based on the changes in cell culture morphology under a light microscope (OLYMPUS, Tokyo, Japan) that were observed at 100× magnification. The CPE of JEV or DMEM at each dilution was evaluated in eight replicates.

Variability of Envelope (E) Genes of Different Passage Viruses
Viral RNA was extracted from the supernatant using an RNA viral kit according to the manufacturer's protocol (Bioer Technology) after a CPE was induced in BHK-21 cells by different passage viruses (P5, P10, P20, and P30). The purified RNA was used as the template for cDNA synthesis using the SuperScript TM III first-strand synthesis system (Invitrogen) with the reverse primer JWR1 5 ′ -AGATCCTGTGTTCTTCCTCACCACCA-3 ′ , according to the manufacturer's instructions. The envelope (E) genes of different passage viruses were amplified using the paired primers shown in Table 4 and sequenced. The alignment of the four E genes was performed using MEGA 7.0.

Statistical Analysis
The data were analyzed using SPSS 17.0 software. Group comparisons for viral titers were carried out by using the t-test for independent means. Experiments were performed in triplicate and a minimum of three independent experiments were evaluated. The value of P < 0.05 was considered statistically significant.

GenBank Accession Numbers
The

Metaviral Sequencing and the Virome of Mosquitoes
To obtain the clean data of the virome of the mosquitoes, contaminating host sequences and barcode DNA were eliminated. A total of 22,964,305 reads were acquired by Illumina sequencing, with averaging read lengths of 186 nt ( Table 1). The amount of viral sequences was 17.00, 0.24, and 13.85% in sample I, II, and III respectively and the viruses could be clearly classified (Figure 2). Among them, 87.4% were  vertebrate viruses, 5.4% were plant viruses, and 1% were fungal viruses; insect viruses and phages accounted for 0.8 and 0.6%, respectively (Figure 2A). Most viral sequences in Sample I were from C. tritaeniorhynchus, the predominant JEV-transmitting vector in China, followed by A. obturbans. A. albopictus is the predominant ZIKV-transmitting vector. Viral sequences were classified at the family level ( Figure 2B), and Flaviviridae was the predominant viral family, together with certain other important viral families such as Circoviridae, Anelloviridae, Parvoviridae, and Peribunyaviridea. Additionally, numerous unclassified viral sequences, presumably belonging to novel viruses not previously identified in mosquitoes, were identified and are worthy of further study.

PCR Verification of the Metavirome Results
Viral sequences were assembled into contigs using SOAPdenovo. From the 13,645 assembled viral contigs, 32 DENV-like contigs with a read coverage of 184× (236-1408 nt), 13 ZIKVlike contigs with a read coverage of 34× (194-652 nt), and 26 JEV-like contigs with a read coverage of 52× (306-1,564 nt) were obtained. The DENV-like contigs shared 92-98% nt identity with known DENV sequences. The ZIKVlike contigs shared 95-99% nt identity with known ZIKV sequences. The JEV-like contigs shared 96-99% nt identity with known JEV sequences. To confirm the outcome of metavirome sequencing, specific primers were designed and synthesized to amplify the identified viruses ( Table 4). The following were the results of PCR identification in mosquito samples.

PCR Amplification of Dengue Virus
Dengue virus belongs to Flaviviridae, and 3,044 cases of dengue were reported in China in 2005-2012, among which 134 cases were reported in Yunnan province. Currently, there is no DENV-5 case reported in China. There have been DENV 1-4 outbreaks in China repeatedly and alternately. However, DENV-1 is predominant. (Xiong and Chen, 2014). Based on the alignment results of viral contigs, viral PCR amplicons from samples II shared a high identity with the dengue virus (DENV) genes, indicating that they belong to insect viral species. The results were verified by nested RT-PCR, which identified a 1,231 nt sequence (DENV-China/YN2016-1) encoding a partial nonstructural protein 4A, 4B, and 5 gene; a 491 nt sequence (DENV-China/YN2016-2) encoding a partial E protein, and a 939 nt sequence (DENV-China/YN2016-3) encoding a partial nonstructural protein 5 gene from DENV. Phylogenetic analysis showed that the newly identified DENV sequences were genotype IV. In addition, DENV-China/YN2016-2 was highly homologous to a DENV sequence from Thailand that was isolated in 1991, sharing 95.3% nt identity. This suggested that the viruses represented a DENV variant ( Figure 3A). Interestingly, the partial E gene in DENV-China/YN2016-2 may represent a new genotype IV gene in the DENV variant strain.

PCR Amplification of Zika Virus
The first imported case of ZIKV, pertaining to Flaviviridae, was confirmed in China in February 2016 Zhang et al., 2016). A total of 28 cases were reported in Zhejiang province, China in September 2016 (Deng et al., 2016;Li et al., 2016). In addition, 19 cases were reported in Guangdong province, China in 2016 (Sun et al., 2017). In line with the alignment results of viral contigs, viral PCR amplicons from sample II shared a high identity with ZIKV, also belonging to insect viral species. These results were also confirmed by means of nested RT-PCR, which identified 694 and 470 nt sequences, both encoding partial nonstructural protein 3 from ZIKV. These sequences were named ZIKV-China/YN2016-1 and ZIKV-China/YN2016-2, respectively. Phylogenetic analysis showed that the newly identified Zika virus sequence ZIKV-China/YN2016-1 was from an Asian genotype. The nucleotide identity analysis revealed that it shared the highest nucleotide sequence identity (97.1%) with a ZIKV sequence from Thailand that was isolated in 2004 ( Figure 3B).

PCR Amplification of JEV
Since 1951, 33,900 cases of JEV, which belongs to Flaviviridae, were reported in China. The most obvious feature for the territory distribution of JE in China in the 1970s or 2010 is that a small number of cases were reported in west and north China, whereas a high number of cases occurred in east and southwest China. Yunnan province is still a highly epidemic area for JEV. First reported in 1949, JEV has existed in China for more than 60 years. The prevalence of JEV is genotype I and III in China, and genotype I is predominant in Yunnan province (Zheng et al., 2012). According to the alignment results of viral contigs, viral PCR amplicons from sample I shared a high identity with JEV genes, indicating that they belong to insect viral species. The results were verified by means of nested RT-PCR, which identified a 1,500 nt sequence (JEV-China/YN2016-1) encoding a complete E protein. Phylogenetic analysis showed that the newly identified JEV sequence, JEV-China/YN2016-1, was genotype I (Figure 4). Nucleotide identity analysis revealed that it shared the highest nucleotide sequence identity (98.8%) with a JEV sequence from Laos that was isolated in 2009.

Viral Identification of JEV and Phylogenetic Analysis of its Complete Genome
After inoculation into BHK-21 cells, the obtained JEV strain induced an obvious CPE. Monolayer cells appeared rounded, with bulging, aggregation, shedding, and incomplete destruction of a single loose layer, with gaps between cells. The severity of the CPE was time-dependent. JEV was propagated in BHK-21 cells for 24 h. The levels of intracellular virus were measured by IFA using an E-specific monoclonal antibody. As expected, mockinfected cells showed no E expression, while JEV-infected cells were positive for E expression. Viral particles were observed using negative-stain electron microscopy. The particles were rounded with a diameter of ∼30-40 nm. Moreover, there were tiny protrusions on the surface, similar to JEV particles ( Figure 5). The phylogenetic analysis of the complete JEV genome showed that JEV-China/YN2016-1 still shared the highest nucleotide sequence identity (97.6%) with a JEV sequence from Laos that was isolated in 2009 (Figure 6).

DISCUSSION
Mosquitoes are the intermediate hosts of numerous viruses that infect humans, animals, insects, plants, and other species, and play a significant role in the prevalence of many infectious diseases Fisher et al., 2015;Shi et al., 2015). Before 2007, virus identification by traditional virological methods required several decades and involved laboratories worldwide (Krah, 1991). Compared with the traditional method, metagenomic sequencing by Illumina sequencing combined with a high throughput analysis is more efficient, resulting in the identification of many viruses from different species in recent years by several laboratories (Wilson et al., 2015). Moreover, metagenomic sequencing has revealed the high abundance of viruses in mosquitoes (Bolling et al., 2015;Fauver et al., 2016;Shi et al., 2017). Metagenomic sequencing has also contributed to the discovery of novel viruses and their identification and characterization (Carissimo et al., 2016;Frey et al., 2016). In the present study, more than 50 viral families were detected by Illumina sequencing in mosquito samples. However, there was a considerable difference in the percentage of viral sequences obtained from the three mosquito groups. It could be related to the biology of the vectors (Atoni et al., 2018;Xia et al., 2018;FIGURE 6 | Phylogenetic trees of the complete JEV genome. Phylogenetic trees based on the complete genome of JEV. The trees were constructed using the p-distance-based neighbor-joining method in MEGA 7.0 software. Bootstrap values were calculated with 1,000 replicates. Black solid circles indicate the complete genome of JEV isolated in this study. Zakrzewski et al., 2018). Still, new findings of viral distribution and evolution in Yunnan province were obtained. RT-PCR and nested-PCR were used to confirm the results of Illumina sequencing. The detection of DENV in the mosquito samples showed that the virus still exists in Yunnan province, and many strains of DENV have been uncovered by previous studies Shi et al., 2014;Zhang et al., 2014;Wang et al., 2016). However, the detection of new DENV strains is important, as it suggests a new prevalence of DENV in this area. Moreover, the identification of the Zika virus genes is significant, suggesting that increased environmental surveillance is required because of its strong infectivity to humans, especially pregnant women and infants. Meanwhile, the presence of unidentified viruses in the Illumina sequencing results may ascribed to inadequate mosquito sampling and limited collection sites.
ZIKV and DENV were detected simultaneously, indicating that these viruses co-circulate in Yunnan province and suggesting that co-infection by the viruses is possible. Zika virus pathogenesis can be enhanced by preexisting dengue virus (Bardina et al., 2017). The simultaneous detection of the genes of two viruses may represent a new challenge to the control of ZIKV. We performed the isolation of DENV and ZIKV viruses in Vero and C6/36 cells besides BHK-21 cells as well as the suckling-mouse brain was also used in viral isolation. It was unfortunate that we could not isolate DENV and ZIKV.
The detection and isolation of JEV indicated that the virus still exists in Yunnan province. Analyzing the variant sites in the E gene and the consequences for the sequence of the encoded protein, we discovered that amino acid 138 in the E proteins of viruses at P5 and P30 is K (lysine). However, amino acid 138 in the E proteins of viruses at P10 and P20 is E (glutamic). Lysine at position 138 in the E protein is an attenuated viral site (Ni et al., 1994). The results showed that JEV-China/YN2016-1 has undergone mutations when passaged in BHK-21 cells, suggesting the possibility of increased threat to human health induced by JEV variation. The phylogenetic analysis based on the complete genome and the E gene of JEV-China/YN2016-1 showed that it shared the highest nucleotide sequence identity with a JEV sequence from Laos that was isolated in 2009, indicating the possible threat of JEV transboundary transmission.
This is the first report of the detection of Zika virus genes in Yunnan province, China. The newly identified DENV-China/YN2016-2 is a novel genotype IV gene in the DENV mutant strain. Metagenomic analysis of the viromes of mosquitoes in Yunnan province, combined with nested PCR, revealed that the distribution of viruses is dependent on both the mosquito species and the geographical location. The present work explored only a small portion of the virome of mosquitoes prevailing in this region of China. A much larger study on mosquitoes in other provinces in China and from other countries is needed to assess the diversity of mosquito-borne viruses. In conclusion, our findings provide a useful insight into the viral isolation and the characterization of Zika and dengue viruses and could hence be applied to larger studies in the future.

AUTHOR CONTRIBUTIONS
PX, HJL, and NJ conceived and designed the experiments. PX, JH, CL, SW, YL, MW, and NJ performed the experiments. PX, HL, JR, and NJ analyzed the data. YZ, XG, HZ, HJL, and MT contributed reagents, materials, and analysis tools. PX and NJ wrote the paper. All authors read and approved the final version of the manuscript.


text:
Impact of DNA Sequencing and Analysis Methods on 16S rRNA Gene Bacterial Community Analysis of Dairy Products

Validated methods are urgently needed to improve DNA sequence-based assessments of complex bacterial communities. In this study, we used 16S rRNA PCR amplicon and gDNA mock community standards, consisting of nine, dairy-associated bacterial species, to evaluate the most commonly applied 16S rRNA marker gene DNA sequencing and analysis platforms used in evaluating dairy and other bacterial habitats. Our results show that bacterial metataxonomic assessments are largely dependent on the DNA sequencing platform and read curation method used. DADA2 improved sequence annotation compared with QIIME 1, and when combined with the Ion Torrent PGM DNA sequencing platform and the Greengenes database for taxonomic assignment, the most accurate representation of the dairy mock community standards was reached. This approach will be useful for validating sample collection and DNA extraction methods and ultimately investigating bacterial population dynamics in milk- and dairy-associated environments.

This approach has enabled the simultaneous identification of the majority of bacteria in complex microbial communities. Although analysis of 16S rRNA gene diversity has provided significant new perspectives on bacterial habitats, there remain challenges to sample preparation, DNA sequencing, and data analysis approaches for ensuring accurate measurements of bacterial populations.
DNA sequencing platforms, including 454 pyrosequencing, Illumina, Ion Torrent, and Pacific Biosciences have also been shown to cause variation in bacterial community assessments (4,(18)(19)(20)(21). Moreover, data analysis methods, especially read clustering approaches (i.e., generating representative sequences), are known to have a significant impact on the interpretation of bacterial composition (21)(22)(23)(24)(25)(26)(27). De novo sequence clustering can result in unstable operational taxonomic units (OTUs) between projects that are composed of different sequences with each clustering iteration (28,29). Reference-based sequence clustering tends to result in fewer sequence variants than de novo methods (21,23), but can still lead to overestimation of bacterial community diversity caused by insufficient read quality control and error filtering (27). As a result, there is now an effort to move away from OTU-based methods toward DNA sequences that represent single nucleotide variation (30,31). One of the amplicon sequence variant (ASV) clustering methods is DADA2 (Divisive Amplicon Denoising Algorithm 2), which builds a quality-based model for filtering error and identifying variation in 16S rRNA gene sequences (26).
Herein, we sought to compare different DNA sequencing, read assembly, and data analysis strategies for the capacity to accurately detect the composition of a mock bacterial community consisting of nine species commonly found in milk. To eliminate biases introduced by sample type and DNA extraction method and focus on close examination of the biases introduced by PCR, sequencing, and bioinformatics analyses, we employed two different mock communities consisting of either organism-specific PCR amplicons or purified genomic DNA (gDNA) (Fig. 1). This approach allowed us to compare the performance of two popular benchtop DNA sequencers (Illumina MiSeq and Ion Torrent PGM), paired-end read assembly of Illumina MiSeq, OTU (QIIME 1 open reference)/ASV (DADA2) analysis methods, and reference taxonomy databases (Greengenes and RDP). Our results showed that the combination of DADA2 and the Greengenes analysis pipeline, paired with Ion Torrent PGM sequencing, results in the most accurate representation of the mock communities.

RESULTS
Comparison of representative sequence analysis of 16S rRNA V4 region reads generated with the Illumina MiSeq and Ion Torrent PGM. A mock community was prepared by combining 16S rRNA V4 region PCR amplicons from nine bacterial strains (Table 1) in equimolar quantities prior to DNA sequencing on the Ion Torrent PGM and Illumina MiSeq instruments. Sequences were either assembled (Illumina MiSeq) or maintained as single-end reads (Illumina MiSeq and Ion Torrent PGM). A low percentage of reads were identified as chimeras (0 to 1.4%) (see Table S1 in the supplemental material), and the remaining reads were analyzed using QIIME 1 (UCLUST) following the open-reference pipeline at a 97% threshold or DADA2 pipelines for OTU or ASV identification using the Greengenes (version 13.8) and RDP (version GOLD for QIIME 1 and version 11.5 for DADA2) reference databases. Total reads after quality filtering for each sequencing and read assembly method are shown in Table S1.
Illumina MiSeq paired-end assemblies. QIIME 1 analysis of Illumina MiSeq pairedend assembled reads with recommended parameters (32) resulted in total OTU numbers that were at least 4.2-fold greater than the expected nine OTUs encompassing strains included in the mock community (Table 2). When the Greengenes database was used for OTU alignment, 85 OTUs were identified. The majority of those OTUs (i.e., 65) were assigned to taxa included in the mock community. Although the numbers of OTUs varied for each taxon, a single OTU representative encompassed the majority (Ͼ 60%) of reads for each of the mock community members (Table 2). For example, out of nine Staphylococcus OTUs identified with Greengenes, 99% of the reads were represented by one OTU. The remaining 20 OTUs identified with Greengenes were either designated as taxa that were not included in the mock community or were designated only to the order level. When the RDP database was used as the reference database for QIIME 1 analysis, the total OTUs decreased to 70, despite the increase in spurious "other"  assignments (Table 2). Nevertheless, the overall OTU number remained considerably higher than the expected nine OTUs based on the mock community composition. The numbers of ASVs identified with DADA2 from Illumina paired-end assemblies were lower than the numbers of OTUs assigned with QIIME 1 ( Table 2). Because DADA2 assigns ASVs independently from taxonomic reference databases, total ASV numbers were the same using both Greengenes and RDP. The only distinction was that one Clostridiaceae ASV and both Escherichia ASVs identified using Greengenes were designated as Clostridium and Enterobacteriaceae in RDP (Table 2).
Illumina MiSeq unassembled single-end reads. Without read assembly, the QIIME 1 pipeline resulted in 68 and 36 total OTUs with the Greengenes and RDP databases, respectively (Table 3). These OTU numbers were lower than the paired-end assemblies ( Table 2). DADA2, on the other hand, resulted in a slightly higher number of ASVs (40 ASVs) than the paired-end assemblies (38 ASVs) (Tables 2 and 3). More reads were regarded as "other" taxa, and this result was most likely due to the shorter lengths of  the unassembled, single-end MiSeq reads. Between the two reference databases, Greengenes resulted in more accurate taxonomic assignments with DADA2. Two Bacillales ASVs and one Clostridiales ASV that were included in the "other" ASV category by RDP were assigned as Bacillus and Clostridiaceae with Greengenes. Moreover, the Escherichia/Shigella ASV in RDP was unambiguously allotted to the Escherichia genus by Greengenes (Table 3).
Ion Torrent PGM reads. The application of QIIME 1 with the Greengenes database to the Ion Torrent reads resulted in the highest number of OTUs out of any of the methods applied (Table 4). The use of RDP in QIIME 1 also yielded high OTU numbers, comparable to those found for the paired-end Illumina MiSeq assemblies (Table 4). Conversely, DADA2 resulted in only 13 ASVs (Table 4). The 4 additional ASVs compared to the expected 9 ASVs were created due to errors in the homopolymer regions (see Fig. S1 in the supplemental material), and the 13 ASVs were distributed across the 9 bacterial taxa included in the mock community, with the exception of one ASV with ambiguous taxonomy (Bacillales) identified using RDP, which was identified as Bacillus using Greengenes. Greengenes also improved the assignment of Escherichia coli to the genus level, as opposed to the family level in RDP (Table 4).
For each of the three DNA sequencing/read curation methods tested, DADA2 assigned fewer ASVs per taxon and resulted in fewer spurious ASVs than QIIME 1 (UCLUST) assigned OTUs, except in Illumina single-end results analyzed with the RDP database. DADA2 taxonomic identification was more specific with the Greengenes than the RDP database. Therefore, the combined DADA2/Greengenes approach was used for the subsequent analyses described below.
Assessments of the gDNA mock community were altered depending on DNA sequencing platform. A gDNA mock community was prepared by mixing equal quantities of gDNA from the nine milk-associated bacterial species prior to barcoded 16S rRNA V4 region PCR amplification (Fig. 1). The PCR products were then used for sequencing on either the Illumina MiSeq or Ion Torrent PGM, followed by analysis with the DADA2/Greengenes method. More chimeras were found for the gDNA mock community (ranging from 0.3 to 4.6%) than the PCR amplicon mock community (Table S1), indicating amplification errors arose from multitemplate PCR. However, except for the known variation in platform-dependent read lengths, nucleotide sequences of the most abundant ASVs assigned to each of the nine mock community species were identical between the Illumina MiSeq (single and paired ends) and Ion Torrent PGM platforms (see Illumina MiSeq paired-end assembly in Fig. S2, Illumina MiSeq single-end reads in Fig. S3, and Ion Torrent PGM reads in Fig. S4 in the supplemental material). Nucleotide sequence alignments of those ASVs to the corresponding ASVs identified from the PCR amplicon mock community also showed 100% nucleotide sequence conservation ( Fig. S2 to S4). For both Illumina MiSeq assembled and unassembled (single-end) reads, the gDNA mock community resulted in high numbers of ASVs (Table 5). These numbers were higher than those found for the PCR amplicon mock community (Tables 2 and 3) and were primarily due to the higher quantities of spurious ASVs (e.g., Clostridiales, Lactobacillus, and Oscillospira) present at low proportions (0.02 to 3.85% of total reads for each ASV) (see Table S2 in the supplemental material). As a result, the Shannon index of the gDNA mock community was elevated compared to the PCR amplicon mock community for both paired-end and single-end Illumina MiSeq results (Fig. 2), and these values were significantly increased compared to the expected ␣ diversity based on mock community composition. Interestingly, the same number of 13 ASVs was found for the gDNA and PCR amplicon mock communities when the Ion Torrent PGM was used (Table 5), and the Shannon index of the gDNA mock community resembled the PCR amplicon mock community expected value (Fig. 2).
Ion Torrent PGM sequencing with the DADA2/Greengenes method resulted in more accurate representations of the gDNA and PCR amplicon mock communities. DNA sequencing approaches were next compared for their capacity to yield the expected ␤ diversity and proportions of bacterial taxa included in the two mock communities. According to UPGMA (unweighted pair group method using average linkages) hierarchical clustering of Bray-Curtis dissimilarity metrics, results from the three DNA sequencing approaches (Illumina MiSeq paired-end assembly, single-end, and Ion Torrent PGM) were all in different clusters compared to the expected bacterial composition, independent of the gDNA or 16S rRNA PCR amplicon community type (Fig. 3A). Conversely, UPGMA of the weighted Unifrac distance metrics clustered the sequences according to mock community type (Fig. 3B). These comparisons showed that the gDNA mock communities sequenced with the Ion Torrent PGM were the most similar to theoretical (expected) proportions. No single method was found best suited for representing the PCR amplicon mock community (Fig. 3B). To assess whether the use of DADA2/Greengenes influenced this outcome, the other data analysis methods were compared, and it was found that the DNA sequencing platform used was consistently influential on mock community ␤ diversity (e.g., QIIME 1 with the Greengenes database is shown in Fig. S5 in the supplemental material). Examination of the relative abundances of individual taxa across the three DNA sequencing approaches showed that for the 16S rRNA PCR amplicon mock community, the proportions of most bacterial species were mostly not significantly altered compared to expected theoretical values. Exceptions to this finding were the reduced proportions of Enterococcaceae and Enterococcus found for Illumina paired-end assemblies and Streptococcus for both Illumina MiSeq methods as well as the Ion Torrent PGM platform (Fig. 4). For the gDNA mock community, the proportions of Escherichia and Streptococcus were significantly different from expected for all three DNA sequencing platforms. The proportions of Pseudomonas were also significantly lower than expected for the single-and paired-end assemblies from the Illumina MiSeq, and the proportions of Bacillus and Lactococcus were also altered for the paired-end assemblies. Lastly, there were higher proportions of "other" taxa for both Illumina MiSeq methods (Fig. 4), especially in the gDNA mock community. Overall, even though DNA sequencing with the Ion Torrent PGM combined with DADA2/Greengenes analysis did not completely provide the expected bacterial composition, this approach resulted in the most accurate representations of the bacteria and their proportions in both gDNA and PCR amplicon mock communities.

DISCUSSION
By comparing DNA sequencing methods, analysis algorithms, and reference databases using dairy relevant bacterial DNA (PCR amplicon and gDNA) mock communities, we found that the DADA2/Greengenes data analysis methods with the Ion Torrent PGM yielded the most accurate interpretations of the 16S rRNA V4 variable region relative to the other methods (Illumina MiSeq, QIIME 1, RDP) tested. This conclusion is notable considering that DADA2 was developed for analysis of Illumina DNA sequence reads (26). Although successfully applied for that purpose (25,33), our findings show that the DADA2 algorithm is compatible with the Ion Torrent reads and error profile. Moreover, our study also offers new and detailed 16S rRNA data comparisons on single-versus multitemplate PCR and single-versus pair-end assembled Illumina reads, which can be broadly informative to benchmark bioinformatics workflows and to the study of bacterial diversity and composition in other microbial habitats besides dairy products.
Application of DADA2 and QIIME 1 (UCLUST) analysis pipelines to the same 16S rRNA gene data showed that DADA2 assigned fewer total and spurious OTUs/ASVs than QIIME 1 even with stringent filtering (32). Because read length was kept consistent FIG 4 Relative abundance of taxa in the 16S rRNA PCR amplicon and gDNA mock communities. Relative abundances of expected taxa are labeled with the corresponding taxonomic level from sequencing results. "Amplicon" represents the 16S rRNA PCR amplicon mock community, and "gDNA" represents the gDNA mock community. The results shown were analyzed following the DADA2 pipeline with the Greengenes database. Each bar represents the mean Ϯ SD from three replicates. Proportions for each community were compared to expected proportions using ANOVA with Bonferroni's multiple-comparison test. P values of Ͻ0.05 were considered to be significantly different from the expected values and are indicated by an asterisk above each bar plot. within each DNA sequencing and data assembly platform, platform-specific differences in OTU/ASV numbers were mainly derived from the core algorithms used for filtering and clustering representative sequences. The DADA2 core algorithm includes errorrate-based denoising, isBimeraDenovo chimera identification, and ASV inference (26). In QIIME 1, the core analysis includes the USEARCH chimera identification and OTU picking strategy (34). Comparison of OTUs and ASVs using the QIIME 1 and DADA2 pipelines, respectively, also showed that the DADA2 analysis pipeline was able to assign ASVs to more specific taxonomic levels (genus) than QIIME 1. This could be the result of the different taxonomy classifiers employed by DADA2 (RDP's naive Bayesian classifier) and QIIME 1 (UCLUST classifier) (35).
OTU and ASV taxonomy assignments were also compared with consideration to 16S rRNA gene reference databases. Results from the different combinations of analysis methods and reference databases showed that the majority of OTUs/ASVs detected were representatives of bacterial taxa included in the 16S rRNA PCR amplicon and gDNA mock communities. Each bacterial species was represented by at least a single OTU/ASV. Additional OTUs/ASVs were largely due to low-abundance sequence variants. DNA sequences of the predominant ASVs/OTUs were 100% identical between gDNA and PCR amplicon mock communities, further supporting the precision of the technique. When QIIME 1 was applied, the RDP_GOLD database (36) yielded lower numbers of total OTUs than found with Greengenes 13.8, independent of whether the Illumina MiSeq or Ion Torrent PGM was used to generate the DNA sequence reads. However, the RDP_GOLD database has not been updated since 2011 (36) and could potentially be missing many bacterial sequences, leading to less differentiation between OTUs. With the DADA2 pipeline, ASVs were inferred prior to taxonomy assignment (26), resulting in the same total ASV numbers for both the RDP 11.5 and Greengenes 13.8 databases. However, assignments of DADA2 ASVs were still influenced by reference databasespecific taxonomic nomenclature and DNA sequences (37), such that Greengenes provided deeper, more accurate taxonomic assignments than those found with RDP.
The Illumina MiSeq and Ion Torrent PGM methods also clearly impacted the outcomes of our mock community analyses. The Illumina MiSeq is well established and known for its low error rate, high-volume read outputs, and low sequencing cost per Gb (38,39). Although Illumina MiSeq reads had higher Phred quality scores, for both single-end and paired-end assembled Illumina MiSeq reads, greater numbers of unexpected taxa and OTUs/ASVs were observed compared to the Ion Torrent PGM. This finding could be the result of differences in library preparation methods, external contamination, index switching (40), and/or substitution errors (41,42) specific to the Illumina MiSeq. To reduce misassigned reads, previous studies have suggested using a dual-index strategy (43) and stringent filtering at the index region (40), as well as sequencing of negative controls for in silico removal of contaminant reads (44). In contrast, the use of the Ion Torrent PGM with our read trimming parameters resulted in the lowest numbers of DADA2 assigned ASVs. At 13 ASVs for both mock communities, this number was only slightly greater than the nine predicted. All 13 ASVs were repeatedly assigned to members of the mock communities, except for one lowabundance ASV when RDP was applied. The four additional ASVs were the result of read errors in the homopolymer regions, a common Ion Torrent error model (20,38,39) that still passed the DADA2 filtering with recommended parameters (https://benjjneb .github.io/dada2/faq.html#can-i-use-dada2-with-my-454-or-ion-torrent-data). This error model could be further reduced by increasing the homopolymer error penalty value. Interestingly, the Ion Torrent PGM reads resulted in the highest numbers of OTUs when QIIME 1 was used to analyze the data. This might have been due to the higher number of erroneous reads that were passed by QIIME 1 filtering, but were identified as sequence chimeras and artifacts by DADA2.
For the PCR amplicon mock community, bacterial diversity analyses based on the DADA2/Greengenes pipeline showed that the results from the Illumina MiSeq were similar to Ion Torrent PGM and the in silico expected values. However, the gDNA mock community relative abundances of certain bacteria in the gDNA mock community were significantly altered compared to expected proportions according to the paired-and single-end Illumina MiSeq methods. This was particularly the case for Streptococcus and Lactococcus. Because Streptococcus and Lactococcus have similar 16S rRNA gene sequences, variation in their relative abundances could be caused by the accumulation of substitution errors, a common error that occurs with the Illumina MiSeq instrument (41,42). No spurious taxa were found in either the PCR amplicon or gDNA mock community in Ion Torrent results analyzed with the DADA2/Greengenes pipeline. In contrast, the gDNA mock community contained 9-fold and 5-fold more "other" spurious taxa compared to the PCR amplicon mock community in Illumina paired-and single-end results, respectively. To this regard, the majority of these taxa (and proportion of reads, Ͼ74%) were assigned to bacterial orders, families, and genera that are highly related to the species included in the mock community (e.g., Clostridiales, Lactobacillus, Oscillospira, and Turicibacter). This, together with the lower numbers of ASVs found for the corresponding PCR amplicon community and single-end results, indicates that errors resulting from paired-end Illumina MiSeq assembly are augmented by combining multitemplate PCR with joining forward and reverse reads. This issue can be mitigated by using single-end reads (as shown by the data here), fewer PCR cycles (33), and increasing the denaturing time (45).
By the use of bacterial DNA standards from nine dairy-relevant bacterial species, we found that DNA sequencing and analysis pipelines contributed significant variations to OTU/ASV distributions and observed bacterial diversities. Moreover, PCR biases and errors from multitemplate DNA amplifications are not entirely filtered with the Illumina MiSeq method. Overall, the Ion Torrent PGM DNA sequencer combined with the DADA2/Greengenes pipeline led to more accurate OTU/ASV assignments and bacterial diversity measurements of the PCR amplicon and gDNA mock communities under our study conditions. The Ion Torrent PGM method is recognized for shorter run times, lower instrument cost, and flexibility in sequencing scale per run by the use of different sequencing chips (38,39). Therefore, this platform could be of particular use to study dairy and other food products with short shelf life times. Moreover, with DADA2 being wrapped in the QIIME 2 platform, we agree with the QIIME 2 developers that new sequencing results should be analyzed using QIIME 2 with a standardized analysis pipeline (e.g., DADA2) instead of QIIME 1 (UCLUST) (46). Further improvements might be reached by refinements to taxonomy classifiers (35), updating reference databases to emphasize bacteria found in different environments, such as dairy foods, and/or testing other reference databases, such as SILVA (37,47,48). Lastly, we recognize that upstream sample processing and DNA extraction protocols can introduce significant biases into assessments of bacterial community composition (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15). Therefore, the data analysis methods applied here should be tested using whole-cell mock communities containing different proportions of bacteria as well as on complex environmental samples. Moreover, to increase reproducibility, consistent methodology and inclusion of negative and positive controls in each run/project are recommended (49). The findings here and the continued development of microbial diversity analysis methods should result in even more reliable comparisons within and between bacterial habitats.

MATERIALS AND METHODS
Bacterial strains and culture conditions. Bacterial strains representing species commonly found in bovine milk were used to construct a mock bacterial community (Table 1). Each bacterial strain was grown in standard laboratory culture medium with negative controls for that species and harvested at early stationary phase by centrifugation at 13,000 ϫ g for 2 min. The laboratory culture media were as follows: Bacillus subtilis, Pseudomonas fluorescens, and Escherichia coli, LB (Lennox broth; Thermo Fisher Scientific); Enterococcus faecalis and Streptococcus agalactiae, brain heart infusion broth (Thermo Fisher Scientific); Staphylococcus aureus, tryptic soy broth (Becton Dickinson); Corynebacterium bovis, tryptic soy broth (Becton Dickinson) with 0.1% Tween 80; Lactococcus lactis, M17 broth (Becton Dickinson) with 0.5% glucose; and Clostridium tyrobutyricum, reinforced clostridial broth (Becton Dickinson). All strains were incubated at 37°C, with the exception of B. subtilis, L. lactis, and P. fluorescens, which were incubated at 30°C. B. subtilis, C. bovis, E. faecalis, E. coli, and P. fluorescens were grown under aeration (250 rpm).
Genomic DNA extraction and PCR amplification. Genomic DNA was extracted using the MagMAX Total nucleic acid isolation kit (Thermo Fisher Scientific, Vilnius, Lithuania) according to the manufacturer's protocol with the repeat bead beating method on a FastPrep-24 instrument (MP Biomedicals LLC).
The DNA concentration was measured with the Qubit 3.0 fluorometer using the Qubit double-stranded DNA (dsDNA) HS assay kit (Life Technologies, Eugene, OR). PCR amplification was performed using Ex Taq DNA polymerase (TaKaRa, Otsu, Japan) and primers F515 and R806 (50) with a random 8-bp barcode on the 5= end of F515 for sample multiplexing (51,52). PCR was initiated at 94°C for 3 min, followed by 35 cycles of 94°C for 45 s, 54°C for 60 s, and 72°C for 30 s, with a final extension step at 72°C for 10 min. Negative controls were run for each barcoded primer. No PCR product for the negative controls was observed on a 1.5% agarose gel. PCR products were pooled and then gel purified with the Wizard SV gel and PCR clean-up system (Promega, Madison, WI).
Preparation of the mock communities. A schematic experimental design for preparing the mock communities is shown in Fig. 1. For the gDNA mock community, 100 ng gDNA isolated from each of the strains was pooled in three separate replicates. The proportion of each bacterial strain in the gDNA mock community was determined by taking into account the genome size and 16S rRNA gene copy number (Table 1). To construct the amplicon mock community, gDNA of the nine bacterial strains was amplified in triplicate by using three different barcoded PCR primers. Amplicon concentrations were measured with the Quant-iT PicoGreen dsDNA assay kit (Life Technologies, Eugene, OR) prior to pooling at equal molar concentrations.
DNA sequencing. For Illumina sequencing, the KAPA HTP library preparation kit (KK8234, Kapa Biosystems, Pittsburgh, PA) was used for the ligation of NEXTflex adapters (Bioo Scientific, Austin, TX) to the 16S rRNA amplicons prior to 250-bp paired-end sequencing (with 7% PhiX control) on an Illumina MiSeq instrument at the University of California, Davis, Genome Center (http://genomecenter.ucdavis .edu/). For Ion Torrent sequencing, non-barcoded Ion A and Ion P1 adapters were ligated to the pooled amplicons, followed by templating, enrichment, and sequencing on the One-Touch 2 and One-Touch ES systems and Ion PGM using the 400 sequencing kit and a 318 v2 chip (Life Technologies, Carlsbad, CA).
16S rRNA gene sequence analysis. An in silico mock community, termed "expected," was created using the 16S V4 amplicon sequences from published genomes and reference genomes for the specific bacterial species (Table 1). In addition, the expected 16S V4 region copy numbers were normalized based on the genome size and 16S rRNA gene copy numbers.
Illumina MiSeq sequencing outputs were trimmed with the fastx_tools (53) to keep the first 245 and 170 bases for the forward and reverse reads, respectively (for quality profiles, see Fig. S6 and S7 in the supplemental material). The Ion Torrent sequence output BAM file was converted to FASTQ format using BEDTools (54), and reads shorter than 200 bp were also removed. The first 280 bases of the Ion Torrent reads were kept for analysis (for quality profiles, see Fig. S6 and S7 in the supplemental material).
The FASTQ files were then analyzed with QIIME version 1.9.1 and DADA2 1.6.0 (26,55). In QIIME 1, Illumina reads from the two orientations (forward and reverse) were analyzed either with or without assembly where the join_paired_ends.py (fastq-join method) (56) script was used with minimum 100-bp overlap and 1% maximum difference between overlapping sequences. Ion Torrent single-end and paired-end assembled Illumina FASTQ files then had the barcode (8 bases) and primer regions (forward primer, 21 bases; reverse primer, 20 bases) removed and were demultiplexed using the split_libraries-_fastq.py script with no barcode error and quality filtered at Q30. Chimeric sequences were identified using USEARCH (34,36) with both the de novo and reference-based methods against the Greengene database version 13.8 (57,58) via the identify_chimeric_seqs.py command with default parameter values. Sequences from both Illumina and Ion Torrent as well as the in silico mock community with expected proportions were merged as one fasta file for operational taxonomic unit (OTU) clustering using the pick_open_reference_otus.py script with recommended parameters (32) and the UCLUST method at 97% similarity thresholds. The Greengenes version 13.8 (57,58) and RDP_GOLD (36) databases were used as references for OTU assignments. Archaea, chloroplasts, and low-abundance (0.005%) OTUs were removed from the OTU tables (32).
In DADA2, for single-end analysis, the truncated Illumina and Ion Torrent FASTQ files after barcode (8 bases) and primer sequence (forward primer, 21 bases; reverse primer, 20 bases) trimming were demultiplexed using split_libraries_fastq.py script with no barcode error and no quality filter (-r 999, -n 999, -q 0, -p 0.0001). Since the single-end reads were already quality trimmed, no additional truncation was performed in DADA2 to be consistent in read length with QIIME 1 analysis. For paired-end analysis, in order to get matched sequence files, raw Illumina reads were demultiplexed in pairs using the idemp tool (59) with no barcode error. Barcode and forward and reverse primer regions were then trimmed with fastx_tools (53). The resulting reads were truncated in DADA2 to keep the first 196 bases of the forward reads and 121 bases of the reverse reads, which were later merged after ASV inference with no error allowed and a 51-bp minimum overlap to be consistent with the QIIME 1 method in resulting read length. For reads from both Ion Torrent and Illumina MiSeq, the error model learning [learnErrors()], dereplication [derepFastq()], and ASV inference [dada()] were performed in R with the DADA2 default parameter, except for added parameters for Ion Torrent [dada(HOMOPOLYMER_GAP_PENALTYϭ-1, BAND-_SIZE ϭ 32)]. Chimeras were identified and removed after sequence clustering via the removeBimeraDenovo() function with the "consensus" method and the isBimeraDenovoTable() function default settings.
Taxonomy was assigned to the resulting amplicon sequence variants (ASVs) using RDP database version 11.5 (60) and Greengenes database version 13.8 with the minimum bootstrap confidence at 80 (57,58). Ion Torrent and Illumina single-end and paired-end assembled reads were merged with the in silico mock community using the phyloseq package in R (61), and singletons and low-abundance (0.005%) ASVs were removed to be consistent with QIIME 1 analysis. Sequences of spurious ASVs were further aligned with sequences in the NCBI nr/nt database using BLASTn (62) with default settings.
Statistics. OTU/ASV counts were rarefied at 5,483 sequences per sample to retain all samples for downstream analyses. Significant differences in the observed mock community composition (␣ diversity and taxonomic distribution) were determined by analysis of variance (ANOVA) with the Bonferroni's multiple-comparison test. A P value of Ͻ0.05 indicates significance. The significance of sample clustering was indicated by permutational multivariate ANOVA using the adonis function from the vegan package in R (63) with a P value of Ͻ0.05 through 9,999 permutations.
Accession number(s). Joined-and single-end DNA sequences after quality filtering and trimming have been deposited in the Qiita database (https://qiita.ucsd.edu) under study ID no. 11351 and in the European Nucleotide Archive (ENA) under accession no. ERP104377.


A meta-analytic review of the association between two common SNPs in miRNAs and lung cancer susceptibility

Background MicroRNAs (miRNAs) are involved in many biological processes, including tumor suppression. Multiple studies have shown an association between the miRNA-196a2 rs11614913 and miRNA-146a rs2910164 polymorphisms and cancer risk. However, the implications of the reported data are debatable and inconclusive. Materials and methods Relevant articles were retrieved from the PubMed, EMBASE, China National Knowledge Infrastructure, and WanFang databases from January 1, 2007, to April 30, 2017. Studies were assessed based on designated inclusion and exclusion criteria, and data were manually extracted from relevant studies by two investigators. Pooled odds ratios (ORs) and 95% confidence intervals (CIs) were calculated to explore the association between two single-nucleotide polymorphisms (SNPs) in miRNAs and lung cancer susceptibility. Results Nine eligible articles were included, consisting of 3,101 cancer cases and 3,234 controls for miRNA-196a2 rs11614913, and 3,483 cases and 3,578 controls for miRNA-146a rs2910164. For studies evaluating miRNA-196a2 rs11614913, significant associations with lung cancer risk were discovered. Overall, the pooled analysis showed that miRNA-196a2 rs11614913 was associated with a decreased cancer risk (CC vs TT: OR = 1.25, 95% CI: 1.09–1.44; CT vs TT: OR = 1.26, 95% CI: 1.03–1.53). For miRNA-146a rs2910164, only the CC genotype was found to be associated with high lung cancer risk (OR = 1.30, 95% CI: 1.13–1.49). Subgroup analyses based on ethnicity, source of control group, and country indicated that there were strong associations between miRNA-146a rs2910164 and cancer risk. Conclusion The results indicated that lung cancer risk was significantly associated with miRNA-196a2 rs11614913 and miRNA-146a rs2910164. These two common SNPs in miRNAs may be potential biomarkers of lung cancer.

Introduction
Lung cancer is the most commonly diagnosed type of cancer and has the highest mortality rate, leading to approximately 1.38 million deaths in the world per year. 1 Studies show that smoking is the primary risk factor for lung cancer. However, although more than 80% of lung cancer patients have smoked, only 15% of smokers will eventually be diagnosed with lung cancer, indicating that other factors such as genetic predisposition, contribute to the development of lung cancer. 2 MicroRNAs (miRNAs) are an abundant category of endogenous small noncoding RNAs of about 18-25 nucleotides in length, which negatively regulate their target mRNAs via posttranscriptional gene silencing. 3 Studies have shown that about 70% of mammalian submit your manuscript | www.dovepress.com Dovepress Dovepress 2420 Xiao et al miRNAs are located in transcription units, mostly in intron regions. 4 miRNAs are present in various forms, and mature miRNAs exert a broad range of biological effects by binding to the mRNA of the target gene. Research suggests that miRNA regulates up to 30% of human genes, 5 and accumulating evidence supports a significant association between miRNAs and the risk and prognosis of lung cancer. 6,7 Takamizawa et al reported that lower expression of the let-7 miRNAs aggravated lung cancer, 8 while Hayashita et al found that upregulation of the miRNA-17-92 may decrease lung cancer risk. 9 Single-nucleotide polymorphisms (SNPs) can change the properties of miRNAs, thereby influencing an individual's susceptibility to cancers. [10][11][12] The miRNA-146a rs2910164 13 and miRNA-196a2 rs11614913 polymorphisms, 14 in particular, have been found to be relevant to lung cancer susceptibility. However, the conclusions of the published studies are often in conflict. 15,16 For example, Volpe et al found that miRNA-196a2 regulated the HOX gene family, which is associated with lung cancer metastasis. miRNA-196a2 rs11614913 polymorphisms regulated the expression of miRNA-196a2 itself and influenced the risk of lung cancer. 15 However, Vinci et al did not find any connection between miRNA-196a2 rs11614913 polymorphisms and the risk of lung cancer, although they did find that miRNA-146a rs2910164 polymorphisms had a role in enhancing the risk of lung cancer. 16 Therefore, in this study, a meta-analysis of previously and newly published data was performed to further explore the correlation between SNPs of miRNA-146a rs2910164 and miRNA-196a2 rs11614913 and lung cancer risk.

Literature retrieval method
We searched the relevant published literature from January 1, 2007, to April 30, 2017, in English or Chinese, using the PubMed, EMBASE, China National Knowledge Infrastructure, and WanFang databases, with the following keywords: "lung cancer", "miRNA-146a", "miRNA-196a2", "polymorphism", "rs2910164", and "rs11614913". References of the retrieved articles were published in the primary literature and had no distinct overlap of population with other studies. The inclusion criteria were as follows: 1) case-control study; 2) case group for the diagnosis of lung cancer patients and control group for non-lung cancer population; and 3) data on the miRNA-196a2 rs11614913 and miRNA-146a rs2910164 polymorphisms. We excluded reports with the same data or overlapping data by the same authors.

Data extraction
Two reviewers extracted data from eligible studies in duplicate with a standard data collection form and reached a consensus on each item. The following information was extracted for each study: first author, publication date, country, ethnicity, source of control, genotype of miRNA-196a2 rs11614913, and miRNA-146a rs2910164 polymorphism.

Statistical analysis
Odds ratios (ORs) and 95% confidence intervals (CIs) were used to measure the degree of correlation between the polymorphism of rs2910164 or rs11614913 and the genetic susceptibility to lung cancer. All analyses used the wild genotype as a control. The combined OR values for CC vs TT, CC vs GG, CT vs TT, and CG vs GG, respectively, were calculated.
The combined OR value was calculated using a z-test with a significance level of P = 0.05. The Q value and I 2 value were used to evaluate the heterogeneity between the results of the study. The calibrated Q-test based on the chi-square test was used to evaluate whether the distribution of effects of each study had significant heterogeneity. With α = 0.10 as the demarcation standard, the heterogeneity of the effect of the study was not statistically significant at P . 0.1, and the study results were merged by the Mantel-Haenszel method. When P # 0.1, a random-effects model (inverse-variance method) was used to combine the results. I 2 was used to describe the proportion of the variation between the studies accounting for total variation. Publication bias was assessed by visual inspection of the Begg's funnel plot. If the funnel plot was symmetrical, it was considered that there was no publication bias in the literature being considered; otherwise, it was considered biased. Our Begg's funnel plot presented a symmetrical and funnel-like distribution of points, indicating that the publication bias risk was low in this analysis. Subgroup analyses were conducted to investigate the association of the subgroups and to identify any heterogeneity.

Study characteristics
We searched the online databases using search criteria related to miRNA-196a2 rs11614913 and miRNA-146a rs2910164 SNPs and the risk of lung cancer. Thirteen published articles were initially collected. Four studies were excluded after duplicate removal, primary screening by reading titles and abstracts, and secondary screening by intensive reading. One of these was written by Tian and its data overlapped. 17 One only reported survival times, 18 one performed a metaanalysis, 19 and another did not report sufficient related data. 20 In total, nine studies met our inclusion criteria. Of these, seven were about miRNA-196a2 rs11614913, with Six studies were about miRNA-146a rs2910164, with 3,483 and 3,578 participants in the case and control groups, respectively (Table 1).

Meta-analysis statistic results of miRNA-196a2 rs11614913
In the CC vs TT group, the heterogeneity test demonstrated the following: χ 2 = 10.03, P = 0.12. The heterogeneity was not significant. The OR was 1.25, with a 95% CI of 1.09-1.44, and P = 0.002 using the fixed-effects model. The results showed that people with the CC genotype had an increased risk of lung cancer, and that the difference was statistically significant. For the CT vs TT group, the heterogeneity test demonstrated the following: χ 2 = 13.07, P = 0.04, I 2 = 54%, with moderate to high heterogeneity. The OR was 1.26, with a 95% CI of 1.03-1.53, and P = 0.02 according to the random-effects model. The results showed that individuals with the CT genotype had an increased risk of lung cancer, and that the difference was statistically significant ( Figure 1). Heterogeneity: χ 2 = 10.03, df = 6 (P = 0.12); I 2 = 40% Test for overall effect: Z = 3.12 (P = 0.002) Yin et al 21 Shang et al 22 Yin et al 23 Vinci et al 16 Hong et al 24 Kim et al 25  Heterogeneity: τ 2 = 0.03; χ 2 = 13.07, df = 6 (P = 0.04); I 2 = 54% Test for overall effect: Z = 2.29 (P = 0.02) Yin et al 21 Shang et al 22 Yin et al 23 Vinci et al 16 Hong et al 24 Kim et al 25 Tian et al 26

Meta-analysis statistic results of miRNA-146a rs2910164
For the CC vs GG group, the heterogeneity test showed that χ 2 = 2.64, P = 0.75, and that the heterogeneity was not significant. The OR was 1.30, with a 95% CI of 1.13-1.49, and P = 0.0002 according to the fixed-effects model, and the difference was statistically significant. The results showed that individuals with the CC genotype had an increased risk of lung cancer, and the difference was statistically significant. For the CG vs GG group, the heterogeneity test showed that χ 2 = 4.05, P = 0.54, and the heterogeneity was not significant. The OR was 1.08, with a 95% CI of 0.96-1.22, and P = 0.18 according to the fixed-effects model (Figure 3). A Begg's funnel plot was constructed; again, there was no observed publication bias in this study (Figure 4).

Results of subgroup analysis by ethnicity
Stratified analyses significantly reduced the heterogeneity of the subgroups. The results of subgroup analysis by ethnicity are shown in Table 2. There was only one study of ethnicity

2423
SNPs in miRNAs and lung cancer susceptibility for Caucasians, and so no heterogeneity test was performed for Caucasian subgroups.
A significant increased risk of lung cancer based on ethnicity was found. For miRNA-196a2 rs11614913, the CC and CT genotypes were associated with increased susceptibility to lung cancer in the Asian population compared with the TT genotype, with a statistically significant difference; however, in the Caucasian population, no polymorphism of this locus was observed to be connected with susceptibility to lung cancer.
For miRNA-146a rs2910164, compared with the GG genotype, the CC genotype was associated with an increased susceptibility to lung cancer in the Asian population, with a statistically significant difference. Based on the current sample size, there was no correlation in Caucasians. The CC genotype was not associated with susceptibility to lung cancer in the Asian population, but CG genotype of Caucasians was associated with an increased risk of lung cancer, with a statistically significant difference.

Results of subgroup analysis by country
The results of subgroup analysis by country are shown in Table 3.
For miRNA-196a2 rs11614913, the CC genotype was not associated with an increased susceptibility to lung cancer compared with the TT genotype in either the Chinese or the Italian subgroup, whereas in the Korean subgroup, an association was found. The CT genotype was not connected with lung cancer risk in any country subgroup.
For miRNA-146a rs2910164, the CC genotype was associated with an increased susceptibility to lung cancer compared with the GG genotype in both the Chinese and Korean subgroups, and the difference was statistically significant. In the Italian subgroup, the CC genotype was not associated with lung cancer susceptibility. However, the CG genotype was associated with an increased susceptibility to lung cancer only in the Italian subgroup.

Results of subgroup analysis by source of control group
The association between SNPs and genetic susceptibility to lung cancer may also depend on the source of the control group. The results of the subgroup analysis according to the source of control group are shown in Table 4.
For miRNA-196a2 rs11614913, the CC genotype was not associated with an elevated susceptibility to lung cancer

2424
Xiao et al compared with TT genotype in either the community or hospital subgroups. The CT genotype was associated with an increased susceptibility to lung cancer in the hospital control subgroup, but not in the community subgroup, and the difference was statistically significant. For miRNA-146a rs2910164, the CC genotype was associated with an increased susceptibility to lung cancer compared with the GG genotype in both the community and hospital control groups, and the difference was statistically significant. However, the CG genotype was not connected with an increased susceptibility to lung cancer in either subgroup.

Discussion
Lung cancer is the most commonly diagnosed cancer, as well as the chief cause of cancer deaths globally. In addition to environmental exposure, SNPs can make an individual more susceptible to lung cancer. Therefore, there is a need to seek out valid biomarkers to predict the risk of the disease.
miRNAs are endogenous small RNAs of about 18-25 nucleotides in length. Previous studies have shown that miRNAs have important roles in key biological processes and are considered a key factor in tumorigenesis. 29 Many researchers have aimed to understand the role of SNPs that occur in precursor and mature miRNAs, and their effects on gene function or expression, cancer susceptibility, and progression of various diseases. 30 Identifying SNPs is important in predicting the risks of individuals and understanding the pathogenesis of cancer. Many studies have explored the common SNPs of miRNAs and their associations with risk of various cancers, including lung cancer. [31][32][33][34][35][36] However, the mechanisms of carcinogenesis are complex, and the investigations on single genetic polymorphisms and studies with small sample sizes that have been performed so far have not provided a means to reliably predict the risk of cancer. In addition, the conclusions of previous studies have sometimes been inconsistent and conflicted with each other. 16,22,37,38 To address this, we carried out a meta-analysis by gathering all of the studies that reported a connection between the miRNA-196a2 rs11614913 and miRNA-146a rs2910164 SNPs and the risk of lung cancer.
For the miRNA-196a2 rs11614913 locus, we found that the risk of lung cancer in the CC and CT genotype individuals was higher than in the TT genotype individuals. For the miRNA-146a rs2910164 locus, individuals with the CC genotype exhibited a higher risk of lung cancer than individuals

2425
SNPs in miRNAs and lung cancer susceptibility with the GG genotype. Thus, the miRNA-196a2 rs11614913 and miRNA-146a rs2910164 loci can serve as potential candidate sites for lung cancer susceptibility screening.
To identify any differences between individuals of Caucasian and Asian ethnicities, we carried out a subgroup analysis. This showed that different ethnicity, to some extent, was linked to differences in the association between genetic polymorphisms and lung cancer risk. On this basis, a further subgroup analysis by country was performed -this was, to our knowledge, the first attempt to explore whether the susceptibility to lung cancer varies across different countries. We considered that the dietary and health-related habits of different countries might influence the risk of lung cancer.
As there was only one study of Caucasian individuals and the participants in this study were Italian, the Italian subgroup was identical to the Caucasian subgroup. Our results indicated that the association between the miRNA-196a2 rs11614913 polymorphism and the risk of lung cancer only could be found in Asian populations. Our initial results indicated that the predictive value of miRNA-196a2 rs11614913 polymorphism was valid in Asian people. However, in the subgroup analysis of correlations between polymorphisms and countries, this relationship was found only in Korean individuals. Differences owing to ethnicity have also been found in other types of cancer. 39 Wang et al 39 found that the miRNA-196a2 rs11614913 CC genotype was significantly associated with breast cancer risk in Chinese and Indian individuals, but not in Caucasians. For the miRNA-146a rs2910164 locus, only the CC genotype was related to an increased susceptibility of lung cancer compared with the GG genotype in the Asian population, whereas the CG genotype of the Caucasian population was associated with an increased susceptibility of lung cancer compared with the GG genotype. After subgroup analysis based on the country, the CC genotype was associated with an increased susceptibility to lung cancer compared with GG genotype in both the Chinese and Korean subgroups. As there was only one report involving Caucasian population, 16 and the sample size of this study was relatively small, further research is needed before any firm conclusions can be drawn regarding this population.
According to the subgroup analysis of the control source, the miRNA-196a2 rs11614913 CT genotype was associated with an increased susceptibility compared with the TT genotype in the hospital subgroup. By contrast, the miRNA-146a rs2910164 CC genotype was associated with a higher risk of lung cancer compared with the GG genotype in both community and hospital subgroups.
There were some advantages of the comprehensive analysis that we have performed here, such as the use of the funnel plot to evaluate publication bias, which showed good funnel plot symmetry. In addition, all of the nine articles in our study were of high quality. However, some limitations still need to be addressed. First, and most importantly, the results of the meta-analysis were limited by the available literature, such as the number of documents, quality, and sample size, which may affect its validity. For example, there were too few Caucasian studies available, and the heterogeneity among different studies cannot be ignored. Furthermore, we did not explore the association of SNP-SNP interactions. Such gene-gene interactions in cancers have been reported in the literature. [40][41][42][43] In addition, lung cancer is a result of gene-environment interactions. In this study, only the influence of genetic factors was taken into account, and environmental factors, such as smoking and air pollution, were not explored. Future work examining these gene-gene and gene-environment interactions and their association with lung cancer is needed.
In summary, our meta-analysis suggested that the miRNA-196a2 rs11614913 CC and CT genotypes, and the miRNA-146a rs2910164 CC genotype, may increase the risk of lung cancer, especially in Asian populations. Owing to the low number of studies included and the relatively small sample size of the single Caucasian study, the negative results in Caucasian individuals need to be confirmed by further study. Overall, it will be necessary to confirm all these findings using well-designed and larger sample size studies. In addition, the interaction of genetic and environmental factors should be considered.

OncoTargets and Therapy
Publish your work in this journal Submit your manuscript here: http://www.dovepress.com/oncotargets-and-therapy-journal OncoTargets and Therapy is an international, peer-reviewed, open access journal focusing on the pathological basis of all cancers, potential targets for therapy and treatment protocols employed to improve the management of cancer patients. The journal also focuses on the impact of management programs and new therapeutic agents and protocols on patient perspectives such as quality of life, adherence and satisfaction. The manuscript management system is completely online and includes a very quick and fair peer-review system, which is all easy to use. Visit http://www.dovepress.com/testimonials.php to read real quotes from published authors.