Resident microbes of lactation rooms and daycares

Sample collection

Samples were collected from 28 University of California at Davis Lactation Rooms quarterly, excluding summers, between Spring 2015 and Fall 2016 and 4 daycares located in Davis, California in Winter 2016 and Fall 2016. The Spring 2015 quarter was from 26 March 2015 to 11 June 2015. The Fall 2015 quarter was from 21 September 2015 until 11 December 2015. The Winter 2016 quarter was from 4 January 2016 until 19 March 2016. The Spring 2016 quarter was from 24 March 2016 until 9 June 2016. The Fall 2016 quarter was from 19 September 2016 until 9 December 2016. Access to daycares for sampling was only available during the Winter and Fall 2016 quarters. All included rooms were within UC Davis property, so formal permission was not required. At the time, Shirley German managed the lactation rooms, Christine Pretti was the Director of Hutchison Child Development Center, Shelby Faria was the Director at La Rue Park, Janet Thompson was the Director of the Early Childhood Lab School at the Center for Child and Family Studies, and Tanya Chordas was the Director of the Russell Park Child Development Center. To collect samples, sterile swabs were first dampened with sterile phosphate buffered saline, then rubbed over the sampling location to collect any microbes present. Five samples per lactation room per sampling visit were collected from each lactation room; one each from a chair, a table, a pump, the floor, and the door to the room. Lactation rooms could be divided into categories based on other uses of the room, including sole purpose lactation rooms, a resident room, women’s restrooms, restroom lounges, single toilet unisex restrooms, and a shower room. There were five locations targeted for sample collection in daycares: a chair, a table, a toy, the floor, and the door to the room. After collection, swabs were placed in sterile 1.5 mL cyrovial tubes and stored at −80 °C until DNA extraction.

In addition, for a subset of ten lactation rooms, sensor nodes were placed near the door to enter the room to monitor room occupancy and select environmental conditions 24/7. Each custom-built sensor node (microcontroller, Cypress Semiconductor, PSoC5, CY8C5267AXI-LP051) contained a temperature and humidity sensor (Honeywell, HIH6130-021-001), a beam-break sensor (PIR Mini Sensor, Parallax 28033), a microSD card (Kingston 8 GB Class 4 MicroSDHC Card Flash Memory with SD Adapter SDC4/8 GB), and a battery (LiPo, 3.7v 2,500 mAh; Hunan Sounddon New Energy Co., Xiangtan City, China). The temperature and humidity data was collected every 15 min and stored on the microSD card. The beam-break sensor was used to determine when the room was in use by storing a time-stamp on the microSD card whenever the door was opened or closed. The microSD cards were retrieved from the sensor nodes to review the stored data. The sensor nodes were designed to be in an ultra-low-power state when measurements were not being taken to extend the battery life as long as possible, permitting measurements to be taken every 15 min for the entire time the sensors were in place. These sensor nodes were in place from early Fall 2015 to Fall 2016.

DNA extraction and 16S rRNA gene sequencing

The base of the swabs were cut off using sterile scissors, and DNA was extracted using the Zymo Research Quick-DNA™ Fecal/Soil Microbe Miniprep Kit (Zymo Research, Irvine, CA, USA) according to the manufacturer’s instructions. For a negative control, seven swabs not used to collect samples were processed interspersed with the samples. An additional 24 kit control samples were extracted by completing all extraction steps without using a swab sample. Kit control samples were done using the same kit type and methods, but different lot numbers of reagents as kit controls were extracted, sequenced, and analyzed after the original sample extraction and sequencing were complete.

The barcoded F515 (5′-NNNNNNNNGTGTGCCAGC MGCCGCGGTAA-3′) and R806 (5′-GGACTACHVGGGTWTCTAAT-3′) primers were used to amplify the V4 region of the 16S ribosomal RNA gene, as previously described (Huda et al., 2014). After amplification, samples were purified using the Qiagen QIAquick PCR Purification Kit (Qiagen, Hilden, Germany) and pooled. The pooled library was then sequenced by the UC Davis Genome Center DNA Technologies Sequencing Core using 250-bp PE reads on the Illumina MiSeq DNA sequencing system (Illumina, San Diego, CA, USA). At least one negative control sample was included with every 95 samples.

The sequencing results were demultiplexed using Sabre (Joshi, 2011) and imported into QIIME2-2017.12 (Caporaso et al., 2010), which was used for trimming, denoising, taxonomic classification, and phylogenetic tree building. Bases before basepair 22 and after basepair 240 were trimmed from the forward read with bases before basepair 23 and after basepair 249 were trimmed from the reverse read. The trimmed reads were processed with DADA2 (Callahan et al., 2016). Taxonomy was assigned to representative sequences for each amplicon sequence variant (ASV) using the 99% Greengenes (McDonald et al., 2011) pre-generated naïve bayes classifier, and a phylogenetic tree was built using fasttree (Price, Dehal & Arkin, 2009) and midpoint rooted. The resulting ASV table, phylogenetic tree, representative sequences, and taxonomic assignments were exported from QIIME2 for additional analysis in R. These same steps were completed for the kit control samples, and any reads from ASVs found in both the samples and the kit controls were excluded from additional analysis.

Statistical analysis

All statistical analysis was completed in R 3.4.3 statistical software (R Core Team, 2016). Samples with less than 1,000 reads were excluded from all analysis, as the rarefaction curve suggested that the number of species detected per sample reached a plateau for most samples at a read depth of approximately 1,000 (see Fig. S2). Prior to analysis, all samples were rarefied to a read depth of 1,040 reads, as this was the read depth of the sample with the fewest reads after excluding samples with <1,000 reads. For alpha and beta diversity analyses, samples were stratified by quarter (Spring 2015, Fall 2015, Winter 2016, Spring 2016 and Fall 2016). Only the earliest successfully sequenced sample from each sampling site (e.g., sample collected at a specific location in a room) in each quarter was included in the analysis. The vegan 2.5.1 package (Oksanen et al., 2017) was used to calculate the Shannon index. Kruskal–Wallis test was used to test for differences in alpha diversity by room type and collection location, followed by Dunn’s test if the Kruskal–Wallis test was significant to determine which room types and collection locations differed from each other. A two-tailed p < 0.05 was considered significant for the Kruskal–Wallis test and a Bonferroni adjusted two-tailed p < 0.05 was considered significant for Dunn’s test. For beta-diversity analysis, the weighted and unweighted UniFrac distance between samples was calculated using the GUniFrac 1.1 package (Chen et al., 2012). The vegan package implementation of non-metric multidimensional scaling (NMDS) was used to generate ordination plots, if stress was less than 0.2, a 2 dimensional plot was used. Otherwise the number of dimensions was increased to obtain a plot that was a reasonable representation of the data. Beta dispersion was calculated using the betadisper function in the vegan package (Dixon, 2003). ANOVA was used to test for differences in beta dispersion by room type and by sample collection location. If there were no significant differences in beta dispersion, PERMANOVA as implemented by adonis was run to test for statistical significance in beta diversity by room type and sample collection location, with the interaction term included in the model. As the study design was unbalanced for all quarters, PERMANOVA was not run if the beta dispersion differed significantly by room type or sample collection location if the group with smaller sample size had greater dispersion as this will tend to make PERMANOVA too liberal (Anderson & Walsh, 2013).

To explore the relationship between temperature and humidity and the built environment microbiome, the collection day and room of all successfully sequenced samples was compared to the days and rooms with temperature and humidity sensor. For each sample collected while a sensor was present in the room, the average temperature and humidity for the day of sample collection was calculated. To test for an association between alpha diversity and temperature and humidity, linear mixed effects models were run using the lme command in the nlme 3.1-137 R package (Pinheiro et al., 2014) with the Shannon index as the outcome variable and the average room temperature and the average room humidity as fixed variables. Room was included as a random variable. For beta diversity analysis, first the same set of first samples was used in an NMDS based on the weighted and unweighted UniFrac distance measures. Then the association of the ordination with the continuous average temperature and continuous average humidity were tested in separate ordisurf models (as implemented in the vegan package).

Finally, we sought to understand if greater changes in temperature or humidity were associated with greater changes with the built environment microbiome. For any room/location combinations that had multiple samples successfully sequenced with available sensor data, the Jaccard distance between the samples between the repeated measurements. Then the absolute value of the difference between the average temperature and the average humidity at the two time points was calculated (magnitude of change in temperature and humidity). The association between Jaccard distance and magnitude of change in temperature and humidity was calculated using a linear mixed effects model (lme command) with room included as random effect to account for multiple sampling locations within the same room.

Samples

There were a total of 460 samples collected from lactation rooms and daycares for this study, seven negative controls, and 24 extraction kit controls. The kit controls had 11 different ASVs present, five of which were present in the other samples. These five ASVs were identified as Escherichia coli, an unknown species of Pseudomonas, Pseudomonas fragi, Pseudomonas veronii, and an unknown species of Lactococcus. These ASVs represented 1.1% of all reads from samples included in the final analysis, ranging from 0% to 6.7% of reads on a per sample basis, with a median of 1% of reads removed from each sample included in additional analysis. One negative control had an unexpectedly high number (>42,000 reads) of reads after sequencing, suggesting contamination in the set of samples either extracted or amplified at the same time as the negative control, as a result all samples handled in a batch with the contaminated negative control were excluded from additional analysis, leaving 388 samples. Of these samples, 230 had a read depth greater than 1,000 reads. Excluding samples with less than 1,000 reads removed all additional negative controls, as after excluding the likely contaminated samples, the negative control with the next most reads had only 48 reads. Of the samples excluded due to low read depth, half had a read depth at or below that observed in the negative controls. An additional 17 samples were excluded from the analysis stratified by time because sampling was completed twice in the same rooms at the same collection locations in some quarters, leaving 213 samples for inclusion in the diversity by quarter analyses. A total of 55 samples were collected from rooms with sensor data available on the day of collection, 39 of these were collected from unique room and collection locations and were included in the alpha and beta diversity analyses. An additional 16 samples were collected from the same room and location as an earlier sample, resulting in 32 samples available for the change in temperature and humidity analysis, and a total of 221 samples included in one or both analyses. Table 1 summarizes the samples included in each analysis. Demultiplexed sequencing data of all samples included in the analysis are available at SRA database with BioProject accession number PRJNA547638. Table 2 summarizes the results by quarter discussed below.

Spring 2015

There were a total of 79 samples successfully sequenced from 25 lactation rooms in Spring 2015. The majority of reads in most samples, regardless of room type or sample collection location, belonged to family Moraxellaceae (Fig. 1). A single sample, the chair sample from room 23, had markedly higher levels of Enterobacteriaceae than all other samples at this time point. There was no significant association between alpha diversity as measured by the Shannon index and sample collection location (Kruskal–Wallis, p = 0.075), or between Shannon index values and room type (Kruskal–Wallis, p = 0.091) (Fig. 2).

There is no visible separation of points on the NMDS plot based on unweighted UniFrac for the Spring 2015 samples by either room type or sample collection location (Fig. 3). There were no significant differences in beta-dispersion by room type (p = 0.59) or sample collection location (p = 0.87). There were significant differences in unweighted UniFrac beta diversity by sample collection location (PERMANOVA, p = 0.043) but not by room type (PERMANOVA, p = 0.059). The room 23 chair sample with high levels of Enterobacteriaceae skewed the weighted UniFrac analysis, with all other points appearing on top of each other in the weighted UniFrac NMDS, therefore, the weighted UniFrac results were not analyzed further. Because remaining quarters generally showed a similar pattern, figures from other quarters can be found in Supplemental Information.

Fall 2015

There were a total of 42 samples successfully sequenced from 13 lactation rooms in Fall 2015. As with Spring 2015, most samples were dominated with reads mapping to Moraxellaceae, however, there were some samples with relatively low levels of Moraxellaceae and more samples with substantial levels of other families present (Fig. S3). There was no significance association between alpha diversity as measured by Shannon index and room type (Kruskal–Wallis, p = 0.11) or sample collection location (Kruskal–Wallis, p = 0.14) (Fig. S4).

There was no visible separation of points on NMDS plot based on unweighted UniFrac for the Fall 2015 samples by room type, however, there did appear to be some separation between the floor samples and the samples collected at other locations (5). There was no significant differences in beta-dispersion by collection location (p = 0.79), but there was a significant difference by room type (p = 0.0072). As the room type with the larger sample size had greater dispersion, PERMANOVA will tend to be too conservative (Anderson & Walsh, 2013) so we proceeded with the PERMANOVA despite the differences in dispersion. As in Spring 2015, there were significant difference in beta diversity as measured by unweighted UniFrac by collection location (p = 0.002), but now there was also a significant difference by room type (PERMANOVA, p = 0.036) as well. Unlike the unweighted UniFrac ordination, the weighted UniFrac ordination did not appear to have any separation of points by room type or collection location (Fig. S6). There were no significant differences in beta-dispersion by room type (p = 0.10) or by collection location (p = 0.46). Despite the lack of visual separation on the NMDS, there was a significant difference in beta diversity by room type (PERMANOVA, p = 0.001). However, there was not a significant difference by sample collection location (p = 0.11).

Winter 2016

There were a total of 40 samples collected and successfully sequenced from 14 lactation rooms and eight samples collected and successfully sequenced from two daycares in Winter 2016, for a total of 48 samples. As with previous quarters, Moraxellaceae was the dominant taxa in most lactation room samples. Three of the 8 day care samples were also dominated by Moraxellaceae, but the remaining samples appeared more diverse (Fig. S7). Consistent with this result, there was a significant difference in alpha diversity as measured by the Shannon Index by room type (Kruskal–Wallis, p = 0.0068) (Fig. S8). Post-hoc comparison confirmed that the significant differences were between daycare samples and dedicated lactation room samples (Dunn’s test, p = 0.0037) and between daycare samples and unisex restroom samples (p = 0.013). There was no significant difference in Shannon’s index by sample collection location (Kruskal–Wallis, p = 0.90) (Fig. S8).

Daycare samples tended to visually separate from samples from other room types on the unweighted UniFrac NMDS (Fig. S9). There were no significant differences in beta-dispersion by room type (p = 0.50) or sample collection location (p = 0.18). However, unlike in previous quarters, there was no significant difference in unweighted UniFrac beta diversity by sample collection location (PERMANOVA, p = 0.55), although there was a significant difference by room type (p = 0.002). As with the unweighted UniFrac, the weighted UniFrac NMDS showed some visual separation of the daycare samples from the other room types. There was no visible separation in weighted UniFrac distance by sample collection location (Fig. S10). While there was no significant difference in beta-dispersion by sample collection location (p = 0.46) for the weighted UniFrac analysis, there was a significant difference by room type (p = 0.0078). Furthermore, the room types with the fewest samples had large dispersion and as a result, PERMANOVA testing for differences has the potential to be too liberal (Anderson & Walsh, 2013), and so room type was not included in the weighted UniFrac PERMANOVA model. There was no significant association between beta diversity as measured by weighted UniFrac and sample collection location (PERMANOVA, p = 0.325).

Spring 2016

There were only 11 total samples from 6 lactation rooms successfully sequenced from Spring 2016. As with other quarters, the majority of the samples were dominated by Moraxellaceae (Fig. S11). Because most room types and most sample collection locations had fewer than 3 samples per category, additional statistical analysis was not completed for this quarter.

Fall 2016

There were 26 samples from 10 lactation rooms and 8 samples from 4 daycares that were successfully sequenced in Fall 2016. Moraxellaceae remained the dominant taxa in most lactation rooms, and unlike in Winter 2016, Moraxellaceae was now the dominant taxa in all daycare samples (Fig. S12). Because there was only a single sample from a single stall women’s restroom collected in this quarter, this sample was excluded from additional analysis. There were no significant differences in alpha diversity as measured by the Shannon index by room type (Kruskal–Wallis, p = 0.83) or by sample collection location (p = 0.19) (Fig. S13).

Despite the more comparable levels of Moraxellaceae between daycares and lactation rooms, the daycares did appear somewhat visually separate from the lactation rooms on the unweighted UniFrac NMDS. There also seemed to be some separation by sample collection location, particularly between pump and floor samples (Fig. S14). There were no significant differences in beta-dispersion by collection location (p = 0.98) but there were significant differences by room type (p = 0.0078), but the room types with more samples had larger dispersion making PERMANOVA conservative. There were significant differences in beta diversity as measured by unweighted UniFrac by both room type (PERMANOVA, p = 0.009) and by sample collection location (p = 0.006). Using weighted UniFrac distance, there was no clear separation by either room type, but the pump collection location did appear to be different from other collection locations (Fig. S15). There was no significant difference in beta-dispersion of weighted UniFrac distance by room type (p = 0.48), however, there was a significant difference in beta-dispersion by sample collection location (p = 0.027) but the location with the greatest dispersion also had the most samples, tending to make PERMANOVA conservative. There was no significant difference in beta diversity based on weighted UniFrac distance by room type (PERMANOVA, p = 0.97), but there was a significant difference by collection location (PERMANOVA, p = 0.042).

Moraxellaceae genera and species present

Because Moraxellaceae was dominant in many sampling locations and time points, further description of the ASVs identified as belonging to this family is warranted. Of all the reads in the rarefied ASV table, 84% were identified as belonging to family Moraxellaceae. The majority of these reads were found in two ASVs. The first ASV was an unknown species of Acinetobacter found in 97% of all samples, and represented 73% of all reads in the rarefied ASV table and 86% of all Moraxellaceae reads. The second ASV was identified as A. johnsonii found in 93% of all samples, and represented 9% of all reads in the rarefied ASV table and 11% of all Moraxellaceae reads. The remaining ASVs were primarily identified as Acinetobacter, but a few ASVs were identified as belonging to other genera (Table 3).

Alpha and beta diversity association with temperature and humidity

The mean average temperature of lactation rooms on the day of sample collection was 20.9 °C (range 18.3–22.4°C) for the samples included in this analysis. The mean average humidity of lactation rooms on the day of sample collection was 35.6% (range 21.2–68.3%) for the samples included in this analysis. There was no significant association between Shannon index values and temperature (p = 0.22) or humidity (p = 0.77).

The unweighted UniFrac NMDS ordination had two dimensions and a stress of 0.0962. There was no significant association between the ordination results and the average room temperature on day of collection (p = 0.16) or the average room humidity on day of collection (p = 0.41) by ordisurf (non-significant, resulting figure not shown). The weighted UniFrac NMDS ordination had two dimensions and a stress of 0.0413. There was no significant association between the ordination results and temperature (p = 0.109) or humidity (p = 0.75) by ordisurf (non-significant, resulting figure not shown).

Community change association with temperature and humidity

A Jaccard distance of 0 indicates identical samples and a Jaccard distance of 1 indicates no overlap of replicon sequence variants between repeated samples from the same location. The Jaccard distance between repeated samples varied widely, with a median value of 0.232 (range 0.040–0.99). There was no significant association between the Jaccard distance and the magnitude of humidity change (lme value = 0.0206, p = 0.083) or between the Jaccard distance and the magnitude of the temperature change (lme value = −0.124, p = 0.51).