Impact of water heater temperature setting and water use frequency on the building plumbing microbiome.

Hot water plumbing is an important conduit of microbes into the indoor environment and can increase risk of opportunistic pathogens (for example, Legionella pneumophila). We examined the combined effects of water heater temperature (39, 42, 48, 51 and 58 °C), pipe orientation (upward/downward), and water use frequency (21, 3 and 1 flush per week) on the microbial composition at the tap using a pilot-scale pipe rig. 16S rRNA gene amplicon sequencing indicated that bulk water and corresponding biofilm typically had distinct taxonomic compositions (R2Adonis=0.246, PAdonis=0.001), yet similar predicted functions based on PICRUSt analysis (R2Adonis=0.087, PAdonis=0.001). Although a prior study had identified 51 °C under low water use frequency to enrich Legionella at the tap, here we reveal that 51 °C is also a threshold above which there are marked effects of the combined influences of temperature, pipe orientation, and use frequency on taxonomic and functional composition. A positive association was noted between relative abundances of Legionella and mitochondrial DNA of Vermamoeba, a genus of amoebae that can enhance virulence and facilitate replication of some pathogens. This study takes a step towards intentional control of the plumbing microbiome and highlights the importance of microbial ecology in governing pathogen proliferation.

Introduction

There is growing interest in the microbiome inhabiting the built environment and its potential implications to human health (Hanski ; Kelley and Gilbert, 2013; Oberauner ). Building plumbing systems deliver water, along with its resident microbes, into homes and other buildings. Most attention has been directed towards opportunistic pathogens (OPs), such as Legionella pneumophila, which are particularly challenging to control as they are native to the fresh water environment. Collectively, OPs residing in building plumbing have become the top contributor of waterborne disease in the U.S. and other developed countries (Hilborn ; Ashbolt, 2015; Beer ), imposing an economic burden over $1 billion per year in the United States alone (Falkinham ). Thus, understanding the composition of the microbiome at the tap, and its contribution of pathogens and general microbiota to built environment and human microbiomes is of great importance.

Hot water plumbing systems are of special interest as a typical route of OPs exposure is inhalation of aerosols and/or direct contact during showering (Feazel ). Elevating the temperature of hot water storage tanks is known to be an effective measure to control L. pneumophila (Darelid ; Buse and Ashbolt, 2011; ANSI/ASHRAE, 2015; Falkinham ), and 60 °C has been demarcated as a regulatory threshold (VHA Directive 1061; Department of Veterans Affairs, 2014). However, motives to reduce energy usage and avoid scalding are in conflict with higher temperature settings (Lévesque ; Rhoads ) and long-term effectiveness of thermal eradication has been questioned (Farhat ). In contrast to a number of studies directly focusing on Legionella in hot water systems (Zietz ; Borella ; Zacheus and Martikainen, 2004; Mathys ; Bargellini ), very few have more broadly examined the microbiome as a whole. Prior studies included field surveys on microbial communities collected from shower head/hose biofilms (Feazel ; Proctor ) and bulk water samples from a hospital hot water system, where the effectiveness of monochloramine disinfection was examined (Baron ). However, a systematic examination of the combined effects of water heater temperature and system design/operational conditions on the actual microbiome composition at the tap has been lacking.

Although hot water storage tanks may remain at elevated temperature settings, the reality at the distal tap is that the biofilm experiences only brief heat shock. Water use frequency at a given tap is thus of particular interest, as it dictates the frequency of heat shock as well as the length of time for the biofilm microbiota to recover from the perturbation, which could be further influenced by pipe orientation and associated convective mixing conditions (Rhoads ). In the bulk water phase, in-building stagnation was found to have a marked effect on the microbial composition of cold water lines (Lautenschlager ; Ji ). Moreover, warmer indoor temperature together with overnight stagnation could increase microbial metabolic activity level and alter the functional profile (Zhang ). We thus hypothesized that different water use frequencies would alter the microbiome composition of both the biofilm and bulk water.

The relationship between the microbial composition of the biofilm and bulk water is also widely debated and likely to be influenced by operating conditions. Biofilm is often perceived to be a reservoir for OPs, where protist hosts can facilitate their replication and, along with the biofilm matrix, protect them from unfavorable conditions such as chlorine/monochloramine disinfection (Pryor ; Cooper and Hanlon, 2010) and low nutrients (Lau and Ashbolt, 2009). However, microbes entrained in the biofilm must be released into the bulk water for relevance to human exposure, which can occur long after the initial entrapment/attachment of the microbe to the biofilm (Howe ). Interestingly, OPs in the two phases do not necessarily respond the same way under a given condition. For example, at room temperature, stagnation predisposed bulk water to Legionella proliferation (Wang ), but resulted in lowest Legionella concentrations in biofilm relative to those under turbulent and laminar flow conditions (Liu ). Such questions call for a mechanistic understanding of the factors shaping the biofilm and bulk water microbiomes of domestic hot water systems.

To overcome complexities encountered in field studies, while maintaining relevance to real-world conditions, we constructed two identical recirculating hot water plumbing rigs with distal taps operated in triplicate under three water use scenarios (21, 3 and 1 flush per week) and two orientations (upward and downward flow). Water heater temperature was elevated in a step-wise fashion from 39 to 58 °C, while an identical control rig remained at 39 °C throughout the 15-month study period. Illumina amplicon sequencing of 16S rRNA genes was applied to map the bacterial and archaeal community composition of influent, recirculating line and distal taps with paired biofilm and first-flush bulk water samples. Metagenomic data was predicted using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt 1.0.0) to characterize hypothetical functionalities of resident microbes. Overall the study takes a step towards identifying to what extent building owners/operators may be able to intentionally control the microbiome composition at the tap through practical measures, and addresses a critical knowledge gap at the intersection of water infrastructure, energy saving and public health.

Materials and methods

System description

Two identical lab-scale hot water plumbing rigs were designed for paired comparison of the impact of water heater temperature and use frequency (Supplementary Figure 1), as reported previously in studies focused on the occurrence of Legionella (Rhoads , 2016b). After 5-month acclimation at 39 °C, the control rig was maintained at 39 °C, while the experimental rig was increased from 39 to 42, 48, 51 and 58 °C at intervals of 2–3 months. Both rigs were fed with Blacksburg (VA, USA) tap water pre-treated with Pentek GAC-BB granular activated carbon cartridges (Pentair Residential Filtration, Milwaukee, WI, USA). Factors examined within each rig included: water heater temperature setting (°C), time (5, 8, 10, 13 and 15 months), sample location (influent, recirculating line and distal tap), system (control rig and experimental rig), pipe orientation (upward and downward flow), phase (biofilm and bulk water), and water use frequency (high, medium, low: 21, 3, 1 flush/week respectively) in triplicate distal pipes. Notably, distal taps achieved room temperature within 25 min after flushing, in contrast to the recirculating line, which was consistently maintained at the temperature setting (Rhoads ). If not specified, ‘temperature' denotes the water heater setting, while ‘effective' temperature indicates the average in the distal pipe during stagnation (that is, across the flush/stagnation cycle).

Sample collection and processing

Biofilm swab (~65 cm2) and first-flush bulk water (~500 ml) samples were collected after 2- to 3-month acclimation time at each condition, at the end of the stagnation cycle. For biofilm sampling, sterile cotton-tip applicators (Fisherbrand, distributed by Fisher Scientific Company LLC, Pittsburgh, PA, USA) were used to swab the interior surface of the influent, recirculating line, and ends of the distal tap pipes accessible from threaded union connections. Bulk water samples were concentrated onto 0.22 μm pore-size sterile mixed cellulose ester filters (Millipore, Billerica, MA, USA). Filters and swabs were transferred to lysis tubes and subject to DNA extraction using FastDNA Spin Kits (MP Biomedicals, Solon, OH, USA) according to manufacturer protocol. Illumina 16S rRNA amplicon sequencing was performed following the Earth Microbiome Project 16S rRNA Amplification Protocol Version 4_13 (http://www.earthmicrobiome.org/emp-standard-protocols/16s/). Minor changes were as previously described (Ji ), including using Molecular Biology Grade Water (Quality Biological, Gaithersburg, MD, USA) and pooling 200 ng of amplicon from each sample before clean-up using QIAquick PCR Purification Kit (QIAGEN, Valencia, CA, USA). Universal primers 515f and barcoded 806r were used for PCR amplification. Illumina 16S rRNA amplicon sequencing (paired-end 250 bp, MiSeq platform) was performed by the Virginia Bioinformatics Institute (Blacksburg, VA, USA).

Microbiome data analyses

Paired-end reads were stitched using PANDAseq (Masella ) and filtered based on quality score (>0.90, range 0–1) and sequence length (252–255 bp) following the EDAMAME course materials (https://edamame-course.github.io/, accessed 01/24/2016). A de novo Operational Taxonomic Unit (OTU) picking strategy was adopted in QIIME 1.8.0 (pick_de_novo_otus.py, Caporaso ) at similarity cutoff of 0.97 (uclust; Edgar, 2010) using Greengenes 16S rRNA gene database v13_8 (McDonald ). Mitochondrial and chlorophyll sequences were removed from both the OTU table and the phylogenetic tree before downstream analysis. After removal of singleton OTUs, the sequencing depth of the 427 samples ranged from 25 266 to 304 049 with a median value of 82 078, which translated into 62 494 OTUs in total. All samples were rarefied to the smallest sequencing depth (that is, 25 266) before downstream alpha and beta diversity analyses (rarefied 10 times and 100 times, respectively) in QIIME 1.8.0 (Caporaso ). Shannon index was selected for the alpha diversity matrix to evaluate both evenness and richness components, which emphasizes richness and rare OTUs relative to the Simpson Index. Weighted UniFrac distance matrix was selected for beta diversity to capture phylogenetic distance while taking relative abundance into consideration (Lozupone and Knight, 2005). Principal Coordinate Analysis was further applied to the beta diversity and visualized in EMPeror (Vázquez-Baeza ).

Metagenomic data was predicted from the OTU table generated by closed-reference OTU picking strategy (pick_closed_referfence_otus.py, uclust similarity cutoff at 0.97, Greengene v13_8) in QIIME 1.8.0 (Caporaso ) following PICRUSt 1.0.0 online protocol (Langille ). Bray-Curtis distance matrix was calculated from predicted functional profile in QIIME 1.8.0 (Caporaso ). As PICRUSt predictions are based on sequenced reference bacterial genomes, weighted Nearest Sequenced Taxon Index (weighted NSTI) scores were applied to measure average phylogenetic similarities among the OTUs within a given sample and the reference bacterial genomes, with OTU abundance taken into consideration. In addition, the ratios of PICRUSt-compatible OTUs and total sequences were calculated to estimate the predictable proportion of microbial community. All factors were treated as non-numeric variables, and the relative contribution of each factor on the microbial phylogenetic and functional dissimilarities across samples were determined through Permutational Multivariate Analysis of Variance (Adonis, Anderson, 2001) among a priori defined groups in R studio (permutation=999, {vegan} version 2.3-0, Oksanen ). Adonis is less sensitive to dispersion (that is, within-group variation) than alternative methods including ANOSIM (Oksanen ). A complementary multivariate homogeneity of groups dispersions (variances) analysis (betadisper; Anderson, 2006) was applied to account for confounding dispersion effect.

To explore the relationship between biofilm and bulk water microbiota, we first examined if both phases trended similarly in response to various factors through Mantel and Procrustes tests (QIIME 1.8.0, Caporaso ). Second, we explored to what extent biofilm and bulk water shared microbes via core OTU analysis on both sample groups and individual paired samples (169 pairs of bulk water- biofilm distal tap samples). Based on Levene's test of homogeneity of variances, one-way ANOVA with Tukey's honest significant differences (HSD) method for pair-wise comparisons, or Kruskal and Wallis test with Tukey and Kramer (Nemenyi) test was applied ({car} version 2.1-2, Fox and Weisberg, 2011; {stat} version 3.2-2; {PMCMR} version 4.1, Pohlert, 2014). A Venn diagram ({VennDiagram} version 1.6.16, Chen and Boutros, 2011) was constructed for the former. A complementary indicator microbes analysis ({indicspecies} version 1.7.5, De Caceres and Legendre, 2009) was performed to inform the identity of the shared/exclusive OTUs of interest with thresholds for both positive predictive value (A) and sensitivity (B) were set at 0.7 with a 95% confidence interval (P<0.05).

Associations between Legionella spp. and other microbial genera were characterized using Maximal Information-based Nonparametric Exploration (MINE) analysis ({minerva} version 1.4.1, Reshef ) with a false detection rate correction at 0.05 (Benjamini and Hochberg, 1995). Mitochondrial sequences belonging to Acanthamoeba spp. and Vermamoeba vermiformis (née Hartmanella vermiformis) were detected. The ‘relative abundance' of mitochondrial V. vermiformis was calculated as the ratio of the number of mitochondrial sequences to the number of total sequences (no organelle or singletons), which was later incorporated into the taxonomy table at the genus-level before MINE analysis.

Data sharing

All sequences and metadata for this research are retrievable from QIITA Study ID 10502.

Results

Drivers of phylogenetic and functional microbiome dissimilarity

Phase (biofilm versus bulk water) was the major driver of phylogenetic dissimilarity across all samples (R2Adonis=0.246, PAdonis=0.001), followed by water heater temperature setting (R2Adonis=0.172, PAdonis= 0.001), sample location (influent, recirculating line, and distal tap; R2Adonis=0.124, PAdonis=0.001), and time (R2Adonis=0.100, PAdonis=0.001; Table 1). Water use frequency (R2Adonis=0.039, PAdonis=0.001) and pipe orientation (R2Adonis=0.013, PAdonis=0.001) had less pronounced effects on phylogenetic dissimilarity, although their impacts were apparent under certain conditions (for example, biofilm at T=58 °C, where lower water use frequency seemed to enhance the impact of pipe orientation, Supplementary Figure 2). Notably, elevated temperature and time together conferred a greater magnitude of impact on the phylogenetic dissimilarity relative to that of constant temperature and time (control rig) in both recirculating line (experimental rig, R2Adonis=0.759, PAdonis= 0.001; control rig; R2Adonis=0.444, PAdonis= 0.001) and distal tap (experimental rig, R2Adonis=0.249, PAdonis= 0.001; control rig; R2Adonis=0.112, PAdonis= 0.001) microbiomes, with the most striking contrast between the two rigs at 15 months when the experimental rig was at 58 °C (Supplementary Figure 3, distal tap). A statistically significant temporal trend was observed in influent samples (9.6–23.7 °C; R2Adonis=0.160, PAdonis=0.679), which was weaker than that in the control rig (consistently at 40 °C; R2Adonis=0.444, PAdonis=0.001) and experimental rig (incrementally increased from 40 °C; R2Adonis=0.759, PAdonis=0.001) recirculating lines. Temperature appeared to have exerted an influence on the microbiome composition over and above that of succession with time, given that there was no significant difference in the phylogenetic composition between the two rigs during the initial 5-month baseline at 39 °C (PAdonis=0.158).

Table 1

Relative impact of various factors on phylogenetic and functional microbiome dissimilarity across samples (Adonis and Betadisp, permutation=999)

Factor	Phylogenetic			Functional			Samples
	R2Adonis	PAdonis	PBetadisp	R2Adonis	PAdonis	PBetadisp
Phase	0.246	0.001	0.001	0.087	0.001	0.001	All samples (n=427)
Sample location	0.124	0.001	0.001	0.055	0.001	0.003	All samples (n=427)
t	0.100	0.001	0.001	0.398	0.001	0.001	All samples (n=427)
T	0.172	0.001	0.003	0.194	0.001	0.003	All samples (n=427), influent was set at 20°C.
Water use frequency	0.039	0.001	0.039	0.021	0.005	0.040	Distal tap samples (n=348)
Pipe orientation	0.013	0.006	0.716	0.013	0.008	0.971	Distal tap samples (n=348)
Phase	0.311	0.001	0.001	0.104	0.001	0.001	Distal tap samples (n=348)
Rig	0.025	0.158	0.624	0.019	0.267	0.185	Distal tap at 5 months, 39°C (n = 70)
t (at 39 °C)	0.112	0.001	0.043	0.396	0.001	0.244	Control rig distal tap samples (n=175)
Phase	0.466	0.001	0.001	0.138	0.001	0.001	Control rig distal tap samples (n=175)
t (with elevated T)	0.249	0.001	0.015	0.513	0.001	0.398	Experimental rig distal tap samples (n=173)
Phase	0.256	0.001	0.001	0.084	0.001	0.001	Experimental rig distal tap samples (n=173)
Phase	0.126	0.001	0.098	0.104	0.003	0.069	Recirc. Samples (n=45)
Rig	0.124	0.001	0.003	0.072	0.016	0.962	Recirc. Samples (n=45)
t (at 39 °C)	0.444	0.001	0.730	0.681	0.001	0.519	Control rig recirc samples (n=27)
t (with elevated T)	0.759	0.001	0.687	0.608	0.001	0.807	Experimental rig recirc samples (n=28)
t	0.377	0.001	0.443	0.621	0.001	0.455	Influent (n=24)
Rig	0.041	0.001	0.039	0.012	0.011	0.053	Distal tap samples (n=348)
T (<51 °C or ⩾51 °C)	0.175	0.001	0.001	0.183	0.001	0.001	Experimental rig distal tap samples (n=173)

Both time (t) and water heater temperature setting (T) were treated as non-numeric factors. Statistical significance was set at 0.05: italic font indicates non-statistically significant.

When considering the microbiome composition from a functional standpoint, employing PICRUSt, there was no significant dissimilarity as a function of phase or sample location (Table 1). Instead, temperature and time were most influential, though the impact of elevated temperature on function was not as great as that which had been observed on phylogenetic profile (Table 1). For the influent water, the functional microbiome was found to be dissimilar with time (R2Adonis=0.621, PAdonis=0.001). Despite the low PICRUSt-compatible OTU ratio (16.2±6.6%), 81.3±6.3% sequences were PICRUSt-compatible across all samples (Supplementary Table 1). Although complex environmental samples pose a challenge to PICRUSt, accurate metagenome prediction was achieved in a prior study of soil, with NSTI values of 0.17±0.02 (mean±s.d.; Langille ). NSTI scores aid in evaluating PICRUST accuracy, with a negative association with accurate metagenome prediction. Relatively low weighted NSTI scores in the present study indicated that PICRUSt was appropriate (mean±s.d.=0.10±0.02, Supplementary Table 1).

Effect of combined factors on microbiome composition

Grouping samples by phase, sample location, temperature and time resulted in 46 sample groups, which shared the same top abundant phyla (that is, phyla with over 1% of total sequences across all samples), but differed in their relative abundances (Figure 1). For instance, Nitrospirae favored the biofilm phase and were enriched with time in the influent and control rig distal taps, but diminished in the experimental rig distal taps with elevated temperature. In the recirculating line, Nitrospirae remained at a comparatively low relative abundance, with no obvious temporal trend in either rig. Almost all Nitrospirae sequences were classified as Nitrospira spp., which are nitrite-oxidizing bacteria with some members recently reported to be ‘comammox' (that is, complete ammonia oxidizer; Daims ) and to have an optimal temperature range of 35–40 °C. In the recirculating line, higher temperatures (51 °C and 58 °C) were associated with smaller net decrease in ammonia concentration (Supplementary Table 2) and concurrently lower relative abundance of Nitrosomonas spp. in bulk water (Supplementary Table 3), which likely suggests inhibited ammonia oxidization given that influent ammonia level remained stable. Some other ammonia-oxidizing bacteria containing genera, including Nitrosospira (Burrell ), were not found in the recirculating lines.

Figure 1

Taxonomic profiles of sample groups at the phylum level: (a) influent, (b) recirculating line, (c) distal tap. Time is in month (mon) and temperature indicates water heater temperature setting. Top abundant phyla were defined as those that were >1% of total sequences across all samples combined (percentage shown in brackets); the remainder were combined and lumped in a category designated as ‘Other'. The three blank columns correspond to lost samples.

Biofilm versus bulk water: indicator OTUs

The indicator OTU concept provided a framework for examining the preference of individual microbes for the biofilm versus the bulk water. Ideally, an indicator OTU should meet two criteria: (1) 100% of the OTU sequences belong to the samples from the defined group (positive predictive value A=1.0) and (2) the OTU is detected in 100% of the samples from the defined group (sensitivity B=1.0; Cáceres and Legendre, 2009; Cáceres ). With the confidence interval set at 95%, 8038 and 10 424 indicator OTUs were identified for biofilm and bulk water phase, respectively. Currently there is no consensus on the cutoff values for the positive predictive value (A) and sensitivity (B), whereas in general higher values imply stronger indicator OTUs. Applying 0.7 as an arbitrary threshold for both A and B, we observed that biofilm harbored less indicator OTUs (147) than did bulk water (232). Notably, OTUs within the same genus did not universally serve as indicators for the same phase or sample group. For instance, bulk water phase indicators included 11/631 Mycobacterium spp., 34/635 Meiothermus spp. and 1/188 Legionella spp. OTUs. Similarly, different OTUs within the same genus sometimes indicated different phases. Out of 1056 OTUs belonging to Nitrospira spp., 11 served as indicators for the biofilm phase while 7 others were indicators for the bulk water phase. This suggests that it may be inappropriate to assume that members of the same genus uniformly serve as indicators in domestic hot water systems.

Bulk water and biofilm trended with each other in response to various factors including time and water heater temperature (Mantel: R2=0.608, P=0.001; Procrustes M2=0.698). In general, each bulk water and biofilm sample pair shared 8–24% OTUs, which comprised 68–95% of the total sequences between the pairs (Figure 2). Comparison between the two rigs suggested that elevated temperature (T⩾51 °C) induced a subtle increase in the proportion of shared OTUs between the biofilm and bulk water (Tukey HSD; adjusted P51–39°C=0.046, adjusted P58–39°C=0.005), although this trend was not significant when viewed as the proportion of shared sequences (adjusted P51–39°C=0.992, adjusted P58–39°C=0.258) (Supplementary Table 4). At T<51 °C, high water use frequency was associated with the fewest shared OTUs and sequences between paired bulk water and biofilm samples (Figure 2; adjusted P<0.001 for 21/week – 1/week and 3/week – 1/week, Supplementary Tables 5 and 6). Pipe orientation (upward and downward flow) had impact on the relationship between bulk water and biofilm in some situations. For instance, under high water use frequency (21/week), downward pipe harbored higher shared OTUs between phases, except at 10 and 15 months in the control rig (39 °C).

Figure 2

Proportion of OTUs and sequences shared in paired bulk water and biofilm samples at distal taps. Proportion was calculated as either number of shared OTUs divided by number of total OTUs between the paired samples (membership); or the number of sequences belonging to shared OTUs divided by the number of total sequences between the paired samples (relative abundance) across time (t in month, mon) and temperature (T in °C). Color codes indicate water use frequency and shape indicates pipe orientation.

Influence of recirculating line microbes on those at the distal taps: core OTUs

Beta diversity analysis indicated that the greatest contrast occurred between the influent and recirculating line, with the distal taps displaying noticeable overlap with both of these locations (Supplementary Figure 4). Of particular interest was the relative contribution of microbes from the influent and recirculating line to those observed at the tap, which can be conceptualized as survival of allochthonous microbes to the periodically heat-shocked distal tap environment. To facilitate this comparison, control rig samples were subject to core OTU analysis (that is, OTUs shared among all samples within the given group) corresponding to 4 groups including distal tap biofilm, recirculating line biofilm, distal tap bulk water, and recirculating line bulk water (Figure 3). There were 14, 102, 19 and 22 core OTUs confined to each of these four individual groups, respectively, with a total of 19 core OTUs common to all groups. The vast majority of distal tap bulk water core OTUs (82/106 OTUs) were shared with recirculating line core OTUs (both phases included). Of these, 34 OTUs were only shared with recirculating line bulk water and 2 were only shared with the recirculating line biofilm. Interestingly, 20 core OTUs that were shared between the two biofilm locations were also encompassed in recirculating line bulk water core OTUs. This suggests that significant migration occurs from the recirculating line to the tap via physical transport of bulk water. Biofilm at the tap is also a likely source of bulk water OTUs at the tap, with 36 core OTUs shared between these two phases and 2 of which confined to the tap.

Figure 3

Venn diagram of core OTUs among each of the recirculating line and distal tap phases in the control rig across all sampling dates.

Negative association between alpha diversity (Shannon Index) and temperature

Alpha diversity trends were similar for both the biofilm and bulk water phases in that the influent was characterized by the highest Shannon Index, followed by the control rig recirculating line, and the experimental rig recirculating line (Figure 4a). In the biofilm phase, the experimental rig trended with the control rig with time; whereas in the bulk water phase, the experimental rig experienced a more drastic decrease concurrent with the elevated temperature. Linear regression ({stats} version 3.2.2, R Studio) further confirmed a negative association between Shannon Index and ‘effective' temperature (that is, average temperature across the flush/stagnation cycle) in bulk water samples (Figure 4b). Despite this negative association, influent samples, with ‘effective' temperatures below 15 °C, had lower alpha diversity than the fitted model, which contributed to a peak alpha diversity region around 20 to 25 °C. Distal tap samples did not vary much from each other in terms of ‘effective' temperature (26.8±2.2 °C), but displayed a wide range of alpha diversities. Generally, higher temperature (T⩾51 °C) appeared to be conducive to lower alpha diversity.

Figure 4

(a) Alpha diversity (Shannon Index) across influent and recirculating lines in both rigs (note for experimental rig temperature was elevated with time); (b) negative association between alpha diversity and ‘effective' temperature in bulk water samples (gray line is linear correlation). Time is in month (mon) and ‘effective' temperature is in °C.

Patterns in Legionella occurrence and association with other microbes

The relative abundance of Legionella spp. decreased with elevated temperature (Figure 5a), whereas it increased in the bulk water phase of the control rig with time (Figure 5b). Nevertheless, at T=51 °C, 3 out of 5 low water use frequency (1/week) bulk water samples displayed markedly high relative abundance of Legionella spp., 2 of which were from upward oriented pipes. Rhoads also confirmed that this combination of 51 °C and low water use frequency uniquely selected for L. pneumophila at distal taps of the current system using quantitative PCR. Pipe orientation (upward and downward) alone did not show an obvious impact on enrichment of Legionella spp. (Supplementary Figure 5, paired upward/downward pipes).

Figure 5

(a) Relative abundance (RA) of Legionella spp. in all samples under different water heater temperature settings (T in °C); (b) relative abundance (RA) of Legionella spp. in influent and control rig samples across time (t in month, mon). Relative abundance was calculated using Illumina 16S rRNA gene amplicon sequencing. Water use frequency and phase are indicated by color and shape, respectively. All influent samples were assigned a uniform temperature of 20 °C to facilitate comparison.

Mitochondrial sequences belonging to Acanthamoeba and Vermamoeba (V. vermiformis), which are known to serve as hosts for replication and enhanced virulence of L. pneumophila and other pathogens, were detected in 4 and 393 out of 427 samples, respectively. MINE analysis identified a total of 143 statistically significant associations between Legionella spp. and other microbial taxa under a false detection rate of 0.05. The highest maximal information coefficient (MIC) was with mitochondrial V. vermiformis (MIC=0.44). Legionella spp. was also found to be positively associated with Mycobacterium spp., Meiothermus spp., and Xanthomonadaceae family (MIC=0.35, 0.41, and 0.51, respectively); whereas negatively associated with Nitrospira spp. (MIC=0.28).

Discussion

This study utilized pilot-scale, replicated pipe rigs to unravel the complexity inherent in-building plumbing systems and gain insight into the impact of water heater temperature and water use frequency on the microbiome composition at the tap. The findings have important implications for informing selection of water heater settings as well as water system design, which are currently the targets for new design and optimization strategies for conserving energy and water. Such designs could inadvertently shape the numbers and kinds of microbes entering the built environment via tap water.

A total of 62 494 OTUs were identified across this study, about a quarter of which were accounted for by influent samples (15 478 OTUs, both bulk water and biofilm). This number was also approximately one order of magnitude higher than the maximum number of OTUs reported in prior drinking water microbiota surveys (6488 in Lautenschlager ; Pinto ; El-Chakhtoura ; Lührig ; Proctor and Hammes, 2015), which have focused exclusively on cold water. We found that elevated water heater temperature drove significant change with respect to both the phylogenetic composition and predicted functions of the microbiota at the tap. Interestingly, the bulk water composition shifted with a threshold response around 51 °C, whereas the biofilm displayed more of a continuous succession pattern as the temperature was gradually elevated. The findings are consistent with the expectation of the environment selecting microbial populations (the niche theory, Hutchinson, 1957), which could theoretically provide the basis for a prebiotic framework for selection of ‘beneficial' microbes in the built environment. In interpreting the results throughout, it is important to bear in mind that the analysis applied in this study was entirely DNA-based, and therefore cannot directly distinguish live, dead or viable but non-culturable microbes. However, systematic comparison across various interconnected locations of the pipe rigs through various operating conditions and time provided an indication of which microbes were likely most adaptive to each location/condition.

Although there appeared to be a sharp bifurcation between biofilm and bulk water microbial composition, as reported by others using culture and 16S rRNA gene based clone libraries (Martiny ), highly similar functions between the two phases were predicted. This implies functional redundancy of the microbiota, or what has been dubbed as a ‘portfolio effect' in which positive and negative changes in individual taxa observed in the bulk water and biofilm essentially sum to zero (Allison and Martiny, 2008). Conceptually, the distinct dispersal limitations between planktonic (bulk water) and sessile (biofilm) states would allow for the development of functionally-equivalent, yet taxonomically diverse, species adapting to a similar set of environmental conditions (the neutral theory, Hubbell, 2001). In fact, neutral theory has been proposed in conceptualizing biofilms as microbial ‘landscapes' (Battin ). Rather than a niche-neutral dichotomy, recent studies have recognized the importance of both selective and neutral processes in shaping gastrointestinal (Jeraldo ), soil (Dumbrell ) and wastewater (Ofiţeru ) microbiomes.

Fundamental to this study is the fact that the distal tap environment is distinct from that of the recirculating line. Distal tap water samples essentially consisted of recirculating line inflow that had stagnated and cooled to room temperature (within ~25 min) since the previous flush of hot water. Stagnation has been reported to be a significant driving force for shift in potable cold water microbial communities (Lautenschlager ; Ji ). Within the control rig (T=39 °C), the substantial proportion of core OTUs unique to the distal tap bulk water (Figure 3) suggests that stagnation as short as 8 h (high water use frequency 21/week) is sufficient for a shift away from the core microbes defining the recirculating line bulk water or biofilm. It remains unclear, though, at what rate conversion from core to satellite microbiome generally occurs and whether elevated water heater temperature would influence the conversion time. The difference observed in as little as 8 h also demonstrated that the DNA-based approach applied in this study was able to capture microbial community shifts occurring on a short time scale.

The interrelationship between biofilm and bulk water microbes in building plumbing is one that is widely debated and of practical importance. For example, biofilm is widely considered to be a key environment for proliferation of OPs (Temmerman ), but release to bulk water and conveyance via the tap is generally necessary to cause disease (Steinert ). In this study, the high water use frequency was also characterized by the lowest proportion of shared OTUs and sequences between paired biofilm and bulk water sample (Figure 2, Supplementary Table 5). This suggests that the shorter stagnation time translates into less time for microbial interaction between phases to occur and take root. Microbes released from the biofilm to the bulk water would also be flushed out more frequently at a high use tap. Interestingly, the effect of higher water use frequency was diminished by elevated water heater temperature (T⩾51 °C, Supplementary Table 6). Elevated temperature (T⩾51 °C) was associated with a slightly increased proportion of shared OTUs between the two phases (Supplementary Table 4), which could be related to an overall decrease in total OTUs at the tap (Supplementary Figure 6). These findings have important implications for recommended sampling and monitoring protocols and whether biofilm or bulk water is most appropriately targeted over a range of conditions.

As noted above, the effect of gradually elevating the water heater temperature was not continuous, particularly in the bulk water. Notably, 51 °C appeared to be a critical threshold, with phylogenetic dissimilarity among distal tap biofilm samples, alpha diversity in distal tap water samples, ammonia oxidation in recirculating line bulk water, and interactive effects with water use frequency all bifurcating around this temperature. One possibility is that 51 °C demarcates a transition from a mesophilic to a thermophilic environment. For example, in anaerobic digestion, the mesophilic range is defined from 35–37 °C, while thermophilic digestion occurs from 50–57 °C, each imparting distinct dominant pathways for methanogenesis (Ali Shah ). In sum, microbes must be biochemically adaptive, such that proteins, nucleic acids, ribosomes, and lipids can maintain function at hot temperatures. An earlier study of 204 complete genomes and proteomes reported increased purine load index (A+G content) in nucleotide composition as well as increased sum of fractions of IIe, Val, Tyr, Trp, Arg, Gly and Leu amino acids under elevated optimal growth temperatures (−10 to 110 °C, Zeldovich ).

We considered whether the ‘effective' temperature profile (that is, the integral of instantaneous temperature over time between two sampling events) would provide a more meaningful predictor than the water heater temperature setting on observed phylogenetic dissimilarity patterns. Using this approach, distal tap sample compositions overlapped between those of the influent and recirculating samples (Supplementary Figure 4). Further, a negative association between ‘effective' temperature and microbial diversity (Shannon Index) was observed in distal tap and recirculating line bulk water samples. Although we could find no prior domestic hot water system studies for comparison, our observation is consistent with previous surveys of geothermal systems as a natural analog. A negative linear relationship between three temperature groups (62–87 °C) and the Shannon Index was reported for four sediment microbial communities in Great Boiling Spring (Cole ). Similar negative association were documented along a temperature gradient of 52–75 °C in a Japanese hot spring (Everroad ). Importantly, in the present study, the ‘peak' alpha diversity occurred around the ‘effective' temperature of 20–25 °C, which concurs with an extensive survey on 36 geothermal areas with 135 samples (Sharp ), where they observed a bell-shape alpha diversity curve over a temperature range of 7.5–99 °C with peak value at 24 °C. Notably, the indoor built environment is commonly maintained from 20 to 25 °C.

Given that it is impossible to eradicate microbes from water systems, a ‘probiotic', or more specifically a ‘prebiotic' framework has been proposed for desirably manipulating the microbiome at the tap (Wang ). Elevated temperature, in particular, can be highly effective for controlling OPs (Buse and Ashbolt, 2011; ANSI/ASHRAE, 2015). For instance, L. pneumophila has an optimal growth range of 25–42 °C (ASHRAE, 2000), with measurable deactivation at temperatures higher than 50 °C (WHO, 2007). In fact, a recent study has proposed a diagnostic approach to inform potential L. pneumophila risk through closely monitoring temperature profiles at critical control points within the hot water distribution system (Bédard ). On the other hand, elevated temperature is not universally effective, particularly when only temporarily applied. Legionella have been observed to proliferate via necrotrophic growth on heat-killed cells following thermal disinfection (Temmerman ). Thus, the present study makes it clear that temperature settings have much more broad sweeping effects on the microbial composition at the tap, likely exerting indirect microbial ecological controls on Legionella proliferation. Here we observed that water heater temperature settings from 39–58 °C were negatively associated with relative abundance of Legionella, but was less effective for Legionella control after cooling at the tap relative to the conditions of the consistently hot recirculating line. Further, at 51 °C, markedly high relative abundance of Legionella OTUs were observed in the bulk water phase at low-use distal taps, a measurement consistent with high L. pneumophila gene copy numbers determined by quantitative PCR (67 times higher on average than the inflow recirculating line water, Rhoads ). It is speculated that, at 51 °C, the heat shock served to diminish microbes competitive with Legionella (that is, ‘probiotics'), creating an ecological ‘sweet spot' where Legionella could recover and proliferate with the extended stagnation (Rhoads ). In fact, heat pre-treatment (50 °C for 30 min) is widely applied to increase the specificity for isolating Legionellae from hot water systems (Leoni and Legnani, 2001).

Mitochondrial Acanthamoeba spp. and Vermamoeba vermiformis sequences were detected, with the latter occurring in 92% of samples. Mitochondrial 16S-like ribosomal small subunit rRNA was reported to yield genotype-level identification and classification for Acanthamoeba comparable to that of 18S rRNA (Booton ). A recent study (Fučíková and Lahr, 2016) using the complete mitochondrial genomes of Acanthamoeba castellanii and Vermamoeba vermiformis reported higher resolution of intraspecific phylogenetic divergences relative to nuclear rDNA. However, little is known about the primer (515F/barcoded 806R) specificity towards amoebal mitochondrial sequences. A major obstacle is the availability of sequenced whole amoebal mitochondrial genomes, for example, 2 for V. vermiformis in GenBank (accessed 11/14/2016): KT185627.1 and GU828005.1 (identical to NC_013986.1). Each of the two V. vermiformis mitochondrial genome sequences contains only one copy of 16S-like ribosomal small subunit rRNA. Nevertheless, as number of mitochondria per V. vermiformis cell vary as a function of metabolic status, the ‘relative abundance' for mitochondrial V. vermiformis is examined as a relative indicator of V. vermiforms cells.

Associations with Legionella spp. were examined to gain insight into potential antagonistic or synergistic relationships. In particular, a strong positive association was noted between the relative abundances of Legionella spp. and mitochondrial V. vermiformis. This finding suggests the likely prevalent existence of V. vermiformis in both bulk water and biofilm phases of the rigs and potentially a positive association between V. vermiformis and Legionella spp. in terms of abundance. From a microbial ecology perspective, some free-living amoebae are either obligate or facultative pathogens (Brieland ; Lorenzo-Morales ) and can serve as vehicles/reservoirs for L. pneumophila (Wadowsky ; Kuiper ). Positive correlations between V. vermiformis and Legionella have been reported in treated and distributed waters (Valster ; Wang ), yet noted as conditional relationships and hypothesized to be dependent on strain type (Buse and Ashbolt, 2011) and in-building plumbing temperature (Rhoads ). Interestingly, Mycobacterium spp., a genus that contains some OP members, relative abundances were also found to be positively associated with those of Legionella spp. Previous studies have noted co-occurrence of M. avium complex and Legionella in water and/or biofilms in hot water systems ranging from 40 to 54 °C (Bukh and Roslev, 2014) as well as in a drinking water distribution system (Whiley ). Positive correlation between Legionella spp. and Mycobacterium spp. gene copies measured by quantitative PCR has also been reported in a drinking water distribution system (Lu ).

The design of the rig employed in this study uniquely served to simulate building plumbing conditions in a controlled, replicated fashion and isolate the effect of individual variables. The kind of insight gained would not be possible from a field study or simplistic bench-scale experiment alone, as the extreme complexities encountered from building to building make it impossible to identify causal factors. Nonetheless, it is important to note key differences in the design of the present study relative to the field. For instance, much longer distal pipes in the field would likely magnify impacts of cooling from water heater set points, pipe orientation and regrowth on levels of OPs in water. Real-world hot water systems will also have more complex water use patterns and employ a wide range of materials, which together will contribute to greater heterogeneity across distal tap microbiota. Overall this study demonstrates the possibility of implementing practical design and operation measures for intentionally controlling the microbiome composition at the tap, a topic of increasing interest not only for pathogen control, but also as linkages between the human microbiome and that of the built environment become more apparent.