Study population and sample collection

To study whether rectal swabs are a good proxy for faecal samples in infant gut microbiota studies, we used a subset of 131 paired faecal and rectal swab samples from 116 children participating in the Dutch randomised controlled ZEBRA study. The ZEBRA study aims to evaluate the effects of antibiotic treatment indicated for (suspected) neonatal sepsis in the first week of life on the developing gut microbiota. Written informed consent was obtained from both parents. Ethical approval was granted by the national ethics committee in the Netherlands, METC Noord-Holland (committee on research involving human subjects, M014-024, NTR5119). The study was conducted in accordance with the European Statements for Good Clinical Practice.

Rectal swabs were collected using FaecalSwab™ kits (Copan Diagnostics, CA, USA) by trained physicians or research personnel before the start of antibiotic treatment (timepoint 1) and 24–48 hours after cessation of antibiotic therapy (timepoint 2). Faecal samples were obtained at the same timepoints, usually directly after the rectal swab, this being a stimulatory trigger for defaecation, and stored in sterile faecal containers by a nurse during hospital stay, or by the parents if the participant was already discharged at the later timepoint. All material was directly stored at −20 °C before being transferred (<2 weeks) to a −80 °C freezer until further laboratory processing. We only analysed paired samples that were obtained within 24 hours of one another, and in the case of timepoint 1 were also both obtained strictly before the start of antibiotic treatment.

DNA isolation and sequencing

Bacterial DNA was isolated from faecal samples as previously described18. We used approximately 100 μl of faeces, 300 μl of lysis buffer, 500 μl zirconium beads and 500 μl of phenol, and performed an extra phenol/chloroform step. Samples collected on day 1 were presumed to have low bacterial abundance. Therefore, further adaptations were applied as described previously19, with the additional changes of using 150 μl instead of 100 μl of faeces (or 100 μl of material in the case of rectal swabs) and implementing an extra step with wash buffer 1. DNA blanks and a positive control consisting of a mix of up to three random faecal samples were used for quality control. The amount of bacterial DNA was determined by quantitative polymerase chain reaction (qPCR) as previously described19.

After amplifying the V4 hypervariable region of the 16S rRNA, quantification of the amount of amplified DNA per sample was executed with the dsDNA 910 Reagent Kit on the Fragment Analyzer (Advanced Analytical, IA, USA). Samples yielding insufficient DNA after amplification, defined as <0.5 ng/μl, were repeated with a higher concentration of template DNA. A mock control and three PCR blanks were included in each PCR plate. 16S rRNA sequencing was performed on the Illumina MiSeq platform (Illumina, Eindhoven, the Netherlands).

Bioinformatic processing

The samples and their sequences described in this manuscript are part of a larger dataset existing of 2176 samples and controls, and together were processed using our in-house bioinformatics pipeline20. In short, we applied an adaptive, window-based trimming algorithm (Sickle, version 1.33) to filter out low quality reads, maintaining a Phred score threshold of 30 and a length threshold of 150 nucleotides21. Error correction was performed with BayesHammer (SPAdes genome assembler toolkit, version 3.5.0)22. Each set of paired-end sequence reads was assembled using PANDAseq (version 2.10) and demultiplexed (QIIME, version 1.9.1)23,24. Singleton and chimeric reads (UCHIME) were removed. OTU picking was conducted with VSEARCH abundance-based greedy clustering with a 97% identity threshold25. OTUs were annotated using the Naïve Bayesian RDP classifier (version 2.2) and the SILVA reference database26,27. This resulted in an OTU-table containing 18,951 taxa in total. We created an abundance-filtered dataset selecting OTUs present at a confident level of detection (0.1% relative abundance) in at least two samples28, hereafter referred to as our raw OTU-table. The raw OTU-table consisted of in total 730 taxa (0.49% sequences excluded with filtering). Next, we used both the prevalence and frequency methods of the decontam package29 to exclude possible contaminants, discarding 35 taxa, and thus retaining 695 taxa in total. The subset of paired samples studied here contained only 372 of these taxa.

Statistical analyses

All analyses were performed in R version 3.4.330 within RStudio version 1.1.38331 and figures were made using packages ggplot232 and ggpubr33. The alpha diversity of the two sampling methods was compared using the observed species richness and Shannon diversity. When rarefying to a sequencing depth of 25,000 reads after filtering and decontamination (lowest quartile), or even 15,000 reads, the differences found using the raw data remained, so from here the raw (unrarefied) data was used for the comparisons and correlations in alpha diversity. Group differences were tested for using Wilcoxon tests. The correlation in alpha diversity between paired samples was calculated with Pearson.

The effect of sampling method on composition was analysed univariately with PERMANOVA-tests with 1999 permutations (adonis function; vegan package34) for all samples, and also stratified per timepoint, to prevent confounding by repeated measures. To visualise differences in composition we generated non-metric multidimensional scaling plots (nMDS; vegan package34). Ordinations were based on the Bray-Curtis (BC) dissimilarity matrix of relative abundance data with parameter trymax 10,000. To test whether paired faecal samples and rectal swabs (within child comparison) were more similar in microbiota composition than unpaired samples (between children comparison), we calculated the BC similarities (1 – BC dissimilarity) between all samples and compared the level of similarity between the paired and unpaired samples using Wilcoxon. For the paired samples, we also tested the correlation between composition and the difference in collection time and reads between the two sampling methods.

We performed a temporal, multivariable PERMANOVA-test (adonis2 function, vegan package34, 1999 permutations) to test whether sampling method contributes to the variation in overall gut microbiota community composition and how this relates to other known drivers of community composition. First, we tested covariates known to be associated with gut microbiota composition (age, delivery mode, feeding type) univariately with PERMANOVA-tests as described above. Only covariates that showed a significant association in the univariate analysis (age, delivery mode), were included in a multivariable model along with sampling method, whilst setting the strata parameter to individual.

Finally, we evaluated whether the sampling methods correlated on taxonomical level by calculating the Pearson correlation coefficient for all individual taxa based on their relative abundance. We also calculated the delta mean relative abundance between the two sampling methods to show the differences found.

P-values or, where applicable, adjusted p-values calculated using the Benjamini-Hochberg method35, <0.05 were deemed significant.