Characterization of longitudinal nasopharyngeal microbiome patterns in maternally HIV-exposed Zambian infants

Sample processing and storage

NP swabs were obtained from the posterior nasopharynx using a sterile flocked tipped nylon swab (Copan Diagnostics, Merrieta, California). The swabs were then placed in universal transport media, put on ice and transferred to our onsite lab on the same campus, where they were aliquoted and stored at -80°C until DNA extraction. DNA was extracted using the NucliSENS EasyMagG System (bioMérieux, Marcy l’Etoile, France). Extracted DNA was stored at our lab located at the University Teaching Hospital in Lusaka at -80°C. Sample collection, processing and storage were previously described (Gill et al., 2016).¹⁸

16S ribosomal DNA amplification and MiSeq sequencing

For 16S library preparations, two PCR reactions were completed on the template DNA. Initially the DNA was amplified using universal bacterial primers²⁰ specific to the V3–V4 region of the 16S rRNA gene.²¹ Library preparation was performed according to the standard instructions of the 16S Metagenomic Sequencing Library Preparation protocol (Illumina, USA). The 16S primer pairs incorporated the Illumina overhang adaptor (16S forward primer 5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG-3’; 16S reverse primer 5’-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC-3’)

Each PCR reaction contained DNA template (~12 ng), 5µℓ forward primer (1μM), 5 µℓ reverse primer (1μM), 12.5 µℓ 2 X Kapa HiFi Hotstart ready mix (KAPA Biosystems Woburn, MA), and PCR grade water to a final volume of 25µℓ. PCR amplification was carried out as follows: heated lid 110°C, 95°C for 3 min, 25 cycles of 95°C for 30s, 55°C for 30s, 72°C for 30s, then 72°C for 5 min and held at 4°C. Negative control reactions without any template DNA were carried out simultaneously.

The size of the amplicons was then visualized using the 4200 TapeStation (Agilent Technologies, Germany). Successful PCR products were cleaned using AMPure XP magnetic bead-based purification (Beckman Coulter, IN). The IDT for Illumina Nextera DNA UD Indexes kit (Illumina, San Diego, CA) with unique dual index adapters were used to allow for multiplexing. Each PCR reaction contained purified DNA (5 μℓ), 10 μℓ index primer mix, 25 μℓ 2X Kapa HiFi Hot Start Ready mix and 10 μℓ PCR grade water. PCR reactions were performed on a Bio-Rad C1000 Thermal Cycler (Bio-Rad, Hercules, CA) Cycling conditions consisted of one cycle of 95°C for 3 min, followed by eight cycles of 95°C for 30 s, 55°C for 30 s and 72°C for 30 s, followed by a final extension cycle of 72°C for 5 min. PCR products of negative controls were confirmed negative on Agilent TapeSataion (no band observed).

Prior to library pooling, the indexed libraries were purified with Ampure XP beads and quantified using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA). Purified amplicons were run on the Agilent TapeStation (Agilent Technologies, Germany) for quality analysis before sequencing. The sample pool (2 nM) was denatured with 0.2N NaOH, then diluted to 4 pM and combined with 10% (v/v) denatured 20 pM PhiX, prepared following Illumina guidelines. Libraries were then sequenced on the Illumina MiSeq sequencing platform (Illumina, USA) at the Sequencing Core Facility, National Institute for Communicable Diseases (NICD) of the National Health Laboratory Service, South Africa, using a 2 x 300 cycle V3 kit, following standard Illumina sequencing protocols. Negative controls were sequenced as well, resulting in extremely low reads that were not further analyzed.

In addition to using negative controls, all samples were processed at random to account for reagent contamination. Lab technicians were blinded to the timing of sample collection and clinical data.

We assessed the quality of the sequencing data using FastQC v0.11.9.²² Trimmomatic²⁰ v0.39 was used to trim Illumina adapters and remove low-quality sequences. We performed a sliding window trim, cutting once when the average quality score within a window of six bases falls below 15. We removed both leading and trailing low quality or N bases below quality six. All other parameters used the default settings.

Sequencing data were aligned to bacterial genomes and profiled using the PathoScope 2.0 pipeline.²³ All RefSeq representative bacterial genomes available as of November 2, 2018 were used as a PathoScope reference library. We obtained target read counts delineated by National Center for Biotechnology Information (NCBI) unique identifiers (UIDs) and matched them to the NCBI Taxonomy database to retrieve accurate taxonomic hierarchy information. We then aggregated reads by genera. Data were transformed to relative abundances using counts per million (CPM) and normalized using log CPM for subsequent analyses. In most cases, taxa belonging to genera with average relative abundances of less than 1% were grouped as “Other” in analyses. All code has been made available via Zenodo.²⁴

Statistical analysis

All analyses were performed using R Statistical Software (v4.2.1; R Core Team 2022).

Microbial abundances across sample groups were visualized using alluvial diagrams and stacked bar plots using the R packages ggplot2 v3.3.5 and alluvial v0.2-0 . The alluvial diagrams illustrate individual genera as stream fields that change position at different time points. The height of a stream field represents the relative abundance of that taxon. At a given time point, stream fields are ranked from the highest to lowest abundance (top to bottom). These were plotted for infants according to HIV exposure status over several time points. Stacked bar plots were used to visualize the relative abundance of microbes at a given taxonomic level in each sample, represented as a single bar, labeled by time point, and plotted within each HIV exposure status group for separate mothers and infant comparisons. These plotting and diagramming techniques allow for an efficient overview of the types of differences inherently present in the data at the group level.

Generalized estimating equations (GEEs) as described in Liang and Zeger (1986)²⁵ and extended by Agresti (2002)²⁶ have been widely used for modeling longitudinal data,²⁷ and more recently for longitudinal microbiome data.^28,29 For each genus present in the microbial aggregate of samples, we modeled normalized log CPM relative taxon counts, estimating the effects of time point and HIV exposure status and their interaction, while accounting for the underlying structure of clusters formed by individual subjects. We defined the responses Y₁, Y₂,..., Y_n as the collection of infant relative abundances for a given taxon in log CPM, with n = 129. We identified the mean model µ_ij for the i^th infant and j^th timepoint. With regression parameters β_k representing HIV exposure status and time point, and the AR(1) variance structure V_i, we formed the estimating equation:

This then becomes an optimization problem, such that solving U(β) = 0 estimates the parameters β_k. We modeled abundances for all 12 genera and nine of the top species. The link function g was chosen to be a Gaussian link. It was assumed that these abundances are correlated within infants for the various sampling time points. As such, we accounted for this with a first-order autoregressive AR(1) working correlation structure with homogenous variances such that the correlation between adjacent time points was assumed to be more similar. Parameter estimates were collected from each model, along with Wald test p-values. Multiple testing was adjusted for using a Bonferroni correction. Models were created in R using geepack v1.3.10.³⁰

Hotelling’s T² tests³¹ were used to determine whether the microbiome profiles exhibited notable differences or trends across time and groups. We then used t-tests to identify which genera contributed most to these differences. Groups of mothers and infants or HEU and HUU infants are designated as the two sampling units on which the relative abundances of the p most abundant genera were measured. For paired tests, we chose p = 6 variables to ensure that n < p so that singularity could be avoided and T² could be properly computed, where n is the number of measurements in a sampling unit. We generally tested the hypotheses

which, in the paired case given μ_diff = μ_y – μ_x, is equivalent to

to conduct a comparison between groups for the p most abundant genera for testing groups. The population means μ_x and μ_y represent the two sampling units of a given test, which were generally either μ_mothers and μ_infants or μ_HUU and μ_HEU. We treated samples as paired in both cases, such that the mothers are paired with their own infants, or HEU and HUU infants were paired according to the pre-analysis matching carried out previously. When conducting testing on data from HIV positive mothers and HIV negative mothers, we did not pair mothers and instead relied on a standard two-sample test. We assumed that the relevant conditions for testing are met, meaning that the groups are correlated and have a multivariate normal distribution. Normality is met by using microbe abundances in log CPM units.

We conducted three main groups of testing with the T² statistic. The first paired test was between mother-infant pairings at the first or last time point (t ∈ [0, 6]) and for either the mother/infant HIV(+)/HEU or HIV(-)/HUU groups for a total of four tests. These tests offer a clear between-group comparison while accounting for similarity in the mother/infant dyads and calculating genera-specific differences in ruling whether the hypotheses are met. The second testing approach consisted of seven tests. We compared relative abundances of the six most highly abundant genera at each of the seven time points (t ∈ [0, 1,..., 6]) from paired HEU and HUU infants. The final testing approach compared all HIV(+) and HIV(-) mothers’ samples for all 12 genera across mothers, including the “Other” designation as described previously. All samples were unpaired and therefore utilized the unpaired multivariate T² generalization.

Test statistics were calculated to test whether microbiome profiles in paired infants were notably different among the seven time points, and whether HEU mothers and infants had similar profiles at the first and last time points. After calculating the T² test statistic, we calculate an equivalent F_p,n–p test statistic distributed with p and n – p degrees of freedom and use this distribution to calculate the test’s p-value:

Paired or two-sample t-tests were used to distinguish which genera are most important to the differences identified with Hotelling’s T² tests, and α-levels were adjusted after performing the p tests by using a Bonferroni critical value.

Beta diversity using the Bray-Curtis dissimilarity metric was compared using a non-parametric Wilcoxon rank-sum test between HEU and HUU infants, and between HIV(+) and HIV(-) mothers. Tests were conducted separately for each subset of samples for the first and last time points for a total of four tests. The null hypothesis conjectures that the distributions of both populations are equal, under the general assumption that all observations from both groups are independent of each other.

Sample characteristics

Counts of infant and mother samples available at each time point and for each group of interest can be found in Table 2. Although we included 10 HUU and 10 HEU infants and their mothers in the study, there was some small variation in numbers of samples at the different time points, but they were overall close in evenness. Infants were approximately 7 ± 1 (mean ± SD) days old at time of first swabbing, and were on average 22 ± 2 days, 42 ± 3 days, 58 ± 2 days, 74 ± 2 days, 88 ± 2 days, and 105 ± 2 days respectively at subsequent time points 1–6.

Longitudinal differences in the nasopharyngeal microbiome composition

Our FastQC analysis indicated that the overall sequencing quality was excellent, with mean Phred quality scores remaining greater than 25 (99.5% accuracy) for at least 175 bp for both forward and reverse reads. Trimmomatic removed less than 3% of reads in any given sample. The analysis covered 129 infant and 38 mother swab samples, with an average of 124,905 ± 299,120 (mean ± SE) reads per infant sample (max = 3,016,276; min = 13,218) and 70,463 ± 48,730 reads per mother sample (max = 208,029; min = 14,339). The number of reads in infants’ samples was significantly higher than that of mothers’ swabs (p = 0.02), which may in part be the result of mothers’ acquired immunity over time and therefore lower overall NP carriage. In our raw data, we uniquely identified 17 phyla, 647 genera, and 758 total species across all samples. NP microbiome composition separated according to HIV exposure status are depicted in Figures 1A and 1B for infants and mothers, respectively. The alluvial plot in Figure 1C graphically depicts the most abundant genera in log CPM present in HEU and HUU infants. As relative abundances (in log CPM) change over time, the position of a flow stream representing a single genus may change position relative to the other genera. At a given time point, the flow streams are stacked according to the abundance relative to other streams.

Figure 1.

The maturation over 120 days of the NP microbiomes of A) healthy HIV-exposed, uninfected (HEU; n=10) and HIV-unexposed, uninfected (HUU; n=10) infants and B) HIV(+) (n=10) and HIV(-) mothers (n=10). These stacked bar plots reveal variation in the relative abundance of microbes between groups of either infants or mothers clustered at the genus level. Each bar represents a single time point binned by age and is the average of ~10 samples. Genera with an average relative abundance of <1% across all samples are labeled as “Other.” The alluvial plot in C) depicts the changing presence of genera across all infant samples by HIV exposure status. As relative abundances change over time, the position of a flow stream representing a single genus may change position relative to the other genera.

Due to the heavy presence of taxa with relatively low abundance, we labeled taxa with average relative abundances of less than 1% as “Other” at the genera and species level. Post-grouping, the present taxa were limited to three distinctly identifiable phyla, encompassing 12 genera and 87 species. The most abundant phyla were the Firmicutes (~56.5%), Proteobacteria (~22.4%), and Actinobacteria (~16.1%). The remaining 5% of reads were characterized as “Other.” At the genus level, the most abundant groups were the Dolosigranulum (~23.5%), Staphylococcus (17.2%), Corynebacterium (~16.1%), Streptococcus (~13.5%) and Moraxella (~12.4%) genera.

For both the HEU and HUU infant groups, the respiratory microbiome during the first months of life is dominated by Staphylococci and Corynebacteria. Early on, we observed the emergence of more typical respiratory bacteria such as Moraxella and Streptococcus sp. Additionally, the commensal Dolosigranulum emerges as a dominant member of the microbiome within the first weeks of life. We saw the appearance of Haemophilus somewhat later at around four to six weeks. Overall, these transitions occurred in an orderly and stepwise pattern. These transitional patterns align with those found in a longitudinal East Asian infant cohort.³²

While longitudinal trends were strongly apparent across infant groups, the observed differences between the HEU and HUU infants were subtle. Figure 1A illustrates higher amounts of Streptococcus and suppression of Dolosigranulum among HEU vs. HUU infants. Furthermore, Figure 1C suggests increasing amounts of Streptococcus and Haemophilus over time for HEU infants, whereas Staphylococcus largely declines after time point three in both groups. Together these findings support a time-dependent effect on genus abundance taking place soon after birth.

As this is one of the few studies conducted on the NP microbiome in Zambia, there are no conclusive baseline expectations for a typical microbiome composition. However, we have compiled the results of a published longitudinal study that analyzed NP swabs from a cohort of 234 healthy infants from Washington, D.C.³³ (Table 3). Swabs in the study were taken at ~two, ~six, and ~12 months. We have included the proportion of genera present across the aggregate of all swabs, and for the ~two-month sampling (no further information on exact timing of swabs was specified). Slight differences appear to be present, which may be in part due to differences in living conditions, sample size, and choice of analysis database used. The Teo et al. study utilized the Greengenes database, which produces far less sensitive results when compared to other databases like RefSeq and Silva.³⁴ Otherwise, when comparing the Teo et al. cohort at two months to those in our study, we find that the abundance of Staphylococcus is nearly twice as high as in our study for either HEU or HUU infants, although the Corynebacterium, Moraxella, Streptococcus average means are all within a few percentage points for both groups. In our analysis, the most overwhelming difference between the HEU and HUU infants is in Dolosigranulum, for which HUU infants showed a 19.63% higher presence of Dolosigranulum when averaged across time points. Overall, as this study is not a one-to-one comparison to the Teo et al. study, it is difficult to draw exact conclusions, but it appears that the NP microbiome profiles observed here are similar to those seen in healthy infants of a similar age.

Relationship between nasopharyngeal microbiome composition and HIV exposure

GEEs revealed some differences among genera and species for time point and HIV exposure status when adjusting for subject variation and the interaction effect. Models were created for all 12 genera (including the “Other” grouping) and for the nine species that averaged greater than 1% abundance across all infant samples. Of all taxa tested, four of nine species and six of 12 genera exhibited substantial instability across time points (Table 4). One species, Staphylococcus haemolyticus, was highly associated with HIV exposure in infants (p <0.01, adj. p = 0.01). Additionally, Streptococcus mitis indicated a strong interaction effect between time and HIV exposure status (p <0.01, adj. p = 0.04). Figure 2 shows the estimated marginal means of the time and HIV-exposure status effects for S. haemolyticus, S. mitis, Haemophilus influenzae, and the Dolosigranulum genus. We chose to include these microbes in the figure as a representative selection of microbes with significant and non-significant effects. In Figure 2a, differences in the abundance of S. haemolyticus for HEU and HUU infants is cleanly pronounced for several time points, reflecting the Wald test result, but this is not the case for the other microbes in Figures 2b, 2c, or 2d. All of these had at least marginally significant p-values when testing the time point effect but lacked strong HIV exposure status effects.

Figure 2. Plots of estimated marginal means of the abundance of a given microbe accounting for estimated GEE effects.

The time of sampling is on the x-axis, and the abundance of the microbe in log CPM is on the y-axis. The estimated change in abundance over time is illustrated for control infants (in red) and HIV-exposed, uninfected (HEU) infants (in blue). The interaction between time and HIV exposure status is accounted for in the differing group slopes. All microbes had at least marginally significant time effects, but only S. haemolyticus had a very strong HIV exposure status effect.

Hotelling’s T² test statistics were used to compare relative abundances of the top six most abundant genera between paired HEU and HUU infants at each time point, namely Dolosigranulum, Streptococcus, Moraxella, Staphylococcus, Corynebacterium, and Haemophilus. As previously noted, we chose p = 6 variables to ensure that n < p so that singularity could be avoided and T² could be properly computed, where n is the number of measurements in a given sampling unit. Paired infants did not present significantly different microbiome profiles at any of the time points at α = 0.05 (Table 5). The largest differences seemed to occur at time point 5 (p = 0.07). The fluctuating critical F values indicate that uneven sample distribution across time points contributes to the difficulty in identifying clear differences between populations.

Overall microbiome composition varied significantly at the times of first and last sampling for HUU and HEU infants. Among all infant samples at time point 0, a Wilcoxon rank-sum test of the Bray-Curtis dissimilarity revealed a non-significant difference in microbiome composition between HUU and HEU infants (p = 0.68). The same test for all infant samples at time point six was contrastingly showed large differences in beta diversity (p <0.01). Within-group, between-subject tests revealed variation at the first and last time points within the HEU group, and at both first and last time points within the HUU group (all p-values < 0.01). However, large variation was not present between time points within the HEU (p = 0.11) or HUU infant samples (p = 0.78). As a result of heavy inter-subject variability, alpha diversity metrics were not deemed useful for this analysis.

Relationship between infant and mother nasopharyngeal microbiome composition

We also investigated the relationship between the NP microbiome in infants and their mothers. Paired Hotelling’s T² tests were again used to test the log CPM of genera as a measure of relative abundance at the first and last time points for the HIV(+)/HEU and HIV(-)/HUU groups (Table 6). Genera tested were the same as listed for the infant-specific Hotelling’s tests, delineated previously. HIV(+)/HEU mother-infant pairs had marginally significant profile differences at the first time point ($T_{6,4}^2;$ p = 0.10) but were similar enough to not present as having significant profiles at the last time point ($T_{6,1}^2;$ ; p = 0.60). HIV(-)/HUU mother-infant pairs also had marginally significant differences at the first time point ($T_{\mathit{6},\mathit{4}}^{\mathit{2}};$ p = 0.10) but varied notably at the last time point ($T_{6,3}^2;$ p = 0.02). Paired t-tests for the six genera were separately conducted for HIV(-)/HUU mother-infant pairs to distinguish which genera are most important to the identified difference at that time point. The largest difference in abundance occurred for the Haemophilus (t₈, p = 0.02; adj. p-value = 0.10) and Staphylococcus genera (t₉; p-value <0.01; adj. p-value <0.01). HUU Infants were noted as having greater Haemophilus carriage than their HIV(-) mothers, whereas mothers had greater Staphylococcus carriage than their infants.

As part of our inquiry into why HEU infants might be more susceptible to higher morbidity or mortality, we conducted paired t-tests for well-known pathogenic species between mother-infant pairs for Streptococcus pneumoniae, Haemophilus influenzae, and Staphylococcus haemolyticus. Tests compared pairs at the first and last time points, and within either HIV or control subgroups, for a total of four different tests per species (Table 7). Inequality in pathogen carriage was ascertained by the mean of differences in log CPM. The tests for S. pneumoniae were nominally significant at t = 0 for the HIV(+)/HEU group (p = 0.03, adj. p = 0.4) and t = 6 for the HIV(-)/HUU group (p = 0.01, adj. p = 0.1). In the former, the mothers had more pathogen carriage than the infants, but this was the opposite case for the latter. H. influenzae was more highly abundant in HUU infants than in their HIV(-) mothers at t = 6 (p = 0.01, adj. p = 0.1), and S. haemolyticus was likely to be found in the HUU infants at t = 0 in a higher concentration than in their mothers.

Trends in NP microbial community composition between HIV(+) and HIV(-) mothers

We compared HIV(+) and HIV(-) mothers at the time points that we sequenced. We used the Bray-Curtis dissimilarity metric to compare compositional dissimilarity between NP microbiomes. We identified strong dissimilarity between HIV(+) and HIV(-) mothers at time point six (p <0.01), but higher similarity at the first time point (p = 0.70). Additionally, the summed abundances of all twelve genera, including “Other” genera, were used to conduct a two-sample Hotelling’s T² test on unpaired mothers. The differences we observed were marginally significant at the first time point ($T_{12,7}^2;$ p = 0.16), but more noticeable at the latter ($T_{12,5}^2;$ p = 0.04). Unpaired t-tests indicated that the “Other” taxa (p = 0.03, adj. p = 0.37), Alkalihalophilus (p = 0.08, adj. p = 1), and Paracoccus (p = 0.10, adj. p=1) were the largest contributors to this difference. These taxa were all more highly abundant in the HIV(-) mothers.