A universal empirical equation to estimate the abundance of carbapenem-resistant genes during aerobic digestion of wastewater sludge

Wastewater treatment produces sludge, which must be treated to reduce its pathogen and organic matter content before land application. Aerobic digestion is a common treatment method of sludge produced by wastewater treatment plants (WWTPs) up to a capacity of 2 m³/s, performed at mesophilic and thermophilic temperatures (Tchobanoglous et al. 2014). Antibiotic-resistant genes (ARGs), which are a part of DNA, encode resistance against bacterial infections (Davies & Davies 2010). Multiple studies demonstrated that antibiotic-resistant bacterial infections are becoming a major cause of death in North America and Europe (Bancroft 2007; Spellberg et al. 2011; Cassini et al. 2019). ARGs are an emerging group of contaminants in sludge (Mao et al. 2015). Aerobic digestion affects the quantity of ARGs, some of which survive aerobic digestion (Zhang et al. 2021). Therefore, safe design and monitoring of aerobic sludge digesters in terms of ARG removal would need engineers and operators to quantify ARGs.

Carbapenem-resistant genes (CRGs) are an important group of ARGs because carbapenem antibiotics are one of the few remaining options to treat bacterial infections (Sekyere 2016). There are many CRGs (Nordmann et al. 2012), among which bla_GES, bla_OXA and bla_IMP (Liang et al. 2018; Ekwanzala et al. 2020; Wang et al. 2022) have been detected in wastewater and sludge. Therefore, these CRGs are important and need to be the focus of further studies on wastewater sludge treatment. Nevertheless, it has not been investigated how the quantity of these genes changes during the aerobic digestion of wastewater sludge. Moreover, some bacteria are capable of horizontal gene transfer (HGT) using a gene that code for integrase class 1 (Int 1 gene). This gene can help to transmit CRGs (Xiong et al. 2013) and exists in wastewaters (An et al. 2018). Apart from that, 16S rRNA gene is a part of the bacterial genome, and it can be used for analyzing bacterial biomass (Smets et al. 2016). Hence, researchers usually study Int 1 and 16S rRNA variations with other genes in ARG studies (Nordmann et al. 2012; Wang et al. 2022).

To enumerate genes, DNA isolation and quantitative polymerase chain reaction (qPCR) are required (Xiong et al. 2013; Paulus et al. 2019), which are expensive, tedious and error-prone (Miłobedzka et al. 2022). In addition, practitioners working with wastewater sludge do not normally have the background and training for such molecular techniques. Thus, a reasonable method of gene estimation is needed that can be employed by the personnel who design and operate the aerobic digesters of wastewater sludge. Such estimation can be based on the sludge characteristics. Wastewater sludge consists of solids and water, where the former contains DNA. This compound can be chemically oxidized and consists of hydrogen, oxygen, carbon, nitrogen and phosphorus (Brezonik & Arnold 2011). Chemical oxidizability means that DNA is partly responsible for the chemical oxygen demand (COD) of the sludge. With this perspective, it is not surprising that the relationship between ARGs quantity and various sludge parameters has been scrutinized frequently (Tong et al. 2018; Tong et al. 2019a; Jang et al. 2020). The parameters include total solids, volatile solids (VSs) and volatile suspended solids (VSSs) (Jang et al. 2020; Poorasgari & Örmeci 2022a), phosphorus, COD and nitrogenous compounds (Tong et al. 2018, 2019a). Note that these parameters are normally measured for operation and monitoring purposes. VS and VSS encompass intercellular biomolecules such as proteins, lipids and carbohydrates. Also, the mentioned intercellular compounds are partly responsible for COD. A fraction of COD is soluble (SCOD), whose concentration changes during aerobic digestion (Kavitha et al. 2016), and a part of it is regarded as food for bacteria (Tchobanoglous et al. 2014). Ammonia is a nitrogenous compound in sludge whose concentration varies during aerobic digestion (Liu et al. 2012). This compound is a nutrient, but it becomes toxic when its concentration goes beyond certain limits (Gray 2004; Tong et al. 2019b). Associations between ammonia and ARGs have been demonstrated (Tong et al. 2019b; Goetsch et al. 2020). Soluble orthophosphate (SOP) is another nutrient (Tchobanoglous et al. 2014; Baird et al. 2017) whose concentration also changes during aerobic digestion (Poorasgari & Örmeci 2022c).

The associations described above imply that there is an ARGs–solids–nutrients (ASN) nexus in wastewater, which may be mathematically represented by regression modelling. A single regression model establishes a quantitative relationship between a dependent variable and a predictor, whereas a multiple regression model estimates the dependent variable from the values of more than one predictor. Several single and multiple regression models have been reported (Tong et al. 2018, 2019a, b). Nevertheless, the predictors of those models included, among other things, Chao and Shannon indices which are suitable for studying the biodiversity of a microbial community using specific software (Gihring et al. 2012; Thurkral 2017; Willis 2019). As well, quantification of particular parts of the 16S rRNA gene would be required to apply those models for estimating the genes (Větrovský & Baldrian 2013; Dueholm et al. 2022). This would make the models difficult to use by the wastewater practitioner. Based on the ASN nexus, Poorasgari & Örmeci (2022a) established a new modelling approach to predict the abundance of target genes during anaerobic sludge digestion. In that approach, all the predictors were parameters routinely measured during the operation of anaerobic digesters. However, it is not known whether the regression models developed by Poorasgari & Örmeci (2022a) would work for the aerobic digesters of wastewater sludge. Due to nitrification and denitrification during aerobic digestion (Tchobanoglous et al. 2014), nitrate concentration can change (Liu et al. 2012; Poorasgari & Örmeci 2022c). Therefore, factoring in the nitrate concentration may be important for building proper regression models of the aerobic sludge digesters.

In this investigation, it is hypothesized that the abundance of 16S rRNA, Int 1, bla_GES, bla_OXA and bla_IMP can be estimated from sludge total solids (TSs), SOP, SCOD, ammonia and nitrate concentrations. To test the hypothesis, the current study uses the ASN conceptual model to formulate a governing equation through different types of regressions. It is aimed to construct the equation in such a way that it is accurate and comprehensive, yet it can be understood and applied by the engineers and operators of aerobic sludge digesters. Such an equation is meant to estimate the abundance and removal of CRGs during the operation of aerobic digesters, while it could also be used for the safe design of such reactors.

Three batch aerobic sludge digesters were employed in the current investigation. The batch reactors were operated at 35, 45 and 55 °C, named AD.35, AD.45 and AD.55, respectively. The working volume of each reactor was 10 L. The reactor content consisted of 7.4 Kg primary sludge, 1.6 Kg T-WAS and 1 Kg inoculum. More information about the reactors is provided elsewhere (Poorasgari & Örmeci 2022b, 2022c). Sludge samples were taken from the digesters within time intervals of 5 days. The sludge samples were refrigerated at 1–3 °C until further analysis and DNA isolation.

DNA of the sludge samples was isolated using the Power Soil kit (Qiagen). About 0.25 g of sludge was used for each DNA isolation and three isolates were acquired for each sludge sample. A NucleoSpin gDNA clean-up kit (Macherey-Nagel GmbH & Co. KG) was used to improve the quality of the DNA isolates. After the clean-up, the DNA concentration of the isolates was measured using a NanoDrop 1000 spectrophotometer (NanoDrop, USA). Then, the isolates were stored at −20 °C until further analysis.

Five target genes were enumerated: 16S rRNA, Int 1, bla_GES, bla_OXA-48 and bla_IMP-27. Quantification was conducted using qPCR, where a Bio-Rad thermocycler (CC013200) and CFX Mastero software (Bio-Rad) were utilized. The reagent was SsoFast (Bio-Rad), whose ingredients were di-nucleotides, SYBR Green label, DNA polymerase enzyme and proprietary components. The volume of the reaction mix in each PCR tube was 25 μL. The mass of DNA for DNA isolates obtained from the sludge samples was 30 ng DNA per tube. The assay began with preheating at 95 °C for 180 s, and each cycle included denaturation at 95 °C for 10 s, annealing at primer-specific temperatures for 20 s, and polymerisation at 72 °C (see Table 1 for further information). The total cycle number of the CRGs was three more than that of the other two targets because preliminary tests had shown that the assays for the CRGs required more cycles to become quantifiable. Quantities with quantification cycle numbers more than those shown in Table 1 were considered zero.

The amplicon size and the specificity of primer annealing were first verified with the Primer-BLAST tool of the National Centre for Biotechnology Information (https://www.ncbi.nlm.nih.gov/tools/primer-blast/). To further confirm the annealing specificity of the primers and the size of the amplicons, gel electrophoresis (GE) was conducted (see Supplementary material S1 for details). Manufactured positive standards (Thermo Fisher) were used to prepare serial dilutions for calibrating the PCR machine. Standard solutions spanned 10²–10¹¹ copy numbers with an order of magnitude difference between each two subsequent dilutions. Calibration equations (Supplementary material, Table S2) were used to determine the absolute abundance of each target gene in the DNA isolate of the sludge samples in copy number per ng DNA. To determine the relative abundance of Int 1 genes and the CRGs, their absolute abundances were divided by the absolute abundance of their respective 16S rRNA genes.

The predictors of regression modelling included TS, SCOD, SOP, FAN and NO₃-N. SOP, FAN, NO₃-N and compounds responsible for SCOD constitute total dissolved solids (TDS), which is a part of TS. Nevertheless, our solids measurements made it clear that TDS accounted for less than 20% of TS. It was very likely that the remaining 80% contained a large portion of the total DNA, and hence, the target genes. That is why TS was considered a predictor, despite being partly represented by the other four predictors. The overall trends of FAN concentration (Supplementary material, Figure S3) were more meaningful than those of TAN (Poorasgari & Örmeci 2022c). This is the reason for choosing FAN for the ammonia-gene link in the model.

To determine the significance of the differences among the time series data, the two-factor analysis of variance (ANOVA) was conducted. Excel software was used to perform the ANOVA and regression modelling. The significance level of all statistical calculations of this study was 0.05.

Similar trends of the absolute abundance of the target genes (Figures 1(a)–1(f)) suggested the existence of correlations between the CRGs. Linear correlations between the target genes were established, and their R² and P-values are shown in Table 2. All the correlations were statistically significant. The overall strength of the correlations in terms of R² can be put in the following order: AD.55 > AD.45 > AD.35. With these correlations, linear equations can be readily built. Such equations are helpful for the operation of sludge digesters because they enable us to estimate the abundance and removal of all target genes by quantifying only one of the target genes. Correlations between 16S rRNA, ARGs and Int 1 have been reported for influents of WWTPs (Chen & Zhang 2013). In the current investigation, the high abundance of Int 1 and its correlation with bla_GES and bla_OXA-48 imply that sludge bacteria were capable of horizontal transfer. Jang et al. (2019) analyzed the correlations between Int 1 and multiple ARGs by the Pearson correlation coefficient (r or square root of R²). The reported r values were between 0.021 and 0.983, which translates to an R² range of 4 × 10⁻⁴ − 0.966. Also, the same study reported both positive and negative correlations. In the present research, R² values were greater than 0.8, and no negative correlation was observed. The R² values obtained from correlating Int 1 and the other target genes in the current investigation are more in line with the range 0.37–0.99 reported by Jang et al. (2020). The discrepancies could be due to the differences in the process, reactor feed and types of correlated ARGs.

Table 2

R² and P-values (subscripts) for linear inter-gene correlations

View Large

Table 3 displays kinetic parameters and respective coefficients of determination. The rate constants were specific for each target gene. The magnitude of the rate constants followed the order bla_OXA-48 > bla_GES > 16S rRNA. In addition, there was an overall increase in the absolute values of K₁ with increasing temperature. The R² values of the first-order fit for Int 1 in AD.35, AD.45 and AD.55 were 0.212, 0.368 and 0.870, respectively. This range is considerably higher than the range 0.01–0.54 reported by Diehl & LaPara (2010), where wastewater solids were digested at mesophilic and thermophilic temperatures. Notwithstanding, the R² value for the first-order fit of Int 1 in AD.35 was 2.55 times smaller than the corresponding value obtained by Diehl & LaPara (2010) from an aerobic digester operated at 37 °C. On the other hand, Jang et al. (2019) fitted the first-order kinetic equation smoothly to the abundance of Int 1 measured during aerobic digestion at 55 °C, similar to the current investigation. Second-order and CS kinetic fittings were also carried out. The second-order kinetic equation has previously been used to describe temporal variations of the target genes during physical and chemical processes only (Yu et al. 2017; Gao et al. 2020). The current study, however, shows that the second-order kinetics represented a variation rate of Int 1 and bla_OXA-48 in AD.55 more acceptably, i.e., higher R² values than the first-order kinetics did. CS kinetics could not be performed for bla_OXA-48 due to a lack of sufficient data. In the case of the other targets, the CS model yielded more accurate fits than the second-order model.

Table 3

Kinetic parameters and R² for kinetics of the target genes

View Large

It is important to note that the type of kinetics can change the design equations of continuous flow reactors including continuous flow aerobic digesters. Equations for the completely stirred mixed flow reactor (CMFR) and the plug-flow reactor (PFR) are provided below, where τ is SRT. Mathematical operations and assumptions used for deriving Equations (9)–(14) are available in Supplementary material S4.

The predictor ratios of the regression models were obtained from the concentrations of TS, SCOD, SOP, TAN and NO₃-N, which have been provided elsewhere (Poorasgari & Örmeci 2022b, c). The results of L_S and NL_S modelling have been summarised in Supplementary material, Table S5. R² values of up to 0.658 were observed for the L_S models. Also, the SOP concentration non-linearly correlated with bla_GES in AD.55 with an R² value of 0.999. Association has been observed between ARGs and phosphate (Ping et al. 2022). Initial nitrate concentrations were almost zero just before starting the aerobic digestion experiment. Nitrate concentration increased with time, as the reactor was being constantly aerated (Poorasgari & Örmeci 2022c). Chemical oxidation of ammonia requires specific conditions such as catalysts (Johnston et al. 2022). Therefore, the increase in nitrate concentration would have to be due to biological nitrification (Tchobanoglous et al. 2014). After 25 days of operation, there was a significant drop of nitrate concentration in AD.45. This suggests that there was denitrification, which was likely to occur because of the low concentration of dissolved oxygen (Poorasgari 2023) and probable localized anoxic islands deep inside the sludge aggregates. Thus, the observed increase and decrease in nitrate concentration (Poorasgari & Örmeci 2022c) signified the occurrence of nitrification and denitrification, respectively (Tchobanoglous et al. 2014). Bacteria involved in the biological conversion of inorganic nitrogen species may be able to carry ARGs (León & Roy 2003; Agersø & Sandvang 2005; Zhang et al. 2009). In the current investigation, there were strong negative correlations between the nitrate concentration and the target genes (R² = 0.81–1) in a non-linear manner. This is interesting because non-linear negative correlations have also been observed between ARGs and denitrification genes (Sun et al. 2017).

Table 4 shows further details of calculations for L_M and NL_M. The coefficients and exponents of the obtained models were different for each case. An explanation would be the dissimilar interactions between parameters affecting the abundance of genes in various kinds of sludge bacteria during aerobic digestion. Int 1 and CRGs could be nonexistent in some bacteria. Also, not all bacteria have the same copy number of the 16S rRNA (Větrovský & Baldrian 2013). In addition, the temperature at which anaerobic digesters operate impacts the composition of the microbial population (Gagliano et al. 2015). Similar temperature effects could be expected during aerobic digestion. Consumption of substrates and secretion of extracellular material may differ among microbial communities. As expected, the coefficients of the predictors in the L_M models were different from those in the NL_M models. The highest and lowest coefficients were 6.753 for the SOP ratio and −5.316 for the TS ratio, respectively. Both mentioned values belonged to the NL_M model for 16S rRNA ratios during digestion in AD.45. The lowest exponent was −8.262, which belonged to the TS ratio of the NL_M model for bla_OXA-48 ratios during digestion in AD.55. The weight of the parameters can be assessed by the absolute values of the coefficients and the exponents. Overall, the lowest absolute values of the coefficients and the exponents occurred in the FAN ratio in both L_M and NL_M models. This suggests that FAN concentration was not as important as the other predictors. FAN can be toxic at concentrations above 150 mg N/L (Gray 2004). FAN concentrations of AD.45 and AD.55 remained below the mentioned level. However, AD.35 exceeded the toxic level on day 35 (Supplementary material, Figure S3). Thus, FAN should have contributed as a nutrient to the abundance of the target genes in AD.45 and AD.55 during the whole digestion process. In AD.35, on the other hand, FAN may have acted as an inhibitor on the abundance of the target genes. Notwithstanding, the negligible and, in some cases, no contribution of FAN to the target gene ratio suggests that ammonia could have behaved as both nutrient and inhibitor in AD.35 and 45. That is, the antagonistic effects of nourishment and inhibition would have resulted in the diminished impact of the FAN ratio.

The absolute values of β for the NL_M models were smaller than those for the L_M models. Since the intercept is meant to correct the constant error of the model, smaller absolute values of β are more acceptable. In other words, estimations of the model are more reliable when the absolute value of β is smaller. Given the maximum value of 1 for the target gene ratio, an intercept with the absolute value equal to or greater than 1 means that the error compensation largely contributed to the estimated target gene ratio. An extreme case was the L_M model for bla_OXA-48 ratios of AD.55, where β was 2.055. Evidently, applying the exponents of the NL_S models (Supplementary material, Table S5) to Equation (7) improved the reliability of the model by reducing the absolute value of β. The intercept of the NL_M model went down to zero for bla_OXA-48 of AD.55. As well, the exponents’ applications reduced the SSE of multiple regression modelling for AD.35 and AD.45. The most noticeable one was observed in the case of the 16S rRNA of AD.45, where the SSE of the NL_M model was 0.053 below that of its corresponding L_M model. The overall R² and values of the NL_M models of AD.35 and AD.45 were larger than those of the respective L_M models. The most noticeable case was the 16S rRNA of AD.45. In that case, R² and of the NL_M model were 1.07 and 1.28 times higher than those parameters in the L_M model. The multiple regression models of AD.55 had the highest R² and values.

The NL_M models of AD.35 and 45 had lower overall P-values than the respective L_M models. The most striking of all occurred for bla_OXA-48 of AD.45, where the P-value of the NL_M model was 135 times smaller than that of the corresponding L_M model. Also, a 106-fold decrease was observed for bla_GES of the same reactor. The L_M models developed for 16S rRNA, bla_GES and bla_OXA-48 of AD.35, and for all target genes of AD.45 were insignificant (P-value > 0.05). On the other hand, the P-value was less than 0.05 for all multiple regression models of AD.55, and for all NL_M models of AD.35 and AD.45. The good fit at 55 °C is, at least partly, due to the fact that the number of observations and the number of model predictors were the same (Izenmann 2008). Therefore, the L_M and NL_M models for AD.55 need to be tested with other datasets with more observations, before they are considered by engineers in charge of aerobic digesters. Compared to the multiple regression models reported (Tong et al. 2018; Tong et al. 2019a; Jang et al. 2020; Poorasgari & Örmeci 2022a), the NL_M models of the current investigation resulted in a lower range for the P-value and a higher range for the R² value. Nevertheless, the cited investigations developed those models for processes other than the aerobic digestion of wastewater sludge.

It would be interesting to see how the goodness parameters of single and multiple regression models compare. To do so, the SSE values of the multiple regression models were divided by the minimum SSE values of their respective linear regression models, and the values of the multiple regression models were divided by the maximum R² values of their respective single regression models. The results of the former and latter divisions are termed the SSE ratio and the ratio, respectively. The multiple regression model is more accurate than the single regression model when the SSE ratio < 1 and the ratio > 1. On the other hand, the single regression model is more accurate when the SSE ratio > 1 and the R² ratio < 1. The computed SSE and ratios are presented in Table 5. All target genes and reactors had SSE ratios smaller than 1 and zero in several cases. Also, the ratio of AD.35 and AD.45 for all target genes was greater than 1. In the case of AD.35 for bla_GES, the ratio of the L_M model (Table 4) was 1.95 times larger than the highest R² value of the respective L_S model (Supplementary material, Table S5). Note that these calculations were carried out using all the decimals produced by Excel for the goodness parameters of the involved models but the values in Tables 4 and Supplementary material, Table S5 are shown with three decimals. If the R² values of the single (Supplementary material, Table S5) and multiple (Table 4) regression models are compared, the multiple regression models seem even more precise and accurate. From this perspective, multiple regression modelling with the general formula shown in Equation (7) was a stronger method of estimating the target gene ratio than single regression modelling. In the current study, the superiority of non-linear models to linear models is obvious. This may be explained by the effects of toxicity and biokinetics on the relationship between ARG abundances and model predictors (see Section 2.6). The effects of non-linearity on the decisions made by the designers and operators of aerobic sludge digesters would depend on the predictor ratios and the exponent of the predictor ratios. A worked example of using the non-linear form of Equation (7) is provided in Supplementary material S6.

Table 5

Comparison between single (L_S and NL_S) and multiple (L_M and NL_M) regression models

View Large

The single predictor in the new model, indeed, represents all predictors of the NL_M model (Equation (7)). The unified model was applied in both linear (L_S,U) and non-linear (NL_S,U) fashions, in accordance with the explanations in Section 2.6. Basically, the unified model would be more practical for engineering purposes, as its equation is simpler than Equation (7). The unified model was first tested with the experimental data obtained from anaerobic sludge digesters operated at the same temperature and with the same feed. The test was relatively successful as the experimental data fitted the unified models well (Poorasgari 2023). However, the unified single regression performed poorly for the observations of the current study conducted on aerobic digestion. This was quite obvious in the case of the reactor AD.35 with both L_S,U and NL_S,U models (Supplementary material, Figure S7 and Table S8). Finally, parameters such as pH could be incorporated into the regression models. However, including additional parameters would further complicate the models and veer them off the ASN conceptual model.

Experiments and data analyses were conducted to mathematically describe the ASN nexus during the aerobic digestion of wastewater sludge. The aerobic digestion process reduced the absolute abundance of the target genes quantified in the current study. The absolute value of K₁ followed the order bla_OXA-48 > bla_GES > 16S rRNA. Single and multiple regression models were established to predict CRGs as a function of TS, SCOD, SOP, FAN and NO₃-N. Calculations of the SSE and R² values indicated that the multiple regression models, including L_M and NL_M, were more accurate than the single regression models. In multiple regression models, FAN was the lowest impact predictor among the other four. The absolute value of β ranged 1.023–2.055 and 0–0.649 for the L_M and NL_M models, respectively. Consistent with this, the highest SSE value of NL_M modelling was up to 18.7 times, smaller than that of L_M modelling. Also, the NL_M models of AD.35 and AD.45 were superior to their L_M counterparts in terms of R², and P-values. This study, for the first time, demonstrated that the non-linear multiple regression models can estimate the abundance of target genes during the aerobic digestion of sludge. These models can help engineers and operators determine the SRT and digestion temperature required for removing the target genes of interest at a given rate during the aerobic digestion of wastewater sludge.

This research was funded by the Natural Sciences and Engineering Research Council of Canada, Discovery Grant program and Jarislowsky Chair Funds awarded to Professor Banu Örmeci.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.