Miseq reporter v2

Sequencing of the GLA gene by NGS was performed in collaboration with KM Biologics Co., Ltd. The procedures were described in our previous report [18 ]. Briefly, the 13.3-kbp region including the GLA gene was amplified by long-range PCR (Supplementary Fig. S2). Library preparation and sequencing were performed using a Nextera XT Kit (Illumina, San Diego, CA, USA) and MiSeq sequencer (Illumina). After sequencing runs, the data were aligned to the human reference genome sequence (NC_000023.10) using MiSeq Reporter software (Illumina). Sequence data analysis, mapping, and variant calling were streamlined using MiSeq Reporter v2 (Illumina). Visualization of sequencing reads was performed with IGV_2.3.10 (Broad Institute). Variants detected in the GLA gene by NGS were resequenced by the Sanger method.

The raw sequence reads were demultiplexed, and the reads with quality scores of less than 20 were removed by the MiSeq reporter v2.5 (Illumina). The subsequent procedure of sequence assembly, quality check, and assignment of taxonomy was performed using USEARCH v.11.0.667^{31 (link)}. Reads with more than 1.0 expected errors were removed. The joined reads less than 200 bp were also removed. The default USEARCH settings were used for the identification of unique sequences. Using the UNOISE algorithm, chimeric reads were removed, and the sequence reads were clustered into zero-radius operational taxonomic units (ZOTUs). Each ZOTU was taxonomically assigned against the ITS2 database^{32 (link),33 (link)} using the SINTAX algorithm with an 80% confidence level^{34 (link)}. The ZOTUs with the same taxonomy identification were combined for further taxonomic analysis.

The sequences were demultiplexed and trimmed with a read quality score above 20 by the MiSeq Reporter v2.5 (Illumina, Inc.). The sequence data obtained was processed following the Miseq SOP in Mothur^{46 (link)}. Sequences with any ambiguous bases, sequences with more than 8 homopolymers and sequences with lengths less than 200 bp were removed using the screen. Seqs command in Mothur. Putative chimeric sequences were detected and removed via the Chimera Uchime algorithm contained within Mothur in de novo mode. Rare sequences (less than 10 reads) were removed to avoid the risk of including spurious reads generated by sequencing errors^{47 (link)}.High quality sequences were assigned to OTUs at ≥99% similarity level.
Taxonomic classification of each OTU was obtained by classifying alignments against SILVA Release 115 databases^{48 (link)} using the classify command in Mothur at 80% cutoff with 1000 iterations. The Miseq sequence data used in this study are deposited in the MG-RAST server (Meyer et al., 2008) under project ID 17384 (http://metagenomics.anl.gov/linkin.cgi?project=17384).

All sequence data were quality filtered (QF) prior to taxonomic assignment using taxonomy-dependent workflows. Metabarcoding reads recovered by paired-end sequencing were first stitched together using the Illumina MiSeq analysis software (MiSeq Reporter v 2.5) under the default settings. Sequences were then assigned to samples based on their unique index combinations and trimmed in Geneious Pro v 4.8.4⁸³ . In order to eliminate low quality sequences, only those with 100% identity matches to Illumina adaptors, index barcodes, and template specific oligonucleotides were kept for downstream analyses. Sequences were further processed in USEARCH v 9.2^{73 (link)}, which was used to trim ambiguous bases, remove sequences with average error rates >1% and those that were <200 bp in length, dereplicate each sample down to unique sequences, abundance filter the unique sequences (minimum of two identical reads)^{13 (link)}, and remove chimeras. To enable comparison among replicates, the 18 S sequences in each replicate were sub-sampled to 20,000 reads prior to dereplication. Due to lower amplification success of the ITS2 marker, replicates within sites were merged and subsampled to 20,000 and 4,000 reads prior to dereplication for sediment and seawater samples, respectively.

Downstream analysis of sequences was performed with QIIME 22019.1 [51 (link)]. Raw sequence data were quality filtered and denoised with DADA2 [52 (link)]. The number of filtered and non-chimeric sequences were described in Additional file 3: Table S2. The sequences were aligned and positions that were highly variable were masked using the mafft program [53 (link)], which was used for phylogenetic diversity analyses such as unweighted and weighted UniFrac using FastTree [54 (link)]. Alpha-diversity metrics (number of observed OTUs and Shannon’s diversity index) and beta diversity metrics (unweighted and weighted UniFrac) were estimated using q2-diversity after samples were rarefied to 16,000 sequences per sample. Taxonomy classification was determined using the MiSeq Reporter v2.3 based on an Illumina-proprietary classification algorithm and an Illumina-curated version of the Greengenes taxonomy database. Microbial functions of uterine microbiota were predicted based on 16S rRNA genes using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt)-1.1.4 [55 (link)]. The taxonomic and metabolic profiles were analyzed in linear discriminant analysis effect size (LEfSe) [56 (link)] and the Statistical Analysis of Metagenomic Profiles (STAMP) v2.1.3 [57 (link)] in order to identify features that were statistically significant between groups.