3.1. Variation of water environmental factors during carcass decay
All water environmental properties including TN, TOC, NO3-N, NH4-N, DO, CON and pH were significantly affected by treatment (or corpse degradation) (Fig.S1; Table S2, One-way ANOVA, P < 0.01). Compared with control groups, pH, DO and NO3-N in experimental groups decreased significantly (Fig.S1-a, c, e, One-way ANOVA, P < 0.01). However, CON, NH4-N, TOC and TN in experimental groups increased significantly (Fig.S1-b, d, f, g, One-way ANOVA, P < 0.001). In addition, NH4-N, CON and pH in the experimental groups showed a growing trend during corpse decomposition (Fig.S1-a, b, d). On the contrary, DO, TOC, TN and NO3-N in experimental groups showed a downward trend (Fig.S1-c, e, f, g). And only NH4-N, CON and pH were significantly affected by time (Table S2, Two-way ANOVA, P < 0.001). The interaction between corpse degradation and time also had significant effects on the concentrations of CON, TOC and NH4-N (Table S2, One-way ANOVA, P < 0.01).
3.2. Diversity and abundance of ARGs during corpse decomposition
Among all 108 target genes, only 8 MGEs (including two integrase and six transposases) and 56 different ARGs were detected in all water samples. ARGs were classified into 9 types which included vancomycin, sulfonamide, beta-lactamase, aminoglycoside, chloramphenicol, Macrolide-Lincosamide-Streptogramin B resistance genes (MLSB), multidrug, tetracycline, and others (Fig.S3-a). And compared with control groups, the number of tetracycline, beta-lactamase, and vancomycin genes in carcass groups increased (Fig.S3-a). ARGs in all samples were divided into four resistance mechanisms which included cellular protection (43.87%), antibiotic deactivate (36.04%), efflux pump (16.25%), and others (3.84%) (Fig.S3-b).
Notably, the absolute abundance of all MGEs and ARGs was significantly impacted by corpse degradation (Table 1, Two-way ANOVA, P ≤ 0.001). Time or the interaction between treatment and time also posed an impact on absolute abundance of most resistomes, which included vancomycin, tetracycline, sulfonamide, multidrug, MLSB, chloramphenicol, and beta-lactamase (Table 1, P < 0.05, Two-way ANOVA). However, these two factors insignificantly affected on MGEs’ absolute abundance (P > 0.05, Two-way ANOVA). Multidrug resistance genes’ absolute abundance was mostly affected by treatment (i.e. corpse decomposition) (F = 154.190, P < 0.001), time (F = 21.787, P < 0.001) and their interaction (F = 21.785, P < 0.001). Similarly, treatment had a significant effect on the relative abundance of aminoglycoside, vancomycin, tetracycline, and sulfonamide (all P < 0.05, Two-way ANOVA). Time posed a significant impact on the relative abundance of aminoglycoside, beta-lactamase, multidrug, sulfonamide, tetracycline and vancomycin (P < 0.05, Two-way ANOVA). The interaction between treatment and time had a significant effect on the relative abundance of beta-lactamase, chloramphenicol, MLSB, tetracycline and vancomycin (P < 0.05, Two-way ANOVA). Compared with absolute abundance (PERMANOVA, R2 = 0.363, P = 0.001), cadaver degradation had a greater impact on the relative abundance (Table 2, PERMANOVA, R2 = 0.428, P = 0.001). Although time had no significant effect on the whole ARG profiles in all samples (Table 2, P > 0.05), time posed significant effect on the relative abundance of ARG profiles in experimental groups (PERMANOVA, R2 = 0.363, P = 0.029). In addition, time had significant effect on the absolute abundance of 24 ARGs genes (i.e. acrA-05, acrF and so on) and the relative abundance of 19 ARGs genes (such as aadA-01 and ampC-09) in experimental groups (Table S8, Table S9), while insignificant influenced most of ARGs’ abundance (37 ARGs’ relative abundance and 32 ARGs’ absolute abundance). In other words, the abundance of most ARGs in experimental samples were not affected by time and maintained relatively stable in different corpse decay stages.
Thirty-seven shared ARG subtypes among three different copse subgroups persistently existed throughout decomposition process, and the proportion of these ARGs accounting for 68.5% proportion of ARGs in all corpse groups (Fig.S2). Interestingly, the number of the shared ARGs subtypes between control groups and experimental groups was only 14, but later was 29 and 26 at 7th, 15th and 100th day, respectively (Fig.S2). In addition, all the ARG profiles were significantly influenced by treatment at each time point (Fig.1-a, b, Table S3, PERMANOVA, P < 0.05). And the effect of treatment increased firstly and decreased during decomposition, which reached the highest at 15th day (Table S3, PERMANOVA, absolute abundance R2 = 0.607, relative abundance R2 = 0.759, P < 0.05). In addition, no significant difference was found among all control groups at different time points (Table S3, PERMANOVA, P > 0.05), while significant difference was existed on 7KE vs. 15KE (Table S3, PERMANOVA, absolute abundance R2 = 0.255, relative abundance R2 = 0.403, P < 0.05). Notably, the ARG profiles of the three experimental groups became more similar or converged over degradation time (Fig.1-a, b).
Both the absolute and relative abundance of the total ARGs in experimental increased significantly compared with control groups (Fig.2-a, b). In detail, the abundance of total ARGs increased firstly, and then decreased over time, but did not recover to level of control group yet. Additionally, the relative abundance of ARGs was enriched about 2-266 folds in experimental groups compared to the control group (Fig.2-d) while the enrichment multiple of these ARGs’ absolute abundance was even 259-413,640 folds (Fig.2-c). In the experimental groups, some ARGs’ absolute abundance such as tetB-02, aadA-02 and tetA-02 were enriched over 258, 59,194 and 413,639 times than the control groups, respectively. Obviously, most of ARGs in experimental groups displayed higher concentrations compared to the control groups, and remain a relatively stable abundance over decomposition (Fig.S4-a, b). After screening for significantly enriched ARGs of experimental groups, the abundance of ARGs genes such as tetB-01 and tetA-02 were significantly enriched in the cadaver groups in the three time points (i.e. 7th day, 15th day and 100th day) (Fig.3, Fig.S5). In addition, the abundance of MGEs such as intI-1(clinic) and tnpA-01 significantly enriched in experimental groups and slightly increased over decomposition stages (Fig.3, Fig.S6).
It was also shown the different distribution of ARGs and MGEs in corpse groups at 7th, 15th and 100th day of corpse decomposition (Fig.4-a, b, Fig.S7). The abundance of vancomycin genes was enriched at 7th day, 100th day while the abundance of tetracycline genes was enriched at 15th day (Fig.4-a, b). Besides, efflux pump was main resistance mechanism at 7th day and 15th day while cellular protection was main resistance mechanism at 100th day (Fig.S7).
In order to investigate the assembly mechanisms of ARG communities of drinking water, we govern stochastic and deterministic processes’ relative importance during corpse decomposition (Fig.S3-c). Deterministic processes in corpse groups dominantly affected shape of ARG profiles at all the time points instead of stochastic processes (all compositional determinacy > 50%), indicating treatment (i.e. corpse decay) may more influence the shape of ARG profiles than stochastic factors.
3.3. Bacterial community composition and opportunistic pathogens during corpse decomposition
We acquired 1,864,462 high-quality sequences across 30 samples (mean 62,149 reads each sample, SD =19,819, min = 14,843, and max = 95,345), which were clustered into 11,053 unique OTUs with a threshold of 97% sequence similarity.
During the corpse decomposition, Proteobacteria and Bacteroidetes were main phyla (Fig.S8-a). Compared with control groups, the relative abundance of Bacteroidetes and Firmicutes increased significantly in corpse groups at the three time points (Fig.5-a, Fig S9-a, b, c, Table S4, P < 0.001, Two-way ANOVA). The Proteobacteria relative abundance in experimental groups gradually increased over time (Two-way ANOVA, P < 0.001). At genus level, Blvii28 and Novispirillum dominated in experimental groups (Fig.S8-b). Blvii28 and Novispirillum relative abundance increased significantly compared with control groups (Fig.5-b, Fig S9-a, b, c, Table S5, P < 0.001, Two-way ANOVA). In addition, Novispirillum abundance declined over time (P < 0.001, Two-way ANOVA). In addition, the dominant phylum exhibited different distributions at three time points (Fig.4-c). The dominant genera in experimental groups were Sphingomonadales (UG), Nevskia and Azospira at 7th day, 15th day and 100th day, respectively. At OTU level, OTU14020_Blvii28 (UG) and OTU20104_Novispirillum (UG) in corpse groups increased significantly compared with control groups (Fig.5-c, Fig.S9-d, e, f). And OTU14020_Blvii28 (UG) slightly increased firstly but then decreased over time. The dominant OTU in experimental groups at 7th day, 15th day and 100th day were OTU31089_Acetobacteraceae (UG), OTU3597_Azospira (UG) and OTU48291_mitochondria respectively (Fig.4-d).
Importantly, corpse decomposition had a significant effect on observed OTUs and chao1 index in three time points (P < 0.05, One-way ANOVA) with the exception for observed OTUs at the 15th day (P > 0.05, One-way ANOVA, Fig.S8-c, d). Both Chao 1 index and observed OTUs in corpse groups gradually decreased over time (Fig.S8-c, d). Additionally, based on the unweighted and weighted UniFrac distance, there were significant shifts of community structures and members among all groups over time (Fig.1-c, d). Corpse decay significantly influenced the bacterial communities’ structures and this effect became smaller over time (Table S5, PERMANOVA, weighted UniFrac at 7th day R2 = 0.935, 15th day R2 = 0.896, 100th day R2 = 0.682, P < 0.05). Notably, bacterial communities’ dissimilarity of experimental groups became gradually smaller over time.
In addition, we screened 20 opportunistic pathogenic genera according to previous references (Table S6), 12 genera i.e. Comamonas, Bacteroides, Pseudomonas, Brevundimonas, Delftia, Mycobacterium, Legionella, Halomonas, Clostridium, Burkholderia, Aeromonas and Acinetobacter were detected in our dataset. The relative abundance of Bacteroides and Pseudomonas in experimental groups increased significantly compared with control groups (Fig.5-d). In addition, Bacteroides and Aeromonas in experimental groups gradually decreased over time, and Aeromonas gradually decreased to the level of the control groups at 100th day (Fig.S8-e, f, P < 0.05).
3.4 The influencing factors of ARG profiles during corpse decomposition
For the relative and absolute abundance, ARG profiles were significantly impacted by treatment (i.e. corpse decomposition) (absolute abundance R2 = 0.363, relative abundance R2 = 0.428) and beta-diversity (absolute abundance R2 = 0.239, relative abundance R2 = 0.287, all P = 0.001) (Table 2). The relative abundance of ARG profiles (i.e. observed OTUs) was only significantly shaped by alpha diversity (R2 = 0.083, P < 0.01). In addition, seven environmental factors (i.e. TN, TOC, NO3-N, NH4-N, DO, CON, and pH) were also significant impacting factors shaping the ARG profiles (P < 0.01). In addition, we found ARGs’ absolute abundance were significantly related to some MGEs (i.e. tnpA-07, tnpA-05, tnpA-04, tnpA-02, tnpA-01, and intI-1(class1)) (P < 0.01). Meanwhile, ARGs’ relative abundance were significantly correlated to certain MGEs (i.e. tnpA-05, tnpA-01, intI-1(class1), and cIntI-1(class1)) (P < 0.05). Among these factors, we found the most affecting factors related to the absolute abundance were treatment (R2 = 0.363), DO (R2 = 0.359) and intI-1(class1) (R2 = 0.267), while the most influencing factors associated with the relative abundance were treatment (R2 = 0.428), DO (R2 = 0.425) and pH (R2 = 0.317).
And then, we divided these influencing factors into 5 parts (i.e. time, MGEs, environmental factors, microbial community and treatment) (Table S7). For the absolute abundance, treatment (R2 = 0.648), MGEs (R2 = 0.190) and environmental factors (R2 = 0.158) were significantly associated with ARGs (MRM, P < 0.01). And for the relative abundance, treatment (R2 = 0.544), environmental factors (R2 = 0.356) and microbial community (R2 = 0.015) were significantly linked with ARGs (MRM, P < 0.01). Notably, time had no significant correlations with both absolute and relative abundance of ARGs (MRM, P > 0.05). And compared with these important factors based on the PERMANOVA, we found these similar results based on the MRM analysis (Table S7). For the absolute abundance, treatment (R2 = 0.648), DO (R2 = 0.462) and intI-1(clinic) (R2 = 0.274) were three most important factors (MRM, all P = 0.001). And for the relative abundance, treatment (R2 = 0.544), DO (R2 = 0.449) and pH (R2 = 0.200) were three important factors (MRM, all P = 0.001).
3.5 Relationships among ARGs, water physicochemical factors, MGEs, and opportunistic pathogenic genera
Firstly, we detected correlations between MGEs and ARGs (Fig.7-a, Fig.S12, Fig.S13-a). The total ARGs’ absolute abundance was in positively relation to that of MGEs (R2 = 0.445, P < 0.001) while there was no association between the total relative abundance of MGEs and ARGs (Fig.S12, R2 = 0.007, P > 0.05). In addition, network analysis showed that tnpA-01 was positively associated with the absolute abundance of vanC-03 (r = 0.860), tetD-01 (r = 0.852) and mphA-01 (r = 0.874), and the relative abundance of dfrA12 (r = 0.762), tetL-02 (r = 0.710) and aadA-01 (r = 0.631). In addition, tnpA-07 was positively associated with the absolute abundance of cmlA1-01 (r = 0.821), tetE (r = 0.820) and aadA-01 (r = 0.767), and the relative abundance of dfrA12 (r = 0.703) and tetL-02 (r = 0.722) (Fig.S13-a).
In addition, we found correlations between environmental factors and ARGs (Fig.6-a, b, Fig.S10). For instance, DO was negatively associated with tetracycline (tetB-01 and tetG-02) of both the absolute and relative abundance at three time points (Fig.S9, P < 0.01). In addition, CON was positively correlated with aminoglycoside (aadA-01) and tetracycline (tetB-01 and tetA-02) of both the absolute and relative abundance (Fig.6-a, b, P < 0.01). And these correlations only existed at 7th and 15th day (P < 0.01), but were insignificant at 100th day (Fig.S10, P > 0.05). Similarly, network analysis showed the absolute abundance of aadA-01 (r = -0.854), tetB-01 (r = -0.789) and tetA-02 (r = -0.787), and the relative abundance of tetB-01 (r = -0.831), aadA-01 (r = -0.771) and dfrA12 (r = -0.740) were negatively associated with DO. CON was positively related to the absolute abundance of aadA-01 (r = 0.857), tetB-01 (r = 0.779) and tetA-02 (r = 0.819), and the relative abundance of tetB-01 (r = 0.865).
In addition, we also found the relative abundance of key genera (mean relative abundance > 1 %) was associated with the abundance of ARGs (Fig.6-c, d, Fig.7-b, Fig.S10, Fig.S13-b). Opportunistic pathogen like Comamonas was in positive correlation with absolute abundance of tetD-01 (r = 0.895), catB8 (r = 0.861) and floR (r = 0.853), and the relative abundance of tetD-01 (r = 0.897), dfrA12 (r = 0.881) and catB8 (r = 0.844) (Fig.S11). And Brevundimonas was positively associated with the absolute abundance of folA (r = 0.929), tetE (r = 0.865) and vanC-03 (r = 0.786), and the relative abundance of catB8 (r = 0.839) (Fig.S11).