In this multicenter prospective study conducted in two neonatal intensive care units (NICUs) in Madrid, Spain, we analyzed the association between nasopharyngeal and gut microbiota profiles, viral infections, and the development of recurrent wheezing in preterm infants. By integrating epidemiological, clinical, and microbiological data, this study aims to identify microbial and viral predictors of recurrent wheezing and to explore potential microbiota-modulating strategies to mitigate respiratory complications in this high-risk population.35,38
A prospective multicentre study (2019-2022) at two Madrid NICUs (H. Severo Ochoa and H. La Paz), investigated respiratory viral infections and microbiota profiles in preterm infants born at <32 weeks' gestational age and <8 days old at the time of inclusion. Infants were excluded if they had major malformations, chromosomal disorders, immunodeficiencies, or if their parents refused participation. Ethical approval was obtained (ref. PI4676) and written informed consent was signed by the parents or legal guardians.
The main clinical outcome was wheezing, defined as ≥2 episodes of physician-diagnosed wheezing requiring medical attention during the first year of life. Additional outcomes included respiratory admissions, anti-asthmatic treatment prescription, and atopy (defined as dermatitis and/or food allergy).
All samples were collected during the initial hospitalization for prematurity. Two simultaneous NPAs samples were collected during the first week of life (specifically between 3 and 7 days of life). One of them was kept refrigerated at 4 °C, until its transport to the National Centre for Microbiology (ISCIII) for virological study, and the other one, was used for the microbiota study and stored immediately at -80°C until use. Additionally, weekly samples of NPA were collected for virological study.
Faecal samples were also collected between days 3 and 7 of life, following a standardized protocol and using sterile techniques. Samples were immediately stored at -80°C until analyses. Sample collection was performed by trained neonatal nurses.
RNA and DNA from 200 µL aliquots of NPA were extracted using the QIAamp MinElute Virus Spin Kit in an automated extractor (QIAcube, Qiagen, Valencia, Spain) from the NPA samples stored for the National Centre for Microbiology (ISCIII) for virological study. Respiratory virus detection was performed by four independent real-time multiplex PCR (RT-PCR) assays using the SuperScript III Platinum One-Step Quantitative RT-PCR System (Invitrogen®, Waltham, MA). The first assay detected Influenza A, B, and C viruses; the second assay detected parainfluenza viruses 1 to 4 (PIV), hRV, and enteroviruses; the third assay detected RSV types A and B, human metapneumovirus (hMPV), human bocavirus (hBoV), and AdV. Human coronavirus (HCoV) was investigated using a generic RT-PCR that was able to detect human alpha and beta coronaviruses, HCoV 229E/HCoV NL63, and HCoV OC43/HCoV HKU1, respectively. The primers and Taqman probes used in the study had already been reported by the study investigators In addition, detection of SARS-CoV-2 was performed on an extracted RNA from NPAs from 2020 using a real-time RT-PCR assay based on the method designed by Corman et al., for the specific amplification of the E gene using the One-Step RT-PCR Kit (NZYTech, Lisbon, Portugal).
Total DNA was extracted from paired NPAs (500ul) and faecal material (approx. 100 mg) using the automated assisted method based on magnetic beads (Maxwell® RSC Instrument coupled with Maxwell RSC Pure Food GMO and authentication kit, Promega, Spain) following the manufacturer's instructions with previous treatments to improve the DNA extraction. In brief, samples were treated with lysozyme (20 mg/mL) and mutanolysin (5 U/mL) for 60 min at 37 °C and a preliminary step of cell disruption with 3-μm diameter glass beads during 1 min at 6 m/s by a bead beater FastPrep 24-5 G Homogenizer (MP Biomedicals). Purification of the DNA was performed using DNA Purification Kit (Macherey-Nagel, Duren, Germany) according to manufacturer's instructions and DNA concentration was measured using Qubit® 2.0 Fluorometer (Life Technology, Carlsbad, CA) for further analysis.
Amplicon libraries of the 16S rRNA gene were generated with a protocol from Illumina (V3-V4 variable region), were treated as described in the Illumina protocol, and in Cabrera-Rubio R et al. Samples were sequenced on Illumina paired-end (2×250 bp) on a NovaSeq- PE250 Illumina platform (Novogene Bioinformatics Technology Co., Ltd) according to manufacturer instructions. Controls during DNA extraction and PCR amplification were also included and sequenced.
16S rRNA gene sequencing data were processed with a modified pipeline described by Cabrera-Rubio et al. Sequencing data were quality-filtered (removal of low- quality nucleotides at the 3' end in windows of 20 nt, sequences with less than 200nt, adapter removal and deduplication) with fastp programme. After that, sequences were joined with VSEARCH, and were denoised with Deblur, with default parameters. Amplicon Sequence Variants (ASVs) mapping to the human genome (GRCh38) using the Burrow- Wheeler Aligner in Deconseq v0.4.3 were filtered out. The R package Decontam was applied to control for potential contaminants. All samples with less than 5000 reads were eliminated for the analysis. ASVs were annotated with a Naïve-Bayes classifier based on the scikit-learn system and the RDP database. The ASVs were aligned with MAFFT, to then make a phylogenetic tree with FASTTREE, that was then midpoint-rooted.
Since samples were obtained from two different NICUs (HSO and HULP), potential centre-specific variability was taken into account in all statistical analyses. Centre was included as a stratification or grouping factor in the alpha and beta diversity analyses, differential abundance testing, and clinical outcome comparisons (e.g., wheezing). When appropriate, analyses were also performed separately within each hospital to confirm the robustness of the results. Differential abundance analyses didn't reveal distinct microbial signatures by hospital (with statistical significance). This approach allowed us to control for potential intra-hospital dependencies and hospital-related confounding, as recommended in previous multi-centre microbiome studies.
Binary clinical variables such as sex and delivery mode were encoded as numeric values (0/1) prior to modelling, in order to facilitate integration with continuous microbiome predictors in the random forest algorithm, which was used to predict the binary outcomes of interest. While this numeric encoding can partially reduce the known variable importance bias in standard random forests, it does not fully eliminate it, especially when compared to variables with multiple split points or continuous distributions.
To address this, we additionally implemented conditional inference forests using the cforest function from the party package in R, following the unbiased recursive partitioning framework described by Hothorn et al. 2006. This method reduces selection bias by separating variable selection and splitting criteria and is more appropriate when working with mixed-type predictors. Variable importance was calculated using conditional permutation (varimp with conditional = TRUE), and both sets of results (standard and conditional random forest) were compared. In both cases, the random forest models were applied to evaluate the predictive contribution of clinical and microbiome variables. The conditional model confirmed the robustness of the microbial predictors identified in the original analysis, while providing a less biased estimate of the relative importance of clinical covariates.
In parallel, to further characterize predictive microbial signatures, we applied the microbial balance approach using the selbal.cv() function from the selbal R package. This function implements repeated k-fold cross-validation to identify robust balances of microbial taxa associated with a binary outcome. We used 5-fold cross-validation repeated 10 times (n.fold = 5, n.iter = 10) to ensure robustness. At each iteration, a logistic regression model was trained on the training folds, and the area under the receiver operating characteristic curve (AUC) was computed on the held-out fold to assess out-of-sample predictive performance. The optimal microbial balance was selected based on the highest average AUC across all iterations.
Since samples were obtained from two different NICUs (HSO and HULP), potential centre-specific variability was taken into account in all statistical analyses. Centre was included as a stratification or grouping factor in the alpha and beta diversity analyses, differential abundance testing, and clinical outcome comparisons (e.g., wheezing). When appropriate, analyses were also performed separately within each hospital to confirm the robustness of the results. This approach allowed us to control for potential intra-hospital dependencies and hospital-related confounding, as recommended in previous multi-centre microbiome studies.
Multivariate analyses assessing the relationship between microbiota features and wheezing were adjusted for key clinical covariates, including gestational age (continuous), presence of BPD (yes/no), and antibiotic exposure (categorized by timing and duration). These variables were selected based on prior knowledge of their potential influence on both microbiota development and respiratory outcomes.
Alpha-diversity indices (Chao1, Simpson and Shannon, and differences with Anova tests) and beta diversity (Bray-Curtis and Unifrac, with PERMANOVA with vegan R package were obtained using phyloseq R package. Rarefied was not performed, but a filter of 0.01% relative abundance in a minimum of 25% of the samples was established to perform the following analyses. To determine the most appropriate statistical method, we first applied the Shapiro-Wilk test to assess whether the data followed a normal distribution. The data were not normally distributed. Since the results indicated that the data were non-parametric, we employed Spearman's rank correlation coefficient to evaluate the association between variables. Spearman correlations were calculated using Phyloseq in R programme version 4.1.2. Pairwise Spearman's rank correlations were performed to assess potential ecological relationships between bacterial genera and to explore associations with clinical variables. Given the non-normal and zero-inflated distribution of microbiota data, even after log transformation, Spearman's correlation was selected over parametric alternatives. Binary clinical variables (e.g., delivery mode, antibiotic use) were numerically coded (0/1) for inclusion in the analyses. Locally Estimated Scatterplot Smoothing curves were applied to visualize monotonic trends without assuming linearity. This non-parametric approach is widely used in exploratory microbiome analyses, particularly when integrating binary and continuous predictors, and serves as an efficient first-pass screening tool to complement downstream multivariate methods (e.g., RDA, LEfSe, and Random Forest). These analyses were exploratory, unadjusted for covariates or multiple testing, and should therefore be interpreted with caution. All correlation analyses and visualizations were conducted using the ggpubr R package version 0.6.0. Scatterplots included in correlation analyses display trend lines for visual guidance only, and LOESS smoothing curves with 95% confidence intervals were applied to better illustrate monotonic trends. Spearman's rank correlation coefficients are reported, as linearity was not assumed. Associations between microbiota and wheezing were analyzed via balanced selection. In addition, a linear discriminant analysis effect size (LefSe); Segata et al. was performed to discover specific bacterial biomarkers. All P values were adjusted using False discovery rate (FDR), based on the Benjamini and Hochberg (BH) method.
In addition to univariate comparisons, multivariate logistic regression models were used to evaluate the association between individual genera and wheezing development, adjusting for relevant clinical covariates including gestational age, sex, delivery mode, antibiotic use, and respiratory support. These models allowed us to account for potential confounding and assess the independent predictive value of microbial taxa. Odds ratios (OR) and 95% confidence intervals (CI) were calculated to aid interpretation.