To assess the reproducibility and accuracy of DNA extraction and sequencing steps, mock controls and technical replicates were included on each plate (Figure 1). A ZymoBIOMICS™ Microbial Community Standard (catalog number ZR D6300), consisting of eight known bacterial species including both gram-positive and gram-negative bacteria, was included on each plate as a DNA extraction control. Each plate included a ZymoBIOMICS™ Microbial Community DNA Standard (catalog number ZR D6305), which contains already extracted DNA for eight known bacterial species. This served as a sequencing control. Seventeen within-run and 8 between-run technical replicate pairs were also included across the sequencing plates. These included samples that were collected across several different time points during infancy, specifically 0–5 weeks, 3 months, 6 months, 9 months and 12 months.

Twenty six pairs of stool and rectal swab samples were collected at the age of 0–5 weeks (baseline samples), for infants born between 37 and 42 weeks gestational age. Swabs were stored in Primestore solution (PrimeStore® Molecular Transport medium). All stool and swab samples were transferred to a −80°C freezer for long-term storage.

Stool samples were thawed, and half of a pea-sized scoop was collected from the side/centre of the sample. This was placed in a tube with 750 µL of lysis buffer. For rectal swab samples, 400 µL of sample in Primestore was placed in a tube with 400 µL of lysis buffer. These then underwent off-board lysis, using the QT Qiagen bead beater, prior to DNA extraction.

Manual DNA extraction of the stool and rectal swab samples and mock extraction controls was carried out using the Quick-DNA™ Fecal/Soil Microbe Microprep Kit (ZymoBIOMICS catalog number D6012). Prior to carrying out polymerase chain reaction (PCR), Qubit™ (Thermo Fisher Scientific, Waltham, MA, United States) was done to check the starting DNA concentrations. xGen™ 16S Amplicon Panel v2 kits (Integrated DNA Technologies, Coralville, IA, United States) were used in library preparation for sequencing. Kits included primer pairs for amplification of all nine hypervariable regions of the 16S rRNA gene. Additionally, the primers for these kits have dual indices to allow greater numbers of samples to be run together in a single flow cell. Moreover, the xGen kits include Normalase™ which could be used to enzymatically normalise library sizes prior to sequencing. Finally, qPCR was performed following the Normalase step to quantify the final library size prior to sequencing.

Negative controls, including Primestore, Milli-Q water, elution buffer and Tris EDTA/Nuclease free water, were added to each plate (Figure 1) together with prepared libraries from the stool and rectal swab samples. Mock controls and technical replicates, as described above, were also included on each plate. Sequencing was conducted across two sequencing runs on seven plates. The combined 16S library per run was subjected to paired-end sequencing on the Illumina® MiSeq™ platform, employing the MiSeq Reagent v3 kit with 600 cycles (Illumina, San Diego, CA, United States).

Ethics approval for this research was provided by the Human Research Ethics Committees at the University of Cape Town (801/2016 and 557/2020) and at Stellenbosch University (M16/10/041). Following sequencing, preprocessing steps were carried out for quality control and to prepare the data for statistical analysis (Figure 2). Raw sequencing data was run through FastQC (Andrews, 2010) to assess the quality of the reads. Following this quality control step, forward and reverse reads for each sample were processed using the SNAPP-py3 pipeline (Chai, 2021). Sequencing data from run 1 and run 2 were processed separately. Of the four main output files from the pipeline, the lineage table and an adapted taxonomy table were used for further analysis.

The remaining processing and analysis were carried out in R version 4.2.1 (R Core Team, 2022). This stage of processing began with creating a phyloseq object (McMurdie and Holmes, 2015; McMurdie and Holmes, 2013). Using the decontam package (Davis et al., 2018), decontamination was carried out separately for each plate, using plate-specific negative controls. A combined frequency and prevalence approach was used, selecting a threshold of 0.1 for the prevalence component. Phyloseq objects from runs 1 and 2 were then combined into a single phyloseq object for further downstream processing and analysis.

A normalisation step to account for different library sizes was implemented by determining the median library size and normalising each sample accordingly (Balle et al., 2020; The Jackson Laboratory, 2019). Finally, subsetting into various phyloseq objects was done to prepare for downstream statistical analysis. We ultimately had separate phyloseq objects containing the mock extraction controls, mock sequencing controls, within-run repeats, between-run repeats and baseline pairs of stool and rectal swab samples. Excel spreadsheets containing this data, as well as the corresponding code for importing the files into R as phyloseq objects, are provided in the Supplementary Datasheets S1, S3, respectively.

Batch effect correction was done using MMUPHin (Ma, 2022; Ma et al., 2022). This data was compared to data in which no batch effect correction was carried out, to determine the necessity of accounting for batch effects.

Relative abundances for controls, replicates and stool/rectal swab sample pairs were visualised using QIIME2 software (Bolyen et al., 2019). All microbiome analysis was done at a species level in R. Genus-level and phylum-level analyses were additionally included in select steps to provide additional insights.

Performance measures were calculated as outlined by Özkurt et al. (2022) using data for sequences classified to a species-level. Precision and sensitivity were calculated based on the number of correctly identified species expected to be in the mock control [true positives (TP)], the number of expected species that were not detected [false negatives (FN)] and the number of non-expected species classified as being in the control [false positive (FP)]. F-scores could then be calculated based on these values. The calculations used were as follows: Precision = TP / TP + FP Sensitivity = TP / TP + FN

F-score = 2*precision*sensitivity/(precision + sensitivity) (Özkurt et al., 2022).

The percentage relative abundances of each of the eight expected species were determined for mock cell (extraction) and mock DNA (sequencing) controls on each sequencing plate. For each control, the total percentage of sequences that were correctly classified as an expected species, was calculated. Furthermore, observed relative abundances were compared to the theoretical abundances provided by ZymoBIOMICS for each species in the mock controls by calculating Observed/Expected (O/E) ratios (Maki et al., 2023).

Functions from the phyloseq package in R (McMurdie and Holmes, 2015; McMurdie and Holmes, 2013), specifically the estimate_distance and distance function, were used to calculate alpha and beta diversity measures for replicates and stool/rectal swab samples. The alpha diversity measures included are observed richness (Fisher et al., 1943), Shannon’s index (Shannon, 1948) and Simpson’s index (Simpson, 1949). Bray Curtis (Bray and Curtis, 1957) and Jaccard’s (Ludwig and Reynolds, 1988) distances were the beta diversity measures included in our analysis. Plot_richness and plot_ordination functions were used to plot alpha and beta diversity measures, respectively. Paired Wilcoxon tests were used to compare alpha diversity measures between pairs of within-run replicates and to identify differences between pairs of between-run replicates. In order to determine whether beta diversity measures were reproducible between technical replicate pairs, distance-based intraclass correlation coefficients (dICCs) were calculated separately for within-run and between-run replicates (Chen and Zhang, 2022). Paired Wilcoxon tests and dICCs were similarly used to compare pairs of stool and rectal swab samples collected from the same infants.

All eight expected bacterial species were detected in four of our seven mock extraction controls (Table 1). Bacillus subtilis was not detected to the species level in two controls (Table 2; Figure 3A), however classification to a genus level (Bacillus) was achieved (Supplementary Table S1). Listeria monocytogenes was not detected even at a genus level in the control from run 2, plate 1 (Supplementary Table S1). For the three controls in which we were unable to detect all eight species, the total percentage of sequencing data correctly classified as expected mock species was consequently lower (Table 2). Sensitivity scores were over 0.88 for all controls. Precision scores were lower – particularly for the control on run 1, plate 1, which had a score of 0.32. This was driven by a high number of false positive results in this control. A median F-score of 0.84 (range of 0.47–1.00) was obtained for the mock extraction controls (Table 1).

The percentage abundances of species in these controls did not accurately follow the order of theoretical abundances for some species (Table 2; Figure 4A). In particular, the relative abundances of L. monocytogenes were well below the theoretical abundances suggested by ZymoBIOMICS at both a species and genus level. This is emphasised by the low median Observed/Expected (O/E) ratio of 0.32 at a species level. Similarly, Enterococcus faecalis and Staphylococcus aureus had O/E ratios well below the value of 1. The relative abundances of Escherichia coli and Salmonella enterica were greater than expected, with a range of O/E ratios all lying well above the value of 1.

Among the mock sequencing controls, the eight anticipated bacterial species were detected in five of the seven controls (Table 3; Figure 3B). The overall percentages of sequences correctly classified as expected species were slightly lower for these controls compared to the mock extraction controls (Table 4). For the run 2 plate 1 control the prevalence of S. enterica was particularly low, and this was not resolved at the genus level (Supplementary Table S2). B. subtilis again was not detected at a species level for two of these controls (Table 3). Precision scores for the mock sequencing controls ranged from 0.50 and up, while sensitivity was greater than 0.88 for all controls. F-scores had a median of 0.80 (range of 0.67–0.94).

The relative abundances of species in the mock sequencing controls more closely matched the expected abundances than was observed for the mock extraction controls (Table 4; Figure 4B), as seen by the O/E ratios being closer to 1. In these controls the abundance of Limosilactobacillus fermentum was well below the theoretical threshold expected. The O/E ratios for this species and S. aureus were consistently less than 1. Whereas E. faecalis, E. coli and L. monocytogenes had O/E ratio ranges above 1.

Although the inclusion of a batch effect correction step was trialled, substantial differences in the relative abundances of the various species across the controls were observed compared to when no batch effect correction was done (Supplementary Figure S1). Similarly, stricter decontamination thresholds led to poorer reproducibility in mock controls.

Similar patterns in the relative abundance of species could be seen when comparing pairs of within-run and between-run repeats (Supplementary Figures S2, S3). When comparing alpha diversity of within-run repeats using paired Wilcoxon tests, we found no evidence to indicate differences at a species level in the Observed richness [95% confidence interval (CI) (−3.50, 5.50); p = 0.587], Shannon’s index [95% CI (−0.09, 0.12); p = 0.782] or Simpson’s index [95% CI (−0.02, 0.01); p = 0.487] between these pairs of technical replicates (Figure 5). Furthermore, when comparing the beta diversity distance matrices of these within-run technical replicate pairs (Figure 6), we observed a good level of reproducibility for Bray Curtis as seen by a distance-based intraclass correlation coefficient (dICC) value of 0.940 and Jaccard’s distance showed good, albeit lower, reproducibility with a dICC of 0.762. Similarly, we found no differences between these technical replicates when looking at their alpha and beta diversity measures at a genus level.

A comparison of alpha and beta diversity measures for between-run technical replicates similarly found no clear differences between these pairs at a species level (Figure 7). Paired Wilcoxon tests comparing Observed species [95% CI (−2.00, 13.00); p = 0.362], Shannon’s index [95% CI (−0.05, 0.49); p = 0.195] and Simpson’s index [95% CI (−0.004, 0.11); p = 0.148] did not motivate for the existence of differences in alpha diversity measures. Moreover, dICC results showed good reproducibility for Bray Curtis (dICC = 0.899) and moderate reproducibility for Jaccard’s distance (dICC = 0.597) (Figure 8).

To assess whether stool and rectal swab samples may be combined for analysis, we compared 26 sample pairs collected at 0–5 weeks after birth. The relative abundances of species did not display similar patterns across pairs of samples (Supplementary Figure S4) and alpha diversity was found to differ when running paired Wilcoxon tests. Stool and swab samples from the same participants differed at a species level in terms of Observed species (p < 0.001) and Shannon’s index (p = 0.027), with no differences in Simpson’s index (p = 0.394) found (Figure 9). A comparison of beta diversity measures at a species level found that while Bray Curtis measures were moderately reliable between the pairs of stool and swab samples (dICC = 0.684), there was poor reliability when comparing their Jaccard’s distances (dICC = 0.310) (Figure 10).

We were able to obtain relatively good species-level resolution in terms of sensitivity using a nearly complete 16S rRNA sequence, which the SNAPP-py3 pipeline developers refer to as a ‘consensus’ sequence. The eight species expected to be found in the ZymoBIOMICS mock controls were not consistently detected to a species level in all mock controls across the seven sequencing plates, suggesting that there is still room for improvement in terms of reproducibility. The main species which were not detected in all controls are B. subtilis and L. monocytogenes. These species are gram-positive, containing a strong wall of peptidoglycan which can be a challenge to lyse, as has been presented in previous literature (Claassen-Weitz et al., 2020). ZymoBIOMICS have intentionally developed these mock controls to include both gram-negative and gram-positive species to enable researchers to identify inconsistencies and to optimise their lysis protocols (ZymoBIOMICS, 2024). Moreover, research indicates that the primers used in library preparation may have a greater affinity for some species compared to others (Klindworth et al., 2013), which may also contribute to false negatives in some controls. Our results indicate that the methods used are capable of detecting all species, as seen in several of our mock controls, however future optimisation of the protocol will be required.

A previous study comparing the use of different variable regions in Illumina Miseq sequencing, found that using 1-2 variable regions could at best identify 16 of 20 (80%) mock control species correctly (Fouhy et al., 2016). In five mock controls, we detected 7 of 8 mock species (88%), yet for the remaining nine controls we detected all the expected species (100%). Thus, our results would indicate that using all nine variable regions of the 16S rRNA gene improves accuracy compared to methods that look at only a couple of variable regions. Other short-read multiple variable region methods have been able to identify all taxa within mock controls containing a greater number of species (Fuks et al., 2018; Schriefer et al., 2018), however in these studies there were also mock controls in which not all species were identified. Schriefer et al. (2018) suggested that the depth to which sequencing was done may play a role, however an increase of sequencing depth by 10-fold ultimately had no impact on their results. Therefore, the approach we present requires further optimisation to ensure that 100% sensitivity is consistently achieved and to assess whether this approach can perform at a similar level to other multiple variable region tools when using more complex mock controls.

When focusing only on sequences that were classified to a species level, sensitivity was high in both the mock extraction and mock sequencing controls. However, our results indicate that high sensitivity comes at the expense of obtaining poorer precision in some controls. The overall performance for mock extraction and mock sequencing controls were median F-scores of 0.84 and 0.80 respectively. These results are comparable with other research that obtained F-scores greater than 0.80 using a bioinformatics tool for analysing sequencing data, for which the authors described their results as being indicative of good performance (Özkurt et al., 2022). Poor precision scores in our case were often driven by false positive taxa that were present in very low relative abundances. Thus, removing rare taxa could yield even better results. Schriefer et al. (2018) and Fuks et al. (2018), also reported false positives at low abundances for their multiple variable region approaches.

As our measure of precision is based on the presence or absence of expected taxa, it is quite strict in how it penalises false positives. It is not weighted according to the abundances at which these occur. Therefore, it is important for us to look at the precision scores in conjunction with our abundance-based measures, to get a complete picture of the performance of the methodology used in this paper. The R1P1 mock extraction control displayed poor precision due to having a particularly high number of false positives. However, as these occurred at such low abundances, false positives had little influence on our abundance-related analysis. Rather, we found that the inability to detect all expected species (i.e., false negative results) in our controls was the main factor responsible for differences in observed abundances compared to expected abundances.

The overall percentages of sequences that were classified as expected taxa and to a species level, ranged between 79.11% and 99.98% for mock extraction controls, and between 81.80% and 99.58% for mock sequencing controls. In a study by Szoboszlay et al. (2023) the authors found that between 58.9%–68.9% of reads obtained from Illumina sequencing of the V4 region could be correctly classified to a species level. This improved if rare taxa were excluded. Nanopore sequencing of the entire 16S rRNA gene could yield classification of over 81% of reads to a species level (Szoboszlay et al., 2023). Therefore, although we did not remove rare taxa, our protocol could compete with a full-length sequencing approach and perform better than using a single variable region for short-read sequencing. Ultimately, our results indicate that the RDP classifier is able to accurately identify the expected taxa in our controls.

The relative abundances of species in our mock extraction controls sometimes diverged from the theoretical abundances outlined by the suppliers (ZymoBIOMICS, 2022), indicating that there may be some bias introduced during the extraction process. O/E ratios provided a quantitative means of assessing the bias introduced in each control. In particular, they assisted with identifying the species for which observed abundances differed most from theoretical abundances and showed that there was greater bias in the mock extraction controls.

There are various stages in 16S rRNA sequencing where prejudice can arise, favouring certain taxa or altering the relative composition of a sample (Nearing et al., 2021). Studies have found that DNA extraction kits and methods can differ in terms of their ability to extract DNA from gram negative and gram positive bacteria, indicating that the results of a study can be influenced depending on which kit and extraction method are used (Videnska et al., 2019; Yuan et al., 2012). The mock sequencing controls (for which already-extracted DNA was obtained from ZymoBIOMICS) more closely resembled the expected abundances. L. fermentum differed substantially from the theoretical abundance in the sequencing controls – unlike the extraction controls for which this particular species had followed the expected abundance more closely. This implies that, for the most part, there is minimal bias introduced at the sequencing step. However, there may be some bias specifically in sequencing L. fermentum. Our findings complement other studies which have shown that the extraction and amplification of DNA is particularly prone to prejudicing results, while there is less bias introduced at the sequencing stage of microbial research (Brooks et al., 2015; Lee et al., 2012). Similar to our findings, previous research using a short-read multiple variable region approach, also found that the observed abundance of taxa did not always accurately match the expected abundances (Fuks et al., 2018). The authors also suggested that bias may be introduced during the amplification stage of library preparation.

Our findings suggest that using this multivariate analysis approach might provide a means to improve the ability of short-read 16S rRNA sequencing studies to study microbial samples at a species level. This approach would benefit from additional ways to minimise bias – particularly in the DNA extraction process.

We considered including a batch effect correction step in the processing of our data. The goal was to reduce bias introduced due to samples being sequenced on different plates and in different runs. However, when visualising the relative abundances of mock sequencing controls, we found that batch effect correction led to far less consistent results across the controls. Davis et al. (2018) reported that the decontam package corrects for batch effects. As a result, it may be redundant to carry out both decontamination and batch effect correction. Moreover, as our protocol included a plate-wise decontamination step, batch effects would have already been taken into account here. Thus, based on our findings, we considered an additional batch effect correction step excessive and problematic – leading us to exclude this step. Similarly, we found that using a very stringent threshold for decontamination also had a negative impact on our analysis. As such, we would suggest that care needs to be taken when doing additional processing of data, to ensure that bias is not introduced to the data by overcorrecting for contaminants and batch effects. Researchers should take care in deciding whether batch effect correction is necessary if their pipeline includes a decontamination step, which removes contaminants according to batches.

Observed richness provides an indication of the number of different taxa within samples, while Shannon’s and Simpson’s indices additionally account for how equally represented these taxa are within samples (evenness) (Kers and Saccenti, 2021). Based on our results we have no reason to reject the null hypothesis that there is no difference in alpha diversity measures between within-run technical replicate pairs.

Jaccard’s distance quantifies how diversity varies between samples based on whether or not taxa are present, while Bray Curtis factors in the abundance of taxa (Kers and Saccenti, 2021; Schroeder and Jenkins, 2018). More confidence is usually placed in Bray Curtis compared to Jaccard’s distance (Schroeder and Jenkins, 2018). The similarity between these beta diversity measures for replicate sample pairs was assessed in terms of dICC. dICC is a distance measure which has been established specifically for microbiome data, building on the concept of intraclass correlation coefficients, to look at the similarity between replicate samples (Chen and Zhang, 2022). Higher ICC values are indicative of strong similarity between measures (Koo and Li, 2016), in our case diversity measures for microbial sample pairs. An ICC of 0.5, however, shows moderate reproducibility and a value of <0.5 suggests poor similarity (Koo and Li, 2016). Thus, our results in which we get dICC values of 0.940 and 0.762 for Bray Curtis and Jaccard respectively, indicate that we get good beta diversity reproducibility when running samples in duplicate on the same plate. Likewise, the diversity of between-run technical replicates was similar with dICC values of 0.899 and 0.597. This suggests that we do not have significant batch effects across sequencing plates of the same run, or across different runs, when following the protocols set out in this study for sequencing and processing short-read multivariate 16S rRNA data. It should be noted that we only had one pair of replicates across runs and this was grouped with the between-plate replicates in what we referred to as our “between-run replicates.” Collectively our findings indicate that there is no batch effect introduced. The sequencing and processing of infant stool and rectal swab samples is consistent, providing confidence in the analysis of the microbiome datasets generated using this protocol.

The implementation of a new kit and pipeline for analysing both stool and rectal swab samples in this study warranted investigation. Stool and rectal swab samples collected from the same infants and at the same time point were found to differ for several diversity measures. Contrary to our findings, studies have suggested that rectal swab samples can be used in place of stool samples and provide a reliable representative measure (Radhakrishnan et al., 2023; Reyman et al., 2019; Bassis et al., 2017). However, a couple of studies exploring diversity and functional roles of the microbiome have previously reported differences between stool and swab samples (Short et al., 2021; Sun et al., 2021). Short et al. (2021) specifically found that beta diversity differed between these sample types, while Sun et al. (2021) found that rectal swab samples varied from stool in terms of both alpha and beta diversity measures. Moreover, there have been indications that differences in the relative abundances of specific taxa may occur between different sample types (Jones et al., 2018). This would be particularly important to consider when comparing the prevalence of individual taxa between groups, such as when using differential abundance analysis.

Furthermore, a study by Bokulich et al. (2019) found that the length of time between the collection and processing of samples is an important factor to take into consideration when using rectal swab samples in sequencing studies. Although samples in their study were ultimately frozen, swabs were posted to the laboratory, while stool samples were immediately placed on ice until they could be collected and taken to the laboratory. Rectal swab samples received and processed after 48 h following sample collection did not accurately match stool samples – with a higher proportion of Enterobacteriaceae being favoured (Bokulich et al., 2019). The stool and rectal swab samples in our study were collected together and stored at −80°C until DNA extraction and sequencing were done, and therefore we do not face the same challenges as Bokulich et al. (2019) in terms of the length of time for which samples remained at room temperature. However, it may be important to consider factors such as the time for which samples are stored in freezers, and the time between DNA extraction and sequencing in future work involving both stool and swab samples.

The fact that swabs were stored in Primestore, while stool was not, would suggest a substantial difference which requires consideration when using a combination of stool and swab samples in carrying out analysis. Several studies have found that bacterial compositions differ across different regions of entire stool samples (Zreloff et al., 2023; Huson et al., 2017; Gorzelak et al., 2015). Spectroscopic analysis has shown that patterns of metabolites can also vary substantially across different regions of stool (Liang et al., 2020; Gratton et al., 2016). Therefore, using only part of the sample does not provide a good characterisation of the entire sample. These studies highlight the importance of homogenising the entire stool sample in order to carry out sequencing (Zreloff et al., 2023; Liang et al., 2020; Gratton et al., 2016). As sequencing of rectal swab samples is not limited to only a component of the collected sample, this sample would be homogenous. This may explain the differences observed between stool and rectal swab sample pairs.

The work presented included only one mock extraction and one mock sequencing control on each sequencing plate. Future work could include duplicates of controls on each plate to better assess reproducibility of the controls. It would also expand the options available for statistically comparing these controls.

A limitation of this study is the inconsistent detection of all eight species in each mock control. We were able to detect all eight species in several controls, suggesting that the methods used for library preparation, sequencing and processing of the data are all capable of achieving consistent detection. Given that gram positive species particularly are not always detected (Claassen-Weitz et al., 2020), and these are more challenging to lyse, future work to optimise the lysis protocol should be done – for example, assessing whether carrying out mechanical lysis for longer enables the consistent detection of all species. A study involving sequencing of the entire genome, has previously shown that sequencing depth may play a crucial role in whether or not all expected taxa are detected (Pereira-Marques et al., 2019). Sequencing depth is a variable that differed for each control and may explain the inconsistencies in our results. Sequencing was done based on the Illumina guidelines which recommend a library depth of at least 100,000 reads (ILLUMINA. 16S Metagenomic Sequencing Library Preparation, 2024). However, as we are using a pipeline which combines multiple variable region reads into consensus sequences, and as other short-read multiple variable region 16S rRNA studies have utilised greater sequencing depths (Fuks et al., 2018; Schriefer et al., 2018), future work could seek to optimise the sequencing depth required for this protocol in order to better identify all species in mock controls.

In this study we evaluated a single pipeline and protocol for analysing data from multiple variable region 16S rRNA sequencing, however in future this protocol would benefit from comparisons to other tools and pipelines. Moreover, different DNA extraction kits and amplicon sequencing kits could be compared to identify the best kits for carrying out multiple variable region analysis. In this analysis we have worked with compositional data, looking at the relative abundances of taxa. However, in future, it would be useful to explore absolute abundances for samples sequenced using this protocol. Moreover, we could explore whether other databases – such as databases specific to the human gut microbiome – might enable more effective and accurate classification of consensus sequences from the SNAPP-py3 pipeline.

Researchers intending to utilise this sequencing kit and pipeline should ideally strive (where possible) to use data from samples obtained through the same collection approach. In cases where a combination of sample collection approaches is used, analysis should be done to assess whether collection approaches influence microbial composition. If they do, sample collection approach should be controlled for in any analysis in which this data is used.

The goal of this paper was to assess the use of xGen kits and the SNAPP-py3 pipeline, to set the stage for an observational clinical study. Our assessment of this methodology is limited by the fact that we did not have access to a wider variety of datasets sequenced using the xGen kits. Future research should explore the use of this methodology in a greater selection of clinical settings and in the study of environmental microbiome samples. This might result in the approach being of greater relevance and interest to a broader audience.

The results of our 16S rRNA multiple variable region analysis, using short-read Illumina sequencing data from all 9 variable regions, show that there is promise for using this protocol for species-level analysis. We have identified areas for improvement and future work should assess whether this approach is comparable to other bioinformatics pipelines and tools that have been established for analysing multiple variable region short-read data. As we found differences in diversity between stool and swab sample pairs, future analysis of data which has been generated using this protocol should take this into account. Furthermore, our findings indicate that when using new sequencing kits and protocols to study both stool and swab samples, it would be advisable to do analysis to compare pairs of stool and rectal swabs from the same individual to confirm whether these yield comparable results.