Impact of storage condition upon gene targets quantitation

The stability of refrigerated influent wastewater samples was assessed at five time points (day 0, day 7, day 14, day 21, day 28), as it was determined that samples would not be retained longer than this in the refrigerator prior to processing. Indeed, in Figure S1(a-c), (i) no overall trend was observed in stability, (ii) the average standard deviation followed a decreasing order (CrAssphage > E-Sarbeco > N1), and (iii) non-detects did not affect most measurements. These results align well with other studies [24,25,26], as well as with identified good practice during the COVID-19 pandemic, and US NWSS 2024 guidelines [27], which recommended immediate sample processing, and if necessary to store samples at + 4 °C for up to 4 days [27,28,29]. Thus, storage of wastewater for subsequent processing for RT-qPCR analysis might be considered a viable solution and feasible for low-resource environments, avoiding any requirement to freeze samples.

The stability of frozen wastewater samples was assessed over a year (ten time points: day 0, day 1, 1 week, 2 weeks, 4 weeks, 8 weeks, 13 weeks, 26 weeks, 39 weeks, and 1 year) to consider the impact of longer-term storage upon late-processed and archived samples, common practices in WBE. Numerous non-detects and high variability were observed both within and between timepoints, with an absence of clear stability trends. This effect was particularly evident for E-Sarbeco and CrAssphage (observed in Figure S2(a-c)), as biological replicates for these targets showed higher prevalence of non-detects with few detected measurements, independent of the time spent frozen. Nevertheless, the variability at − 20 °C for all targets was greater than that observed following storage at 4 °C. Average standard deviation was higher for CrAssphage, followed by E-Sarbeco and N1. Since this behaviour was observed across all time points, it could be implied that it was not time in the freezer that caused degradation within individual biological replicates, but rather the freeze/thaw process itself.

On the basis of this stability study, wastewater sample storage at − 20 °C could not be considered the optimal choice due to the presumed deleterious effect on target recovery. This finding was in line with other studies [26, 28, 30] and with the detailed 2024 stability study from Williams et al., who found that SARS-CoV-2 N, E, and S genes degraded in wastewater upon freeze/thaw-cycling, showing reduced RNA integrity compared to samples stored at 4 °C [31]; on the other hand, these findings showed opposite results to Hokajarvi et al. which indicated that N2 and E-Sarbeco SARS-CoV-2 targets were stable in wastewater stored at − 20 °C for 58 days [32]. In addition, the stability of wastewater samples was complemented with up to four FT-cycles to further elucidate the stability under frozen condition. The data presented in Figure S3a shows that repeated freeze/thawing substantially reduces SARS-CoV-2 E-Sarbeco recovery which has a median result of a non-detect for all FT stability samples that had undergone at least 1 freeze/thaw sample. Conversely, the SARS-CoV-2 N1 target appeared more resilient to FT cycling in Figure S3b, despite the intra- and inter-timepoint variability observed, aligning with Fig. S2b and previous study by William et al. [31] and Thapar et al. [30]. The median value observed for CrAssphage upon freeze/thaw cycling never exceeds FT = 0 value, which may indicate FT instability. Therefore, FT cycling demonstrably had a negative impact on target stability and should be avoided if possible.

Since refrigerated conditions exhibited a reduced range of variability and lower occurrence of non-detects compared to frozen conditions, refrigeration was recommended as the optimal choice for RNA storage up to seven days for public health monitoring purposes. This recommendation was based on the observed absence of non-detects and lower variability, which was interpreted as an indicator of reliable and potentially safer storage within the conditions tested in this study. However, it was important to note that stability outcomes did not offer any evidence regarding the broader implications of storage conditions on surveillance and its potential consequences for public health. Specifically, this analysis did not account for long-term monitoring, intra- and inter-catchments variability (flow rates and population size) that may underline further implications which are discussed in the next section.

Epidemiological integrity assessment upon storage conditions

Although refrigeration seemed to be preferable to frozen storage as discussed in the previous section, the actual impact of sample storage conditions upon the generation of reliable epidemiological data for public health insights (e.g. early disease detection, assessment of surveillance infrastructures, and identification of policy impacts) required further investigation.

In particular, WBE targets a wide range of pathogens and genetic markers, and it is essential to understand how storage conditions influence the epidemiological interpretability of wastewater measurements and therefore impact the disease outbreak tracking across different targets. To address this, an epidemiological integrity assessment was conducted by comparing prevalence of SARS-CoV-2 in wastewater with the corresponding UKHSA-derived COVID-19 clinical cases for the WRC catchments on the day of sampling. Beyond evaluating SARS-CoV-2 specifically, this analysis aimed to assess how storage conditions may influence the epidemiological performance of wastewater surveillance more broadly, providing insights and an additional control applicable to other pathogen targets within WBE programmes.

Wastewater sample stability at 4 °C was analysed and sensitivity was determined as described in Data processing section. The sensitivity data are provided graphically in Fig. 3a–d below. Under refrigerated storage, the concentration of SARS-CoV-2 N1 and E-Sarbeco targets in wastewater exhibited a linear relationship with UKHSA-derived COVID-19 caseload on the day of samples collection, with Pearson’s r-values of 0.67 and 0.83 for SARS-CoV-2 E-Sarbeco and N1, respectively, increasing to 0.68 and 0.84 when stored no longer than 1 day before processing. Negligible divergence was observed among sites. This linearity indicated good discriminatory capabilities across the measured concentration range; moreover, the WBE methodology exhibited consistent 86% and 74–77% sensitivities for elevated COVID-19 case detection for N1 and E-Sarbeco, respectively, for all durations at 4 °C, using a specificity threshold of 95%. The lower sensitivity for E-Sarbeco aligned with its higher variability in the traditional stability study, so this may be related to the intrinsic characteristics of the gene target itself. Further, the data was well aligned between the four studied catchments indicating that the observed stability of both N1 and E-Sarbeco might reasonably be considered representative of that expected from samples originating from different WRC catchments. Moreover, these data did not indicate the existence of any significant stability issues for wastewater stored at 4 °C for 2+ days, as compared to 0–1 days, since differences remained within acceptable limits and retained sufficient discriminatory power and reliability. These results suggested that the method demonstrated increased accuracy due to a reduced false-omission rate. At the same time, the number of false negatives indicated room for improvement in sensitivity. This improved sensitivity with consistent specificity is critical, particularly in surveillance for low-incidence infectious disease outbreaks: in such an application, even a highly specific methodology may result in a much higher number of false positives than true positives due to the relative rarity of the monitored disease. Overall, the confusion matrix provided a nuanced picture of performance, supporting the method’s general reliability while highlighting the need to optimise sensitivity while safeguarding a high level of specificity. This approach demonstrates the importance of target-specific performance considerations when extending WBE methodologies to additional pathogen targets.

Fig. 3

a–d Scatter plot of SARS-CoV-2 E-Sarbeco (a, b) and N1 targets (c, d) at 4 °C (FT = 0) in wastewater samples for 0–1 day and 2+ days. Cgc (gc/L) is plotted against COVID-19 cpm for that day. Data points are grouped by WRC. In the top and bottom right corners, the performances of the model are reported

Conversely, samples stored at − 20 °C exhibited significantly reduced sensitivity and quantification of both SARS-CoV-2 N1 and E-Sarbeco: within the plots in Fig. 4a, b, this was observed in the clustering of measurements on the left side of the plots. Thus, the quantified gene target from a frozen sample may be expected to fall within a wide concentration range, providing no epidemiological information about the caseload from that day. As the actual caseload could not be distinguished from gene concentration alone, sample freezing eradicates the capacity to use SARS-CoV-2 quantitation to assess prevalence at community level and provide meaningful epidemiological insights (e.g. early detection of viral outbreaks). Moreover, a significant drop in method performance was observed when comparing refrigerated and frozen storage for both SARS-CoV-2 N1 and E-Sarbeco, with precision (the proportion of positive results that were not false positives) falling from 98 to 33% for N1. Based on these observations, storage under frozen conditions has a clear deleterious impact upon the quantification of SARS-CoV-2 RNA; these findings agreed with section ‘Impact of storage condition upon analyte quantitation’ and with the stability study conducted by Williams et al., who found that genetic material in wastewater samples from the N, E, and S genes of SARS-CoV-2 degraded upon freeze/thaw-cycling, showing reduced RNA integrity compared to samples stored at 4 °C [33]. Therefore, frozen samples within this dataset provided minimal information about the COVID-19 case rate. Although the study conducted by Hokajarvi et al. indicated that N2 and E-Sarbeco SARS-CoV-2 targets were stable in wastewater stored at − 20 °C for 58 days [32], it must be acknowledged that 86% of data included in the epidemiological integrity assessment described in this section were stored at − 20 °C for more than 58 days, with an average storage time of 268 days. The impact of storage time at − 20 °C for these (FT = 1) samples could not be determined since no correlation between gene quantitation and caseload could be identified for these frozen samples, regardless of storage time.

Fig. 4

a, b Scatter plot of SARS-CoV-2 N1 and E-Sarbeco recovered from frozen (FT = 1) wastewater samples plotted against COVID-19 cpm for that same day. Data points are grouped by WRC. In the top and bottom right corners, the performances of the model are reported

Characterisation of variability

Typically, the variability for a target within a methodology might be assessed through analysing the performance of non-native process controls [34, 35]. The size of this dataset and the breadth of replicate data enabled an alternative approach, however, whereby the variability was assessed for native gene targets. The approach used for the characterisation of variability within this manuscript was similar to that utilised by Wade et al. for the visualisation of SARS-CoV-2 N1 gene variability between replicates [36], albeit with three innovations: (i) the x-axis used the concentration instead of the CT value, (ii) the y-axis displayed the multiplicative difference between concentrations instead of the (analogous) additive difference between CT values, and (iii) the difference was here considered against the mean, rendering the resultant distribution symmetrical on the y-axis. These differences from Wade et al. allowed for the generation of a distribution of replicate values, and hence the direct determination of variability.

Analytical variability

Analytical variability (AV) is related to the quantification step via RT-qPCR and it was assessed through a comparison of the concentration of each RT-qPCR technical replicate (TR) to the concentration of its corresponding biological replicate (BR) to describe continuous variability through the determination of the confidence range. The confidence range associated with analytical variability for each target is reported in Table 3.

Table 3 Analytical variability data for each target. We report the confidence range and the equivalent coefficient of variation (%)

The %CV results were within the range reported by Taylor et al. for general RT-qPCR variability and below the upper bound of 200% [37]. Although high variabilities may be unavoidable, their impacts could be reduced through the application of a large number of technical replicates: [36] this may be prohibitively expensive, however. Variability typically decreased in the order CrAssphage > E-Sarbeco > N1 between the RT-qPCR targets, with discrete variability (prevalence non-detect %) and continuous variability (CR, %CV) following the same trend. It can be observed in Figures S4-6(c), that the probability of a technical replicate being a non-detect decreased with increasing biological replicate concentration: low sample concentration was a leading cause of non-detects, indicating that non-detects emerged from the same mechanisms that drove continuous variability, further indicating that non-detects were missing not at random, supporting the conclusions of McCall et al. [38] in their previous analysis of the causes for non-detects [38]. The spread of BR:TR deviation ratios in Figures S4-6(a) was largely concentration-independent: this indicated a consistent %CV throughout the range for each target, while the absolute variance differed with the magnitude of the biological replicate concentration. Furthermore, the spread of ratios in Figures S4-6(b and d) was found to approximate a normal distribution: from these values, the standard deviation (σ) and 95% confidence interval (μ ± 2σ) were estimated for all targets. All target datasets displayed in Figures S4-6 were long tailed compared to the normal distribution indicating that any given technical replicate may be significantly below/above the biological replicate concentration. This was likely primarily driven by stochastic variability of qPCR inhibition between technical replicates, in line with discussion of such variabilities by McCall et al. [38]. The long tails presented within all datasets reinforce the decision to use the non-detect-removed geometric mean for data processing within this methodology as an interim non-optimised compromise between retaining information and ensuring robustness to outliers. Such broad variability could not be entirely resolved at the data processing stage, but the use of the geometric mean with non-detect removal was an effective strategy to minimise the impact of outliers.

Processing variability

Processing variability (PV) describes the variability attributable to sample processing, including the variability from sampling, sample pre-treatment, sample concentration, and nucleic acid extraction. PV confidence ranges for gene targets included in this study are collated in Table 4.

Table 4 Processing variability data for each target, reporting confidence range, and the equivalent coefficient of variation (%)

PV yielded confidence ranges slightly lower than those from AV. It was observed that the probability of a biological replicate being non-detect (discrete variability) decreased with increasing DS replicate concentration, in line with Analytical Variability, although there was a relative reduction in the proportion of non-detect biological replicates as compared to technical replicates (Figs. S7-9(c)). This reflected the impact of non-detect removal, as only biological replicates comprising only non-detect technical replicates were defined as non-detects; therefore, the percentage of non-detects must decrease between technical replicates and the biological replicates they comprise. Moreover, the discrete PV typically decreased in the order CrAssphage > SARS-CoV-2 E-Sarbeco > SARS-CoV-2 N1 in line with the relative continuous variation for these three targets (represented by the confidence range and %CV), as described in Table 4 and similar to what observed for AV.

The relatively higher variability of E-Sarbeco as compared to N1 aligned with observations from other studies in the field. Zhang et al. noted that the SARS-CoV-2 E-Sarbeco target showed poor analytical sensitivity in wastewater, particularly as compared to SARS-CoV-2 N1, in contrast to the equivalent high performances of N1 and E-Sarbeco in spiked samples [17, 39]. It has been hypothesised that E gene variability in wastewater may be due to its variable loss in the clean-up procedure, potentially due to binding with viral proteins in the wastewater sample [40].

The high variability of CrAssphage made it unsuitable for use as an internal positive control and population biomarker. While AV for this target is in line with that observed for SARS-CoV-2, the PV was significantly higher. A confidence range of 8.32× indicated that a daily sample with a CrAssphage concentration of 1000 gc/L may exhibit biological replicate concentrations between 120 and 8320 gc/L: a range so wide as to be quantitatively inapplicable. In this context, it was unsurprising that WBE data normalisation with CrAssphage has been found to reduce the value of epidemiological information compared to un-normalised data [41, 42].

Using CrAssphage as a normalisation biomarker without a robust assessment of its variability within a catchment area (as conducted in this study) may introduce increased uncertainty in the resulting measurements, not due to the normalisation procedure per se, but from the inherent variability of CrAssphage [43, 44], affecting the ability to distinguish genuine detections from false negatives. This is crucial aspect to consider within outbreak detection monitoring programmes, where missed detections can lead to delayed responses and underestimation of public health risks; therefore the variability of normalisation approaches should be assessed to obtain a robust measurement to support timely and reliable decision-making.

PV was more likely caused by processing rather than sample collection; indeed, Ahmed et al. focused on assessing sampling variability for CrAssphage, HAdV, and PMMoV: on converting their data into confidence range terms, sample differences were only associated with approximately a 1.1× confidence range for CrAssphage and a 1.25× confidence range for HAdV on the most variable day [34]. Therefore, sample extraction and concentration were more likely the sources of variability: indeed, Wade et al. discussed the high variability of the concentration step, noting that recovery of φ6 phage (a concentration and extraction control) may range from 1 to 50% [36]. Assuming these recoveries were 5th and 95th percentile values, they were analogous to a confidence range of 7.1×, with a similar magnitude of variability as that observed for CrAssphage in this study [36], exceeding that observed for SARS-CoV-2 N1 and E-Sarbeco.

Methodological variability

Methodological variability (MV) describes the overall variability inherent in the method and could be considered as the combined expression of both AV and PV. As previously observed for both AV and PV, the spread of DS:TR concentration ratios in Figures S10-12(a) was largely concentration-independent (consistent %CV and confidence range); the spread of ratios in Figures S10-12 (b and d) approximated a normal distribution from these values; therefore, the standard deviation (σ), 95% confidence interval (μ ± 2σ), and confidence range (2σ) were estimated for all targets. These figures are reported in both Table 5 and Figures S10-12 (a, b, and d).

Table 5 Methodological variability data for each target, reporting confidence range, and the equivalent coefficient of variation (%)

In line with the AV and PV, discrete MV typically decreases in the order CrAssphage > E-Sarbeco > N1. This was broadly aligned with the relative MV for all RT-qPCR targets (represented by the confidence range and %CV), as described in Table 5. Within the total MV, a strong relationship was observed between gene concentration within a sample and the proportion of non-detects: DS with lower concentrations contained more TR and BR non-detects. As shown in Figures S10-12(c), the probability of a technical replicate being non-detect decreased with increasing DS replicate concentration. Although this was in line with both AV and PV, the proportion of non-detect replicates was higher in MV because this metric comprised both AV and PV, and thus the proportion was larger than either of these components individually. Within this analysis, non-detect observations did not seem to represent data missing at random, but rather data missing not at random dependent upon an unobserved gene expression within the daily sample itself [38]. Consequently, the removal of non-detects in the absence of corresponding partner replicates might lead to bias. The use of multiple TR and BR therefore provides a practical approach to partially mitigate bias caused by handling non-detects in RT-qPCR across different gene targets. Beyond our interim application of non-detect removal, further research is required for the development of optimal and pragmatic methodologies for non-detect handling for routine disease outbreak surveillance. The confidence range metric offered a more interpretable and coherent representation of variability for these data. The use of the confidence range was aligned with the best practice described by Taylor et al., wherein statistical analysis was carried out in ‘CT-space’ using normalised log-transformed data [37]. Although this approach was coherent and appropriate for data with symmetrical concentration-independent variability on a logarithmic scale, the use of the confidence range for RT-qPCR data in WBE represented a novel innovation, as most WBE studies rely on %CV for the characterisation of variability in RT-qPCR data. Alongside the confidence range, the percentage coefficient of variation (%CV) was reported for AV, PV, and MV; however, the %CV metric was not well suited to this data since the underlying variability is log-normally distributed, thus dispersion was most appropriately expressed as a multiplicative factor rather than an additive deviation around a mean (Supplementary Information). Since variability is multiplicative, a confidence range expressed on the log-scale could span several orders of magnitude. For example, 2σ equivalent to a 10× confidence range for analytical variability would imply that a TR measurement of 10 gc/µL would correspond to a BR concentration between 1 and 100 gc/µL with ca. 95% confidence, assuming no underlying bias. Under a linear-normal interpretation however, %CV values could readily exceed 100% [37]. The same 10× spread would translate to a %CV of ca. 250%, which misleadingly would suggest that values could fall below zero (e.g. between −40 and 60 gc/µL for a measurement of 10 gc/µL). In reality, under the log-normal model, the plausible range was strictly positive and would span 1–100 gc/µL. For this reason, %CV values, while provided for reference, would be less informative than the confidence range metric for multiplicative variability.

Assessing performance of composite metric for improved disease monitoring

SARS-CoV-2 N1 and E gene targets were both indicative of SARS-CoV-2 viral concentration. Therefore, the total SARS-CoV-2 viral concentration for each sample could be estimated using a weighted average of the available measurements. The authors were not aware of anyone else who had taken the step of merging WBE targets, except for Nagelkerke et al. who discussed a similar method for generating composite data in which data from consecutively analysed standard curves were merged to create a composite metric for the better estimation of SARS-CoV-2 concentration [45]. Nagelkerke found that this approach better mitigated stochastic error in RT-qPCR data for the proposed better prediction of COVID-19 caseload [45]; similar reduction in random error from the contributing gene targets was attempted in this work through the preparation of the Total Estimated SARS-CoV-2 metric. Final weights of 62.2:37.8 for SARS-CoV-2 N1 and E-Sarbeco were used respectively, derived from calculated %CV values for each target (Table 5).

Weights were defined as the inverse of the squared coefficient of variation ((1/%CV)2 = μ2/σ2) reducing bias towards lower concentrations, while giving a greater weight to the less variable N1 target as compared to E-Sarbeco. Though unusual, the use of %CV was in line with best practice as it was grounded in the statistics of the dataset [46]. However, two considerations influenced the weighting of N1 and E-Sarbeco: although the E gene exhibited greater variability and therefore warranted a lower weight, it was also observed at lower concentrations than N1 and would otherwise require scaling for comparability. The weighting approach prioritised adjusting for variability, although this may not be fully optimal given the opposing consideration of differential target expression rates.

The application of the composite metric Total Estimated SARS-CoV-2 reduced the incidence of non-detect BR as compared to E-Sarbeco and N1 (discrete PV). The number of non-detect DS decreased from 19 (N1) and 120 (E-Sarbeco) to 16 for Total Estimated SARS-CoV-2, out of a total of 525 DS (Figs. S13-15). Although this was not a large reduction from N1, including measurements from multiple targets improved robustness to non-detects. Furthermore, combining both N1 and E-Sarbeco TR in the calculation of BR of Total Estimated SARS-CoV-2 reduced the likelihood of non-detect BR within any given DS concentration range (Figures S13-15(c)); for example, in DS within the concentration range 1–10 gc/µL, the proportion of non-detects was reduced from 17.9% (N1) and 31.3% (E-Sarbeco) to 13.8%.

This showed that the composite metric Total Estimated SARS-CoV-2 preserved more information and robustness to non-detects than the individual gene targets alone, where non-detects were removed. However, using the composite metric led to an increase in the continuous processing variability (Figures S13-15a, b, d; Tables 6 and 7) as it was associated with the increase in TRs from one to two targets. As expected, the variability associated with an estimate of viral target concentration must be equal to or greater than the variability associated with the concentrations of its component gene targets. Thus, the increase in the DS:TR confidence range from 4.18× for SARS-CoV-2 N1 to 5.35× for Total Estimated SARS-CoV-2, in line with the confidence range of 5.15× for E-Sarbeco, could be seen as a consequence of the variability between the different gene target survival rates in wastewater, extraction rates in processing and detection rates in analysis. Figure S16 shows the differing performance of the N1 and E-Sarbeco targets at different concentrations. As such, the Total Estimated SARS-CoV-2 concentration for a BR would be between 125 and 165% of the measured value for an E-Sarbeco TR, while between 92 and 71% of the measured value for an N1 TR at low concentrations (1 gc/µL) and high concentrations (1 × 104 gc/µL) respectively.

Table 6 Processing, and analytical and methodological variability values for Total Estimated SARS-CoV-2 composite metricTable 7 Summary of continuous and discrete variability for all targets across AV, PV, and MV

The composite metric Total Estimated SARS-CoV-2 therefore offered mixed benefits when applied to BR. It reduced the proportion of BR (and consequently DS) observations classified as non-detects, while yielding greater continuous variability in DS concentrations, reflecting differences in performance of gene targets. Overall, this composite metric provided a more robust relative measure for daily SARS-CoV-2 outbreak tracking than individual targets, although with reduced precision. Importantly, the structure of this composite metric allows the incorporation of additional gene targets to better track disease dynamics, suggesting potential applications as flexible approach for different targets beyond SARS-CoV-2.

Limit of quantification (LOQ) and limit of detection (LOD)

The application of the ISO-11843 aligned approach to multiplicative variability is illustrated in Figure S17(a-b). Figure S17a shows the approach used within this manuscript visualised as the conventional additive formulation applied to log-transformed concentration data, where variability was treated as an additive normal error around zero (i.e. the LOD, 100 gc/µL) and the LOQ was defined as ‘LOD + 1.64σ’, consistent with constant variability on the log-scale across the measurement range. Figure S17b presents this same approach visualised as the analogous multiplicative formulation applied to the linear concentration data. Here, variability was expressed on the log10 scale and the LOQ could be understood as being obtained through multiplying the LOD by the 95th percentile of the methodological variability (MV95%). Figure S17b displays the same distribution as in Figure S17a but plotted against the linear concentration axis to show how the additive ‘LOD + 1.64 σ’ logic was preserved when variability is multiplicative instead of additive. LOD and LOQ values calculated using this approach are reported in Table 8, alongside comparison LOD′ and LOQ′ values determined additively using the linear-adjusted standard deviation derived from Eq. 5.

Table 8 Limit of quantification (LOQ) for extracted nucleic acid (extr. NA) and wastewater samples for each target

Values reported below the limit of quantification will typically originate from samples whose true concentrations lie within that range. However, assuming log-normal distributions, samples with true concentrations above the LOQ may be under-quantified into this range for < 5% of values, while samples with true concentrations below the LOD may be over-detected into this range for < 5% of values. Given the limited standardisation in WBE for the calculation of LOD/LOQ values from RT-qPCR data, LOD was defined as the non-detect cut-off (bottom calibrator) and the LOQMV95% values reported in Table