Stool sampling and DNA isolation kits affect DNA quality and bacterial composition following 16S rRNA gene sequencing using MiSeq Illumina platform

Many studies correlate changes in human gut microbiome with the onset of various diseases, mostly by 16S rRNA gene sequencing. Setting up the optimal sampling and DNA isolation procedures is crucial for robustness and reproducibility of the results. We performed a systematic comparison of several sampling and DNA isolation kits, quantified their effect on bacterial gDNA quality and the bacterial composition estimates at all taxonomic levels. Sixteen volunteers tested three sampling kits. All samples were consequently processed by two DNA isolation kits. We found that the choice of both stool sampling and DNA isolation kits have an effect on bacterial composition with respect to Gram-positivity, however the isolation kit had a stronger effect than the sampling kit. The proportion of bacteria affected by isolation and sampling kits was larger at higher taxa levels compared to lower taxa levels. The PowerLyzer PowerSoil DNA Isolation Kit outperformed the QIAamp DNA Stool Mini Kit mainly due to better lysis of Gram-positive bacteria while keeping the values of all the other assessed parameters within a reasonable range. The presented effects need to be taken into account when comparing results across multiple studies or computing ratios between Gram-positive and Gram-negative bacteria.


Results
We analyzed stool samples from sixteen volunteers. Each volunteer collected the samples from the same stool sample using three different sampling kits (SK): a stool container (SK1); a flocked swab (SK2) and a cotton swab (SK3). The DNA was extracted using two isolation kits PowerLyzer PowerSoil DNA Isolation Kit (PS) and QIAamp DNA Stool Mini Kit (QS) (see Methods), totaling 96 samples for the analysis.
Evaluation of user acceptance of the sampling kits. The participants were asked to select the best and the worst kit based on their ease of manipulation including the time spent using it. All 16 volunteers selected the stool container as the easiest to use and 13 out of 16 (81.25%) volunteers indicated the flocked swab as the worst sampling kit. We believe that the manipulation with cotton and flocked swabs is uncomfortable due to the small size and the necessity to insert the swab stick back into the tube without touching the tube wall. On the contrary, the stool container is easy to manipulate even for people with reduced motoric skills. In addition, the flocked swab is designed for sampling of liquid samples and the solid stool samples do not adhere on its synthetic fibers.
The effect of sampling and DNA isolation kits on the bacterial gDNA quality. DNA yield, purity and integrity. Significantly higher DNA yields were obtained with the QS isolation kit, regardless of the sampling kit used (q < 0.01) (Fig. 1, Supplementary Table S1). The median values of the A260/A280 ratio (the measure of purity of DNA) were well within the expected range (1.8-2) and did not differ significantly between the DNA isolation kits or between the sampling kits ( Fig. 1, Supplementary Table S1).
The DNA integrity was determined using the GQN measure (on a scale from 1 to 10; low GQN indicates strongly degraded gDNA sample) and the proportion of short fragments (≤1500 bp; the larger the proportion the more degraded gDNA). We observed interaction effects of isolation and sampling kit for both DNA integrity measures. We found significantly lower proportion of short fragments when using the PS isolation kit (Fig. 1, Supplementary Table S1) and this difference was much larger when the stool container was used for sampling. There was no difference in GQN measure between the isolation kits when cotton or flocked swabs were used. However, for stool container samples, the QS kit provided much lower GQN values compared to the PS kit. These results point to worse DNA integrity for the QS isolation kit compared to the PS isolation kit mostly when stool container is used for sampling.
Presence of PCR inhibitors. The presence of PCR inhibitors in the samples decreases the sensitivity of the PCR reaction and even can lead to the impossibility of amplification of the selected region of 16S rRNA. It is usually measured by median efficiency values estimated from inhibition plots. Ideally, the efficiency should be 100%, meaning the template doubles in each cycle. Usually, the efficiency within 90-110% range is considered acceptable, where lower efficiency is caused by non-optimal reagent concentration or lower enzyme quality, while higher efficiency values are caused by the presence of PCR inhibitors. In our data, the efficiency values ranged from 96.7% to 114.0% (Fig. 1, Supplementary Table S2). In each of the isolation/sampling kit combinations, there were minimum two samples which exceeded the efficiency of 110%. The efficiency values of all isolation/sampling kit combinations, except for stool container samples after DNA isolation with the QS kit, were significantly increased compared to control samples without PCR inhibitors (efficiency med = 94.7%). No difference in efficiency values was observed between the isolation kits. The samples from stool containers (regardless of the isolation kit used) contained less PCR inhibitors in comparison to all other sampling/DNA isolation kit combinations (significantly lower efficiency, Supplementary Table S2). We hypothesize that this sampling kit effect is due to the sample dilution step prior to the DNA isolation step.
The effect of sampling and DNA isolation kits on bacterial diversity and composition. Bacterial diversity. In total, 96 stool samples were sequenced. The number of reads after quality filtering and removal of chimeras ranged from 27680 to 67809, with median of 46192. We assessed the bacterial diversity using the number of observed OTUs and the Chao 1 diversity metric ( Fig. 1, Supplementary Table S1). Overall, both diversity measures were independent of the DNA yield in all sampling/DNA isolation kit combinations.
While there was no difference in Chao 1 measure between the isolation kits, the number of observed OTUs was significantly increased after isolation with the PS kit, but only for cotton swab samples (q-value = 0.029). When comparing diversity measures between the sampling kits within each isolation kit separately, the stool container resulted in significantly higher number of observed OTUs in both DNA isolation kits (Fig. 1, Supplementary  Table S1). In addition, we observed significantly higher number of OTUs in flocked swab samples compared to cotton swab samples after DNA isolation with the PS kit (q-value = 0.04) and significantly lower number of OTUs in flocked swab samples compared to cotton swab samples after DNA isolation with the QS kit (q-value = 0.09). For the Chao 1 diversity metric, significant differences were found in stool container samples compared to flocked swab samples in both PS and QS isolation kits (q = 0.04 and q = 0.09, respectively).
Bacterial composition. We identified 12,948 OTUs belonging to 13 bacterial phyla.
In order to quantify the effect of the sampling and isolation kits on bacterial composition, we performed mixed linear regression on each taxon that passed the filtering criteria (maximum abundance across all samples ≥1%) at all the seven taxonomical levels (phylum, class, order, family, genus, species, OTUs) separately. Interestingly, the proportion of taxa significantly affected by isolation or sampling kit differed between taxonomical levels (Fig. 2). The choice of sampling or DNA isolation kit affected 100% of taxa at phylum, class and order levels, and had decreasing trend from family to OTU level. The effects of sampling and isolation kits on the ten most abundant taxa at different taxa levels are summarized in Table 1 (see Supplementary Tables S3-S8 for complete results), the composition of significantly affected families is shown in Fig. 3. Overall, the choice of the isolation kit affected the abundance of more taxa than the choice of the sampling kit. In most of the cases where the taxa was affected by both factors, the p-values associated with the effect of the isolation kit were smaller than those of the sampling kit, indicating a more significant contribution of isolation kit to the overall model.
We hypothesized that the observed effect of the isolation kit was a result of different efficiency of the kit-specific bacterial cell walls lysis procedure. In this case, one of the kits would be more successful in isolating  Table 2 shows the numbers of significantly affected G+ taxa in all taxonomic levels and statistical pairwise comparison of their proportion after both isolation methods and all sampling methods. We found significantly higher proportions of G+ bacteria after the isolation using the PS kit at all the taxon levels. (96.4% to 100%, Table 2), compared to the QS isolation kit (G+ proportion varying from 0 to 44%). Similar observations were made for the effect of the sampling kit ( Table 2), but this trend was not significant on any of the taxa levels except for the comparison of cotton swab (SK2) and stool container (SK1) on the genus level. We hypothesize that these differences are attributed to the dilution of the samples during the preprocessing steps specific to the stool container (see Methods for more details), resulting in lower sample density thus increasing the efficiency of the bead beating procedure. No difference in proportion of Gram-positive bacteria was found between flocked and cotton swabs. Figure 4 shows estimated effect sizes pairwise between the sampling kits on the genus level. Figure 5 visualizes bacteria with significant changes in abundance between isolation or sampling kits, with nodes colored according to Gram-positivity, where we can observe association of Gram-positive bacteria with the PS isolation kit.

Discussion
The gut microbiome seems to be crucial factor influencing human health and to date, a number of different diseases were correlated with microbiome dysbiosis. Understanding the true role of microbiome and fully comprehending its variability will require many cohort studies and, most probably, comparison of their results in large-scale meta-analyses. As with any other scientific domain, the incoherent methodological approaches constitute an important obstacle for such comparisons 44 . In an attempt to elucidate some of the factors determining the success of such studies, we focused on the effects of sampling and DNA extraction methods on a number of relevant variables from DNA integrity to final bacterial composition at different taxa levels. For this purpose, we selected sampling and DNA isolation kits that are the most common and accessible and hence are probably the most relevant for majority of cohort studies.
Our group of sixteen healthy volunteers used three different sampling kits -stool container, flocked swabs and cotton swabs. Without exception, the stool container was indicated as the most acceptable by the volunteers. Moreover, stool in the container can be easily diluted, homogenized and aliquoted for different analyses. Unfortunately, the stool container is inconvenient for sampling diarrhea or baby stool. Importantly, as we discuss below, the pre-processing specific to stool container samples influences both DNA quality and bacterial composition and these effects seem to interact with the DNA isolation kit.
For measuring the effect of different DNA extraction procedures, we used PowerLyzer PowerSoil DNA Isolation Kit (PS) and QIAamp DNA Stool Mini Kit (QS).
While the PS kit cell-wall lysis procedure is based on combination of bead-beating step and enzymatic lysis, the standard protocol of the QS kit comprises only enzymatic lysis. Considering the fact that the beat-beating step leads to higher DNA yield and higher number of observed OTUs from difficult-to-lyse bacteria, we added the bead-beating step also into the QS protocol, as commonly recommended 8,30,34,35,39 .
DNA isolation by the QS kit resulted in significantly higher DNA yields compared to the PS kit (regardless of the sampling kit). Similar results were observed in other studies 30,32 . In agreement with previous studies 30,35,37 , we found no significant correlation between DNA yield and alpha diversity. The proportion of taxa significantly affected by sampling or isolation kit at different taxonomical levels. Proportion of the tested taxa significantly affected by the sampling kit only (green), by the isolation kit only (dark yellow) and by both sampling and isolation kit (brick red). Grey indicates taxa not affected by sampling or isolation kit. The significance level was chosen at FDR < 10%, only taxa that met the selection criteria (maximum abundance >1%) were tested.   www.nature.com/scientificreports www.nature.com/scientificreports/ On the other hand, the PS kit produced DNA of better integrity, even though in the PS protocol we applied more rigorous mechanical lysis (or higher speed of bead beating), which, according to the literature, should result in more degraded DNA 48 . We hypothesize that the observed differences might be caused by another factor, such as the type of the beads (0.1 mm glass in PS vs 0.1 mm zirconia in QS), the buffer composition, or the incubation temperature. Overall, for preparation of the shotgun libraries or sequencing using third generation of sequencers, we consider DNA integrity to be more important factor than the DNA yield, which favors PS kit over the QS kit.
To properly homogenize the samples from the stool container, we included a preprocessing procedure comprising five times dilution. This naturally resulted in lower yields of isolated DNA, but after adjustment for this dilution we obtained higher final DNA concentrations compared to undiluted stool samples from flocked and cotton swabs. It seems that the dilution step also affected the DNA integrity. Compared to the undiluted samples from flocked and cotton swabs, stool container samples resulted in less degraded DNA after isolation using the PS Figure 3. Distributions of relative abundances of significantly affected taxa at family level. Four graphs represent families divided according to third quartile of their abundance. Only taxa that passed the filtering criteria (maximum abundance >1%), significantly affected by isolation or sampling kit are shown. The colored squares below the graph indicate whether the family was affected significantly by the sampling kit only, the isolation kit only or both.  www.nature.com/scientificreports www.nature.com/scientificreports/ kit and, in contrast, in more degraded DNA after isolation using the QS kit. Interestingly, two other independent studies, where different isolation kits were used, showed either a negative 34 or a positive 48 effect of sample dilution on the DNA integrity. This, together with our results leads us to conclude, that the effect of dilution step on DNA integrity is dependent on the isolation kit.
PCR inhibitors persisted in the DNA of the samples after isolation with both kits. Presence of PCR inhibitors could complicate the use of conventional molecular methods for the detection of low abundance or rare www.nature.com/scientificreports www.nature.com/scientificreports/ pathogenic microorganisms 49,50 . The dilution of stool container samples prior to processing has led to significantly lower proportion of PCR inhibitors, hence for some applications, this approach might be preferred.
Both DNA extraction kits isolated preferentially bacterial DNA, independently on the sampling kit used and the amount of human DNA was negligible. From practical point of view, there is no superiority of any of the DNA isolation vs sampling kit combinations with respect to amount of residual human DNA. Some of the studies, however, use these kits to estimate the concentration of human DNA in stool samples as an indicator of inflammation that might predict onset of certain bowel diseases [51][52][53][54][55] . From this perspective, based on our results, we do not consider these kits eligible for human DNA quantification.
As for the alpha diversity, we observed increased number of OTUs after DNA isolation with the PS kit in all sampling kits, but the difference was significant only for cotton swab samples. We observed significant differences in number of OTUs between all sampling kits combinations, with the stool container resulting in the highest number of OTUs. We attribute the observed differences to higher effectivity of bead beating process in the less dense samples (the dilution preprocessing step used for the stool container). This is in contrast with the results of Santiago et al. 34 , who report no changes in alpha diversity after sample dilution. In that study, however, a different isolation kit was used, so the results are not directly comparable.
The final bacterial composition was more affected by the choice of the DNA isolation kit than by the choice of the sampling kit. The preference of the PS isolation kit for Gram-positive bacteria was confirmed by statistical testing on all taxa levels and we believe that it is a result of more effective lysis of the Gram-positive cell wall bacteria when using the PS kit, despite the additional bead-beating step we introduced into the QS protocol. This is in agreement with previously published results 8,26 . It has to be taken into account, that Gram staining not always corresponds with the cell wall structure (e.g. Pseudobutyrivibrio 56 or Deinococcus 57 , which is for many bacteria unknown. The efficiency of the lysis procedure can be as well influenced by atypical composition of the cell wall, presence of S-layer or capsules. The bacterial cell wall type also plays a role in the sampling effect: in our study it was associated with the dilution preprocessing step of the stool container, although less significantly. There is a common belief that the effect of the individual is the most influential on the final bacterial composition 8,32 . Indeed, many metagenomic studies are reporting differences between groups of interest at the OTU level, where the effect of isolation and sampling is less important, as we showed in this study. However, some hypotheses are connecting particular disease with higher or lower bacterial abundance at the phylum or family level. An example is the commonly used Bacteroidetes/Firmicutes ratio [58][59][60][61][62][63][64] . Our results show, that this ratio is very dependent on both the selected DNA isolation method and sampling kit (dilution step). In our study, the PS kit and the dilution step (stool container) led to significantly higher proportion of e.g. Firmicutes (G+) and Actinobacteria (G+) and significantly lower proportion of Proteobacteria (G−) and Bacteroidetes (G−).
Another example of the cell wall structure effect is the Gram-positive genus Blautia. Blautia is a common and highly prevalent bacteria in the gastrointestinal tract, which is connected with healthy gut, since it is an effective short-chain fatty acid producer 65,66 . Lower abundance of Blautia in the gut is associated with many diseases [66][67][68][69][70][71][72][73] . In our study, Blautia was bacteria the most significantly affected by DNA isolation (across all the taxonomic levels). Similar observations were also described as the effect of isolation in other studies 26,34 .
The sampling kit (dilution effect) influenced most significantly the abundance of genus Sutterella, bacteria correlated with many diseases such as celiac diseases 67 , Down syndrome 74 , autism 75 or irritable bowel syndrome 76 . Clearly, the dilution step represents an important batch effect, which raises a question, whether it is related only to the artificial dilution, or this effect could also be observed in diarrheic samples. The effect of stool consistency was described previously as an important factor 12,77,78 influencing the bacterial composition, but this effect was www.nature.com/scientificreports www.nature.com/scientificreports/ not connected with effect of higher water content (dilution), rather with the transit time. As previously recommended 77 , we also suggest to control for the stool consistency as a potential confounding factor to avoid the effect of sample water content in this kind of studies, especially if one of the illness symptoms is diarrhea.
Despite the fact that the significance of the sampling and isolation dependent batch effects is repeatedly reported, no systematic study of these effects was performed yet on samples from larger numbers of individuals. Efforts for standardization of laboratory practices in metagenomics have been made in large international projects such as Metagenomic Research Group (MGRG), Genomic Standard Consortium (GSC), The Microbiome Quality Control Project (MBQC) and International Human Microbiome Standards (IHMS). IHMS recommends a procedure for fecal sample DNA extraction, based on study of Costea et al., where 21 extraction protocols were compared, including protocols similar to ours -protocol 3 (with PowerLyzer PowerSoil DNA Isolation Kit) and 11 (with QIAamp DNA Stool Mini Kit and bead beating step) 39 . They selected the protocol with QIAamp DNA Stool Mini Kit as the best choice for its accuracy and reproducibility. In contrast to our results, both protocol 3 and 11, provide good lysis of Gram-positive bacteria, but protocol 3 was excluded for insufficient DNA quality. The main difference between the studies is that the Costea study was based on the results of whole metagenomics sequencing and only compared bacterial composition annotated at the species level.
All these above mentioned studies and our results confirm that meta-analytical studies are extremely challenging due to the many sources of batch effects that need to be accounted for. Incorporation of a standardized mock community to the sequencing workflow, followed by normalization of the results to these reference values could be solution in future. The increased cost per run and slightly more complex library preparation is a small price to pay for robustness, consistency and comparability of results.

Conclusions
We performed systematic study of effects of DNA isolation and sampling kit on DNA quality and bacterial composition based on sequencing of gene for 16S rRNA on a the largest number of individuals to day (96 samples from 16 individuals).
We found significant effect of both DNA isolation and sampling kits on DNA purity, DNA integrity, alpha diversity and bacterial composition. Overall, the DNA isolation effect was stronger than that of the sampling kit. Interestingly the proportion of taxa affected by isolation or sampling was decreasing with decreasing taxonomical level.
We confirmed previously reported effect of DNA isolation kit on bacterial composition due to bacterial cell wall structure, namely the better efficacy of The PowerLyzer PowerSoil DNA Isolation Kit in lysis of Gram-positive bacteria. In addition, we report that the dilution pre-processing step of the stool container samples favored Gram-positive bacteria, although mostly at the genus level.
Both the choice of isolation and sampling kits significantly affected the Firmicutes to Bacteroidetes ratio. We conclude that the choice of DNA isolation and sampling kit (dilution step, and by extension the stool consistency) is an important batch effect that has to be taken into account mainly when comparing results between studies.

Methods
Sample collection. Stool samples were collected from a group of 16 volunteers. The subjects were 23-65 years old with an average age of 40.9 and none of them suffered from diarrhea during sample collection. Stool samples were collected at home. Volunteers received three stool sampling kits: sampling kit 1 (SK1) comprising 1x stool container (FL Medical, Italy); sampling kit 2 (SK2) comprising 2x flocked swabs (Copan, Italy) and sampling kit 3 (SK3) comprising 2x cotton swabs (SceneSafe, Great Britain). Sampling kits also contained disposable gloves and hand and surface disinfectant wipes for more convenient sampling. Each volunteer was instructed to collect all the samples from the same stool and from the same spot. Stool samples were then stored in a freezer at −20 °C overnight to freeze completely and the next day were transported on ice buckets to the laboratory, where they were stored at −20 °C prior to processing. Each group of samples was processed at the same time and by the same person. Participants filled out a brief questionnaire about satisfaction with individual sampling kits after stool sample collection. The study design is summarized in Fig. 6.
This study was carried out in accordance with the recommendations of the ELSPAC Steering Committee of Masaryk University with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocols were approved by the ELSPAC Steering Committee of Masaryk University. DNA extraction. Stool in the stool container (SK1) was diluted 5x with molecular grade water and homogenized by vortexing with Zirconia beads 2.3 mm (BioSpec, USA) to receive identical aliquots. This step is not necessary for the swabs, since each swab serves as an aliquot itself. Stool suspension (250 µl) was used for DNA extractions. Flocked swabs (SK2) and cotton swabs (SK3) were transferred into 2 ml tubes to be prepared for subsequent DNA extraction. DNA extractions were performed using a PowerLyzer PowerSoil DNA Isolation Kit (Mo Bio, USA) (PS) and QIAamp DNA Stool Mini Kit (Qiagen, USA) (QS) according to the manufacturer's instructions. Deviations from PS protocol: • 750 µl of Bead Solution and 60 µl of C1 Solution were added to swab samples (SK2 and SK3) after defrosting.
Samples were thoroughly vortexed and centrifuged for 4 min at 36,220 RCF. The swabs were then removed. Next, the samples were homogenized using the FastPrep-24 (MP Biomedicals, USA) 45 s 6.5 m/s. Deviations from QS protocol: • A homogenization step with 0.1 mm zirconia beads (BioSpec, USA) was added to the protocol after the third step (i.e. after the suspension was heated for 5 min at 95 °C).

Presence of PCR inhibitors after different DNA extractions.
The presence of inhibitors was tested by qPCR. A primer pair specific for the conservative regions of 16S rRNA gene ( Proportion of human DNA to bacterial DNA after different DNA extractions. The ratio of human and bacterial DNA in samples was tested by qPCR. Bacterial DNA was assessed using a primer pair specific for the conservative regions of 16S rRNA gene and human DNA using a primer pair specific for protein kinase (Table 3). qPCRwas performed on the TOptical Thermocycler (Analytik Jena -Biometra, Ireland) with KAPA SYBR FAST qPCR Kit (Kapa Biosystems, USA). Cycling conditions are described in Table 3. Melting temperature Figure 6. Study design. Flowchart summarizing the study design and methods used.  www.nature.com/scientificreports www.nature.com/scientificreports/ was determined after PCR to verify the correctness of each PCR product. The amount of human DNA to bacterial DNA was calculated as 2 ΔCt . Ct value of 40 was used for all samples under the limit of detection.
PCR amplification and Illumina library preparation. Extracted DNA was used as a template in amplicon PCR to target the hypervariable V3 and V4 regions of the bacterial 16S rRNA gene. The 16S metagenomics library was prepared according to the Illumina 16S Metagenomic sequencing Library Preparation protocol with some deviations described below (for workflow diagram see Supplementary Fig. S1). Each PCR was performed in triplicate, with the primer pair consisting of Illumina overhang nucleotide sequences, an inner tag and gene-specific sequences 79 . The Illumina overhang served to ligate the Illumina index and adapter. Each inner tag, i.e. a unique sequence of 7-9 bp, was designed to differentiate samples into groups. Primer sequences and PCR cycling conditions are summarized in Table 3. After PCR amplification, triplicates were pooled and the amplified PCR products were determined by gel electrophoresis. PCR clean-up was performed with Agencourt AMPure XP beads (Beckman Coulter Genomics, USA). Samples with different inner tags were equimolarly pooled based on fluorometrically measured concentration using Qubit dsDNA HS Assay Kit (Invitrogen, USA) and microplate reader Synergy Mx (BioTek, USA). Pools were used as a template for a second PCR with Nextera XT indexes (Illumina, USA). Differently indexed samples were quantified using the KAPA Library Quantification Complete Kit (Kapa Biosystems, USA) and equimolarly pooled according to the measured concentration. The prepared library was checked with a 2100 Bioanalyzer Instrument (Agilent Technologies, USA) and concentration was measured with qPCR shortly before sequencing. The library was diluted to a final concentration of 8 pM and 20% of PhiX DNA (Illumina, USA) was added. Sequencing was performed with the Miseq reagent kit V3 using a MiSeq. 2000 instrument according to the manufacturer's instructions (Illumina, USA).

Data analysis.
Forward and reverse pair-end reads, that fulfilled the condition of both quality and length filtering, were merged using the fastq-join method within the join_pair_ends.py command in QIIME 1.9.1 80 . Data were demultiplexed and barcodes and primers were trimmed using package Biostrings 81 in R 3.3.2 82 . Operational taxonomic units (OTUs) were constructed by binding sequences into clusters of greater than 97% sequence similarity using QIIME. In the next step, chimeras were detected on the set of representative sequences of each OTU with UCHIME in USEARCH v6.1.544 83 . These chimera OTUs were subsequently excluded from the analysis. Taxonomy was assigned to each OTU based on SILVA 123 reference database 84 . The observed species metric and the Chao1 index were used to estimate alpha diversity for each sample in QIIME. Beta diversity was computed in QIIME using both weighted and unweighted UniFrac metrics 85 . All statistical analysis was performed in R 3.3.2 82 .
The data were treated as compositional (proportions of total read count in each sample, non-rarefied) and prior to all statistical analyses were transformed using centered log-ratio transformation 86 . The analyses were performed on each of the seven taxonomy levels (Phylum, Class, Order, Family, Genus, Species and OTUs) separately and the resulting p-values were adjusted for multiple hypothesis testing using Benjamini-Hochberg procedure. Results were considered significant at FDR = 10%. The adjusted p-values are referred to as q-values.
To estimate the effects of isolation and sampling kits on bacterial composition while accounting for repeated measurements (effect of individual), we applied linear mixed model with sampling and izolation kits as fixed effects and individual as random effect (intercept). Log-likelihood test was performed to detect significance of each of the fixed effects -each time we compared the full model to the model without the fixed effect of interest.
A non-parametric Wilcoxon paired test, was used for comparison of effect of isolation kits on DNA quality. We used Spearman's rank order correlation coefficient to discover the strength of the link between the number of observed species and DNA concentration.
Bipartite networks were used to visualize the influence of different kits on detection of Gram-positive and Gram-negative bacteria. These networks were reconstructed according to Sedlar et al. 87 using R 3.3.2 and visualized in Gephi 0.9.2 88,89 . Communities within networks were extracted using modularity optimization criterion 88 .

Data Availability
Sequencing data were uploaded to the European Nucleotide Archive under accession number PRJEB24411. Read counts per sample at different taxa levels and sample information table are available in Supplementary Files S9-S11.