Characterization of the Spatial Variation of Microbial Communities in a Decentralized Subtropical Wastewater Treatment Plant Using Passive Methods

: Septic tanks (STs), up-ﬂow anaerobic ﬁlters (UAFs), and horizontal-ﬂow constructed wetlands (HFCWs) are cost-effective wastewater treatment technologies especially efﬁcient in tropical and sub-tropical regions. In this study, the bacterial communities within a decentralized wastewater treatment plant (WWTP) comprising a ST, a UAF, and a HFCW were analyzed using high-throughput sequencing of the V3–V4 region of the 16S rRNA gene. Bacterial diversity and its spatial variation were analyzed at the phylum and family level, and principal component analysis (PCA) was applied to nitrogen- and organic-matter-degrading families. The highest percentage of nitrogen removal was seen in the HFCW (28% of total Kjeldahl nitrogen, TKN, and 31% of NH3-N), and our results suggest that families such as Rhodocyclaceae (denitrifying bacteria), Nitrospiraceae (nitrifying bacteria), and Rhodospirillaceae (sulfur-oxidizing bacteria) contribute to such removal. The highest percentage of organic matter removal was seen in the UAF unit (40% of biological oxygen demand, BOD5, and 37% of chemical oxygen demand, COD), where organic-matter-degrading bacteria such as the Ruminococcaceae, Clostridiaceae, Lachnospiraceae, and Syntrophaceae families were identiﬁed. Redundancy analysis demonstrated that bacterial communities in the HFCW were more tolerant to physicochemical changes, while those in the ST and the UAF were highly inﬂuenced by dissolved oxygen and temperature. Also, pollutant removal pathways carried out by speciﬁc bacterial families and microbial interactions were elucidated. This study provides a detailed description of the bacterial communities present in a decentralized WWTP located in a subtropical region.


Introduction
According to United Nations Sustainable Development Goal Number Six, a substantial increase in water treatment and reuse must be accomplished by 2030 to significantly reduce water scarcity worldwide and to protect the natural environment [1]. Centralized wastewater treatment plants (WWTPs) involve higher maintenance and operational costs than decentralized WWTPs, which are also easier to operate [2]. Anaerobic bioreactors (ARs) and constructed wetlands (CWs) are treatment stages that are commonly used in decentralized systems [3]; these units are classified as passive technologies, since they require low energy consumption and maintenance and operational costs [3,4]. Up-flow The objective of this study was to characterize the spatial variation in microbial communities between and within the treatment stages of a decentralized WWTP combining a septic tank (ST), an up-flow anaerobic filter (UAF), and a horizontal-flow constructed wetland (HFCW). Additionally, we investigated the influence of physicochemical parameters on the microbial structure within all three treatment stages. As a contribution to the literature, this work is focused on the study of the microbial communities throughout a multi-stage decentralized treatment system composed of a ST, a UAF, and an HFCW that uses passive methods and has demonstrated to treat domestic wastewater efficiently. This study seeks to gain a greater understanding of the microbial structure and behavior in this WWTP located in a subtropical climate and the microbial response to physicochemical variations within the system. A thorough knowledge of the structure and behavior of microbial communities within the system under study contributes to the development of strategies to enhance microbial removal processes and improve the performance of decentralized WWTPs integrated with anaerobic reactors and constructed wetlands in tropical and subtropical countries. Furthermore, this is one of very few studies to provide a detailed description of the bacterial communities in decentralized treatment plants located in subtropical countries (such as Mexico), where climatic conditions play an essential role in the performance of treatment systems.

Site Description
The study site is a decentralized WWTP that combines an ST, a UAF, and an HFCW to treat high strength domestic wastewater generated at a public R&D center located in the municipality of Zapopan (Jalisco, Mexico). The study site is located in a subtropical region with an average annual temperature of 21.8 • C and annual precipitation of 926.4 mm [30]. The WWTP receives a wastewater load of approximately 7.5 m 3 /day. Untreated wastewater is received in the pump sump, where it is pumped to the ST. At this stage most particulate solids are removed by sedimentation, and anaerobic conditions promote the anaerobic digestion of organic matter [3,4]. The second treatment stage is the UAF, a tank filled with a porous volcanic rock called tezontle. In this anaerobic unit, wastewater flows from the bottom part of the chamber to the outlet at the top. The porous rocks act as attachment media for the bacteria involved in organic matter degradation. Finally, the wastewater flows to the HFCW (Figure 1), which is a shallow pond (with a hydraulic depth of 60 cm) that is also filled with tezontle which acts as a support medium for the root development of the ornamental plant Agapanthus africanus planted at a density of three plants/m 2 . The mean hydraulic retention times of the ST, the UAF, and the HFCW were reported as 2.45, 6.4, and 11.75 days, respectively [3,4]. The WWTP at the site has been described in detail by de Anda et al. [3] and mathematically modeled by Fernandez del Castillo et al. [4]. The configuration of the treatment system is shown in Figure 2.

Water Quality Analyses
Four sampling points were established for monitoring water quality (SP1: at the sump pump; SP2 at the UAF inlet; SP3 at the UAF outlet; SP4 at the HFCW outlet), indicated by the blue boxes in Figure 2. Water samples were taken fortnightly for three months (January to March of 2020) to determine biological oxygen demand (BOD 5 ), chemical oxygen demand (COD), total Kjeldahl nitrogen (TKN), ammonia nitrogen (NH 3 -N), organic nitrogen (ON) and total suspended solids (TSS). These determinations were made by the Analytical and Metrological Services Unit (USAM) of the Centro de Investigación y Asistencia Tecnológica del Estado de Jalisco (CIATEJ), following the methods published by the Federation, W. E., and the American Public Health Association [31]. Samples were delivered to the laboratory less than one hour after they were taken. Nitrates (NO 3 − ) and nitrites (NO 2 − ) were measured by spectrophotometry using multiparametric kits TNT 835 and TNT 839, respectively (DR 5000, HACH, Loveland, CO, USA). System performance was evaluated by calculating the reduction in mass of each pollutant in each treatment stage and in the overall system. Temperature (Temp), pH, electrical conductivity (EC), and dissolved oxygen (DO) were measured at each sampling point using a multi-parameter probe (HI 9828, Hanna) to analyze the influence of physicochemical parameters on microbial communities.

Dna Extraction and High-Throughput Sequencing
Nine sampling points were established for DNA extraction and high-throughput sequencing ( Figure 2) through the entire treatment system. One sampling point was located in the ST, two sampling points were located in the UAF, and six were located in the HFCW. In the ST and the UAF (the anaerobic stages), samples included water mixed with sediments, while samples taken from the HFCW contained water combined with plant roots, substrate, water, and sediments.
Two sampling campaigns were conducted during 2020, the first in January and the second in March, months which are both within the dry season in Jalisco, Mexico. At each sampling point, and for each sampling campaign, three grab samples, each with three biological replicates, were taken, generating a total of 54 grab samples. Thus, three grab samples (and replicates) were taken from the second chamber of the ST each month (Figure 2; indicated as ST1 inside a red square). In the case of the UAF (Figure 2a), three grab samples (and replicates) were taken at the surface (UAF1), and three grab samples (and replicates) were taken at the bottom (UAF2) each month (adding up to a total of 12 samples for the UAF) ( Figure 2; indicated as UAF1 and UAF2 respectively, inside red squares). No composite samples were prepared for the UAF or the ST.
Regarding the HFCW (Figure 2a), samples were taken at a depth of 10 cm (Constructed Wetland Surface-CWS) and at a depth of 50 cm (CWB: Constructed Wetland Bottom) to assess variation associated with depth. Longitudinal distribution was also considered, as samples were taken at the inlet, in the middle, and at the outlet (at the two different depths previously described). The cross-sectional distribution of the microbial community was assessed by taking three grab samples at three points distributed cross-sectionally for each of the inlet, the middle, and the outlet (all at a depth of 10 cm) to prepare composite samples CWS1A, CWS1B, and CWS1C ( Figure 2b). Likewise, three grab samples were taken at three points distributed cross-sectionally for each of the inlet, the middle, and the outlet (all at 50 cm depth) to prepare composite samples (CWB1A, CWB2B, and CWB3C) as shown in Figure 2b. The three grab samples taken at CWS1A were used to prepare a (cross-sectional) composite space sample and the same procedure was followed to prepare composite samples for CWS2B, CWS3C, CWB1A, CWB2B, and CWB3C. The resulting composite samples correspond to each depth (10 and 50 cm) and each longitudinal point (inlet, middle, and output). For the HFCW, a total of six composite space samples were prepared (as 18 grab samples were collected) each month (a total of 36 samples).
Samples were stored at 4 • C while they were transported to the laboratory for processing. Once in the laboratory, samples were centrifugated at 800× g for five minutes to form a sediment pellet and 500 mg of the pellets were placed on a matrix for DNA extraction. Following the manufacturer's specifications, the DNA extraction procedure was performed using the FastDNA Spin Kit for Soil (MP Biomedicals, Solon, OH, USA). DNA extractions were stored at −80 • C until further analysis. High-throughput sequencing of the V3-V4 region of the 16S rRNA gene was performed by Novogene Corporation Inc. (Chaoyang District, Beijing, China) using an Illumina NovaSeq 6000 PE250 (paired-end to generate 250bp paired-end raw reads), obtaining 100k raw reads per sample employing the 341F (CCTAYGGGRBGCASCAG) and 806R (GGACTACNNGGGTATCTAAT) primers.

Bioinformatic Analyses
Bioinformatic analyses were performed using QIIME 2.0 (Quantitative Insights Into Microbial Ecology) software [32] following a standard bioinformatic pipeline. Demultiplexed sequences were denoised using DADA2 (p-trunc-len-f 0, p-trunc-len-r 0, p-trim-left-f 0, p-trim-left-r 0). After denoising, two characteristics tables (FeatureData[Sequence] and FeatureData[Taxonomy]) were constructed using 99% similarity with Greengenes 13_8 99% OTU full-length sequences [33,34]. Afterward, the classifier was trained using the length of the primers by the Naive Bayes classifier method. Finally, the taxonomic classification and the output of DADA2 (denoised sequences) with the trained classifier were aligned with classify-sklearn and the result was visualized as a taxa-barplot [35]. The sequenc-ing run has been uploaded to the Sequence Read Archive of the NCBI with accession number PRJNA700667.

Statistical Analysis
The sequencing depth of the 16S rRNA gene was presented using a rarefaction curve performed in R using the rarefy function, which is based on Hurlbert's formulation [36] and the standard errors proposed by Heck [37]. Barplots of relative read abundance were used to understand the composition of the bacterial community for the phylum and selected microbial families related to nitrogen and organic matter removal pathways. The DESeq2 R package was used to normalize read numbers [38]. Microbial groups with relative abundance <1% were grouped as "others" and unclassified families were denoted by "Un".

Principal Component Analysis
Principal component analysis (PCA) was performed for the bacterial families degrading nitrogen and organic matter in each unit (ST, UAF, and HFCW) to evaluate differences between treatment units as well as vertical variation within the UAF and both vertical and longitudinal variation within the HFCW. The aim of this analysis was to transform the original variables (bacterial families) to a new set of variables, the principal components, that are linear combinations of the original variables, which are uncorrelated and are ordered so that the first few of them account for most of the variation within the bacterial families [39]. Correlation biplots were used to further interpret bacterial community variation using the first two principal components in each case (variations between and within treatment units). In a biplot: (i) a vector represents a variable (bacterial family) and its length is proportional to the variance of the corresponding variable, (ii) angles between vectors reflect the correlation between the corresponding variables; (iii) points (observations) can be projected perpendicularly onto vectors, and the projection is indicative of the abundances of the families represented in the corresponding observations. The origin represents the average value. Projections in the same direction of the vector indicate values above average while projections in the opposite direction represent values below average [40].

Redundancy Analysis
The distribution of nitrogen and organic matter degrading bacterial communities between and within treatment units, and the effect of physicochemical parameters on bacterial communities, were analyzed using redundancy analysis (RDA) and the same analysis was performed on the variation in their composition. RDA is an extension of PCA that explicitly models response variables (bacterial families in this case) as a function of explanatory variables (physicochemical parameters in this case) [40]. The components in RDA are not only a linear combination of the response variables, but also of the explanatory variables. Correlation triplots were used to extract further information from the RDA results. A correlation triplot consists of two superimposed biplots that include quantitative explanatory and response variables (represented by vectors) and observations (represented by points) [40]. The interpretation of correlation triplots is as follows [41]: (i) The angles between two response variable vectors, or between two explanatory variable vectors, or between a response variable vectors and an explanatory variable, reflect their correlations, (ii) Points (representing observations) can be projected perpendicularly onto the response and explanatory variable vectors and indicate their values in the corresponding samples (observations).

Principal Coordinates Analysis
Dissimilarities in the composition of microbial communities between treatment units (ST, UAF, and HFCW) and dissimilarities between these communities at different (longitudinal and vertical) regions within the treatment units (the UAF and HFCW units) were analyzed using principal coordinates analysis (PCoA) based on Bray-Curtis distances. Permutational multivariate analysis of variance (PERMANOVA) (P < 0.05) and analysis of similarity (ANOSIM), using the Bray-Curtis dissimilarity index, were applied to evaluate statistically significant differences in the composition of bacterial communities between and within units (including vertical and longitudinal variation within units) [42]. The ANOSIM statistic R (which compares the mean of ranked dissimilarities between groups to the mean of ranked dissimilarities within groups) lies in the interval [−1,1]. A value of 0 indicates a completely random grouping, while positive R values suggest dissimilarity between groups and values below 0 suggest that dissimilarities are greater within groups than between groups [43]. All statistical analyses (PCA, RDA, PCoA, PERMANOVA and ANOSIM) were performed using the R software version 4.0.2, and the scales [44] vegan [45] packages. Graphics were made using ggplot2 package [44].

System Performance
The reduction in mass of each pollutant and the content remaining in the system effluent (BOD 5 , COD, TKN, NH 3 -N, ON, and TSS) are depicted in Figure 3. The system was found to be highly efficient in organic matter removal as it displayed significant mass reduction for both BOD 5 and COD, with average overall values of 548 ± 117 mg/L and 847 ± 181 mg/L, respectively ( Figure 3; Table S1), corresponding to average overall reductions of 90% of BOD 5 and 91% of COD. Higher degradation of organic matter occurred in the ST and the UAF, representing a combined reduction of 73% of BOD 5 and 72% of COD (mass reduction: 202 ± 60 mg/L of BOD 5 and 331 ± 93 mg/L of COD in ST; 243 ± 105 mg/L of BOD 5 and 345 ± 137 mg/L of COD in UAF). These results indicate that efficient anaerobic degradation occurs in the anaerobic stages (ST and UAF), in agreement with the DO and temperature levels measured in these stages (Table 1), which are optimal for anaerobic reactors [46]. Degradation of organic matter within the HFCW was also significant, with an average mass reduction of 171 ± 37 mg/L for COD, corresponding to a reduction of 18%. The DO concentrations measured in the HFCW were low but were higher than those observed for the anaerobic stages (Table 1). Table 1. Average values of physicochemical parameters measured.

OS (n = 24)
SP1 (n = 6) SP2 (n = 6) SP3 (n = 6) SP4 (n = 6) TSS removal was high for the system overall, reaching 339 ± 73 mg/L, which represents a reduction of 97% of the inlet load. However, the ST was found to make the highest contribution in comparison to the other treatment stages, with an average mass reduction of 258 ± 65 mg/L (74%), which proved that the ST unit performs its function efficiently, preventing the accumulation of solid particles in the succeeding units (UAF and HFCW) and preventing clogging (obstruction). Regarding nitrogen removal, the overall mass reduction was 173 ± 20 and 67 ± 8 mg/L for TKN and NH 3 -N, respectively. These values correspond to overall reductions (the contributions of the three stages combined) of 51% of TKN and 45% of NH 3 from inlet loads. This reduction is significant considering the initial concentration of TKN (337 ± 12 mg/L), as it is comparable to industrial loads; since the wastewater comes from an R&D industrial and biotechnology center, the effluent may have a higher concentration of nitrogen [47]. Similar results were reported by de Anda et al. [3] for this experimental system. As shown in Table S1, the highest reduction in nitrogen mass occurred in the HFCW, with values of 95 ± 26 mg/L (28%), 46 ± 13 mg/L (31%), and 48 ± 13 mg/L (26%) for TKN, NH 3 -N, and ON, respectively. Although the concentration of nitrates and nitrites entering the system was low (2.88 ± 1.69 mg/L for NO 3 − and 0.26 ± 0.06 for NO 2 ), the removal efficiency was high, reaching 93% and 96% for NO 2 − and NO 3 − , respectively. The performance of this experimental system has previously been characterized by Fernández del Castillo [4].

Diversity and Composition of the Bacterial Communities
High-throughput sequencing was performed on 52 samples (out of a total of 54 samples), from which 9,512,339 raw reads were obtained from the hypervariable region V3-V4 of the 16S rRNA gene, with a mean length of 420 bp. Sequencing depth was represented by a rarefaction curve of the 16S rRNA gene ( Figure S1). Two samples corresponding to the ST (one for each month) were not processed further because the DNA did not meet quality requirements. A total of 5,799,224 (61%) (from the 52 samples) were classified as bacterial employing the Greengenes database, and 3,713,115 (39%) reads were described as unclassified.
Proteobacteria was the most abundant phylum in all the system stages and gradually increased in abundance from one treatment stage to the next (ST < UAF < HFCW), with abundances of 34% in the ST, 38% in the UAF, and 60% in the HFCW (Figure 4a; Table S2). The members of this phylum are involved in a variety of metabolic pathways related to the carbon and nitrogen cycles [48][49][50], which explains their dominance in the whole WWTP. Additionally, Proteobacteria are known to be enhanced by the presence of macrophytes in CWS [51], which explains why their highest relative abundance was in the HFCW and which may also be associated with the high nitrogen mass reduction accomplished in this unit ( Figure 3). In addition to Proteobacteria, Firmicutes (26%), Bacteroidetes (14%), and Caldiserica (6%) were the most abundant phyla found in the ST. Firmicutes (21%), Caldiserica (14%), and Bacteroidetes (9%) were the phyla with the highest relative abundance in the UAF (besides Proteobacteria) while Bacteroidetes (9%), Actinobacteria (8%), and Chloroflexi (5%) presented the most significant abundances in the HFCW. As suggested by these results, the physicochemical conditions of each treatment stage affected its microbial composition. ST and UAF displayed a similar composition, as both are anaerobic units, while the bacterial composition within the HFCW, which contains macrophytes that provide exudates and dissolved organic matter, and transfer oxygen [52], presented more differences. The higher abundance of Firmicutes found in the anaerobic units in this study, compared to the HFCW, can be explained, as this phylum is known to degrade complex organic molecules [26,28]. The ST is the first unit in the system and is thus expected to receive a higher load of complex organic molecules, which diminished in subsequent stages and, accordingly, the relative abundance of Firmicutes decreases after each stage. Consequently, this phylum may be closely associated with the high rates of BOD 5 and COD mass removal reported for the first two units (the ST and UAF; Figure 3).
Besides Firmicutes and Proteobacteria, Acidobacteria, Actinobacteria, Bacteroidetes, and Chloroflexi phyla were present in all treatment units. These bacterial phyla have previously been reported in anaerobic bioreactors [24]. They are also known to conduct hydrolysis and acidification reactions in HFCW treating effluents with high organic loads, thereby playing an essential role in decomposing organic matter and nitrogen removal [53]. The phyla Gemmatimonadetes and Nitrospirae were only found in the HFCW, where higher levels of nitrogen and phosphate removal occurs. In consequence, previous studies have linked Gemmatimonadetes to nitrogen [54] and phosphate removal [50]. In this study, the phylum Gemmatimonadetes was only present at a depth of 10 cm in the HFCW and was completely absent in the first units (ST and UAF). Cheng et al. [26] reported that the Gemmatimonadetes prefer drier soils and disappear when the water content increases in a vertical subsurface flow CW, which explains why this phylum is present only in the upper layers of the HFCW of our experimental system, where the water content is lower because this is a CW with subsurface flow.
Only slight variations in abundance, in terms of phyla, were found in bacterial communities at different depths in the UAF (Figure 4b; Table S3). Regarding the effect of depth on the microbial communities found in the HFCW (Figure 4d; Table S4), the phyla Cyanobacteria, Gemmatimonadetes, and Nitrospirae were present only at 10 cm, which could indicate a dependency on higher oxygen concentrations [55,56]. However, it has been reported that the phylum Nitrospirae includes both aerobic and anaerobic microorganisms related to nitrification (nitrite oxidizers) [19] and methanogenesis [56], indicating that members of this phylum may be found in the upper and lower layers of diverse CW units. However, in the HFCW of our experimental system, bacteria of the phylum Nitrospirae were mainly present at a depth of 10 cm, possibly due to a limited accumulation of nitrite in the lower layers, since nitrite is known to be an unstable intermediate [57]. Significant longitudinal variation in the abundance of phyla was also found in the HFCW, where a decrease in Bacteroidetes and an increase in Actinobacteria were observed from inlet to outlet (Figure 4c; Table S5). The phylum Bacteroidetes is known to degrade complex organic compounds and includes members that can perform denitrification [58]. It is to be expected that the relative abundance of this phylum would decrease as substrate availability decreases from the inlet to the outlet of the HFCW. Furthermore, Gemmatimonadetes were present in the middle and the outlet sections but not near the inlet, which may be attributed to the low concentrations of DO found at the inlet. Some members of the Gemmatimonadetes have been regarded as aerobic heterotrophs capable of assimilating sugars [55,59].

Spatial Variations of Nitrogen and Organic Matter Degrading Families
To evaluate the spatial distribution of bacterial communities within the treatment systems, 18 nitrogen degrading families and 29 families that degrade organic matter were selected ( Figures 5-8). PCA biplots were developed to analyze the spatial variation in their relative abundance (Figures 6 and 8, only selected bacterial families vectors were plotted for illustration purposes) and a PCoA biplot ( Figures S2 and S3) was developed to demonstrate dissimilarities between the treatment units and spatial variation within the system stages. ANOSIM was performed to determine if such dissimilarities were significant (Table S14). As also found at the phylum level, the bacterial communities of the anaerobic stages (ST and UAF) were similar ( Figures S2a and S3a). The low levels of oxygen present in the ST and the UAF are necessary for anaerobic degradation and denitrification. In contrast, differences in the composition of the microbial communities between the HFCW and the anaerobic stages (ST and UAF) may be caused by the development of aerobic microenvironments and symbiotic interactions between bacteria and plants that occur in the HFCW [60].

Nitrogen Degrading Bacterial Families
Rhodocyclaceae was the most abundant family in all treatment stages (Figures 5a and 6a), suggesting that the oxygen concentrations found in all units (ST, UAF, and HFCW; Table 1;  Table S6) were sufficiently low for their growth. This family has been reported to participate in denitrification, organic matter degradation, and biofilm formation in high strength wastewater treatment systems under anaerobic conditions [7,61,62]. Comamonadaceae, Propionibacteriaceae, and Pseudomonadaceae were present in greater abundances in the ST in comparison to the UAF and the HFCW (Figures 5a and 6a).
Previous studies have suggested the participation of these families in denitrification is affected by higher concentrations of DO [63,64]. This condition of higher levels of DO can partially be found in the HFCW, since plant roots can infiltrate oxygen to this part of the system [65]. The output of the ST contained the lowest concentration of DO of all the sampling points (0.26 mg/L), which may explain the higher abundance of these families in this treatment stage. The primary function of ST is to separate the sludge, the effluent, and the scum layer from the wastewater, it removes suspended solids by retention and organic matter by anaerobic digestion [4,66], and reduces pathogen concentrations [3].
The Chromatiaceae were found to be present in the UAF at both depths (Figures 5b and 6b). The Chromatiaceae are a group of purple sulfur bacteria, commonly found in WWTPs, that utilize sulfide as an electron donor [67] and oxidize it to sulfate under anoxic conditions [68]. Therefore, the presence of this family suggests that denitrification and sulfate removal occur at the bottom of the UAF, and nitrification may occur in the UAF unit's upper layers. The abundance of the family Bacillaceae was higher in the ST and the UAF than in HFCW (Table S7). This family is known to be involved in nitrification under aerobic conditions but under anaerobic conditions members are able to perform denitrification [26,69], suggesting that denitrification occurs in the ST, based on the low levels found in this unit.
In this study the families Hyphomicrobiaceae, Nitrospiraceae, Rhodospirillaceae, Bradyrhizobiaceae, and Xanthomonadaceae were present at higher abundance at all sampling points within the HFCW unit compared with those in the anaerobic reactors (ST and UAF; Figures 5a and 6a). Aerobic microenvironments may influence this within the HFCW, these being attributed to activity in the plant rhizosphere. Figure 6d shows the vertical variation in the composition of bacterial communities (families) when comparing depths ( Figure S2d, ANOSIM R = 0.4933, Significance = 0.001) within the HFCW. The presence of plant roots and the permeation of atmospheric oxygen to a depth of 10 cm enable aerobic bacteria (such as nitrifying bacteria) to grow [70]. In contrast, at a depth of 50 cm anaerobic conditions allow other families to grow, such as denitrifying bacteria [4]. Within the HFCW, the Nitrospiraceae, Rhizobiaceae, and Xanthomonadaceae were the most abundant families at a depth of 10 cm (Figures 5d and 6d). These microorganisms have been regarded as nitrifiers, suggesting that nitrification occurs in the upper layers and near the plant roots, supporting microbial attachment, oxygen transfer to the rhizosphere, and root exudates [21]. At a depth of 50 cm, the Rhodocyclaceae, Chromatiaceae, and Bacillaceae occurred at a much higher abundance (in comparison to the upper layer, Figures 5d and 6d). Similarly, the family Rhodocyclaceae was more abundant at the outlet of the HFCW (Figures 5c and 6c). The family Rhodospirillaceae displayed an increasing trend from the inlet to the outlet of the HFCW (Table S9) and was more abundant at a depth of 10 cm than at 50 cm (Figures 5c and 6c). This sulfur oxidizing bacteria is capable of degrading organic compounds under anaerobic conditions and can act as a chemotroph under aerobic conditions [50]. However, Meyer et al. [68] reported a higher abundance of the family Rhodospirillaceae (16%) in environments with higher oxygen concentrations during wastewater treatment. These families may be affected when oxygen availability is low, which explains the increase in their abundance in the HFCW, where plants are known to create aerobic microenvironments [71]. Even if the Rhodospirillaceae includes anaerobes, our results suggest that aerobic conditions, partially found in the HFCW, are preferred by most family members. Longitudinal spatial variation within the HFCW were assessed by comparing samples from the inlet, the middle, and the outlet regions of the treatment unit ( Figure S2c; ANOSIM R = 0.2, Significance = 0.001). Variation was found in microbial family abundances from inlet to output, indicating that different metabolic pathways occur in each region of the HFCW. Families Xanthomonadaceae and Nitrospiraceae exhibited a decreasing trend from the inlet to the outlet (Figures 5c and 6c). These families are known to be aerobic bacteria [72] involved in nitrification [73] in significant symbiotic interactions with plants [74]. This decreasing trend could be related to the ammonia concentration which also decreased along the HFCW. The families Bradyrhizobiaceae, Hyphomicrobiaceae, Pseudomonadaceae, Rhizobiaceae, and Rhodospirillaceae presented higher abundances in the middle section of the HFCW than in the inlet and outlet (Figures 5c and 6c). The abundance of these bacterial groups may be related to substrate availability as they require carbon sources for growth and their abundance falls as carbon sources are consumed [68,75]. However, a preference for specific carbon sources could be suggested as the highest abundance of these families were observed in the middle part of the HFCW, where the decomposition of complex molecules may generate a significant amount of product, which can in turn be processed by these specific bacterial groups [76,77]. Characterization of the carbon sources along the HFCW could be useful in coming to understand the relationships between the abundance of bacterial families and the availability of specific carbon products and should be considered in future studies.

Families Degrading Organic Matter
The spatial variation of families that degrade organic matter within and between the treatment stages was analyzed (Figures 7 and 8; Tables S10-S13). The Ruminococcaceae was the most abundant family in the ST, followed by the Bacteroidaceae, Lachnospiraceae, and Porphyromonadaceae (Figures 7a and 8a). These families were also more abundant in the ST than in the UAF and HFCW. The most abundant family in the UAF was the Syntrophaceae, followed by the Clostridiaceae and Geobacteraceae. These families exhibited higher abundances in the UAF and were significantly less abundant in the HFCW. The families Ruminococcaceae and Lachnospiraceae have been related to the hydrolyzation of a variety of polysaccharides [78], acetogenesis [79], and fermentation [55], which are the first phases of degradation of the organic matter [80] present in raw wastewater. Thus, the higher abundance of these families found in the ST indicates that they degrade the more complex organic compounds present in the raw sewage into small molecules, later available for other microbial populations in the treatment units that follow. The Porphyromonadaceae has also been reported to be involved in the degradation of organic matter in anaerobic reactors [81,82]. The Clostridiaceae are commonly reported to be highly abundant in the effluent of anaerobic reactors such as the up flow anaerobic sludge blanket reactor [83,84].
Desulfovibrionaceae, Enterobacteriaceae, Mogibacteriaceae, and Moraxellaceae were present in both the ST and the UAF (the anaerobic stages). In contrast, Acetobacteraceae, Methylococcaceae, and Syntrophobacteraceae were only found in the UAF and the HFCW, and have been referred to as anaerobic fermentative bacteria [85] and have previously been found in other HFCW [86]. Similarly, the families Chitinophagaceae, Cytophagaceae, Methylophilaceae, Saprospiraceae, and Sphingomonadaceae were found exclusively in the HFCW (Figures 7a and 8a).
Van Lier et al. [87] reported that hydrolysis is usually the first and limiting step in the removal of organic matter, as it converts complex substrates into monomeric and dimeric compounds that form the substrates fed to the reactors that follow, and thus hydrolysis determines the overall removal of organic matter in the WWTP [87]. Similarly, Rajagopal et al. [88] reported that a pre-treatment process, such as the ST, is necessary to carry out hydrolysis since the degradation of particulate organic matter is slow and affects the performance of the following processes. The results presented here suggest that the hydrolysis step is performed within the ST, producing substrates that are available for the treatment units that follow, where fermentative and methanogenic bacteria were found.
The presence of Betaproteobacteria class members, such as Desulfobacteraceae and Desulfobulbacteraceae, was observed in all three units, with similar abundance (Figures 7a and 8a). Members of the phylum Proteobacteria, such as Betaproteobacteria and Deltaproteobacteria, have been reported to participate in the sulfur cycle in CW [22] and anaerobic digestion processes [68,89]. In addition, Meyer et al. [68] found that the presence of sulfur, potassium, manganese, total nitrogen, and total phosphorus in the wastewater contributed to variation in the structure and composition of families related to the sulfur cycle. The presence of these families in our experimental system indicates the occurrence of the sulfur cycle and the reduction of inorganic sulfur compounds [86].
The microbial composition differed between the upper (output) and the lower (input) layers of the UAF ( Figure S3b; ANOSIM R = 0.3056, Significance = 0.003). Families of obligate anaerobes were found in the lower layers of the UAF and facultative anaerobic bacteria were found in the upper layers. The abundance of families affected by oxygen, such as Clostridiaceae, Desulfovibrionaceae, Mogibacteriaceae, and Sinobacteraceae, decreased from the bottom to the surface (input to output) (Figures 7b and 8b). Other families, like the Bacteroidaceae, Desulfobacteraceae, Geobacteraceae, and Methylocystaceae, were more abundant at the surface (Figures 7b and 8b), indicating that they may include some facultative anaerobic species, the abundance of this families can be enhanced by a slight increase in oxygen levels. In addition to oxygen availability, it has been suggested that substrate availability and competitive interactions between microbial populations shape the distribution of the bacterial communities from inlet to outlet in anaerobic treatment units such as up flow anaerobic sludge blanket reactors [90,91].
Regarding vertical variation within the HFCW (Figures 7d and 8d), the families Acetobacteraceae, Geobacteraceae, Hydrogenophilaceae, and Syntrophobacteraceae presented a higher abundance at a depth of 10 cm than at a depth of 50 cm ( Figure S3d, ANOSIM R = 0.4575, Significance = 0.001). In contrast, the families with a greater abundance at a depth of 50 cm were the Clostridiaceae, Cytophagaceae, and Sphingomonadaceae ( Figure 8d). As previously mentioned, these changes may be related to the presence of root exudates and higher oxygen availability, which are more available at a depth of 10 cm. In the same vein, Krasnits et al. [92] reported that methanogenic bacteria and archaea were more strongly influenced by depth than by distance from the inlet in a study of the distribution of microbial communities in an HFCW. These authors attribute this tendency to the oxidationreduction potential (ORP) within the HFCW. The ORP parameter was not considered in this study, but an examination of this parameter is suggested for further studies.
Liu et al. [93] reported an increasing trend for the family Sphingomonadaceae along the flow path of an HFCW and a decreasing trend along the flow path of a VFCW, attributing this behavior to the direction of water flow, and highlighted the influence that the type of CW exerts on the structure of the bacterial community. This fact is consistent with our results (Figures 7c and 8c) where Sphingomonadaceae increased in abundance from inlet to outlet along the HFCW. However, in this study, few longitudinal differences were found in the HFCW with respect to the families related to organic matter removal ( Figure S3c, ANOSIM R = 1996, Significance = 0.001).

Effect of Physicochemical Parameters on Bacterial Communities
An RDA was performed to analyze the correlation between microbial families and the physicochemical parameters measured in the treatment stages (BOD 5, COD, TKN, NH 3 -N, ON, TSS, NO 2 − , temperature, DO, pH, EC; NO 3 − was omitted to avoid redundance in the analysis as it presented collinearity with NO 2 − ). The two main redundancy components explained 61% of the total variability in the bacterial community of nitrogen degrading families (a set of 18 families), while the two main redundancy components explained 70% of the total variability in the case of the families degrading organic matter. Figure 9 shows the RDA correlation triplot that explains the correlation between the physicochemical parameters, and selected bacterial families represented by a number (Table 2).  Table 2 presents the bacterial families assigned to each number. As mentioned above, the angles between the vectors of the bacterial families represent correlations between the families and angles between the vectors of water quality parameters represent correlations between these parameters ( Figure 9). Likewise, angles between the vectors of bacterial families and the vectors of water quality parameters represent correlations between them. Additionally, points representing each observation (sequenced sample) can be projected perpendicularly on the vectors of families or the vectors of the physicochemical parameters and give an indication of their corresponding values in such observations. The origin represents the mean value, projections in the same direction as the vector indicate values above average, and projections in the opposite direction represent values below average.
For both nitrogen and organic matter degrading bacterial families, DO, EC, and temperature displayed a higher influence on the bacterial communities in the ST and the UAF than on the communities present in the HFCW (Figure 9a,b). Organic matter removal is performed mainly by ammonification and methanogenesis in the ST and the UAF. The RDA results suggested that temperature affected the nitrogen removal related families within the ST and the UAF. It has been reported that nitrification processes can occur within a wide temperature range (16.5 to 32.5 • C) with an optimal range of 20 to 25 • C [73]. As shown in Table 1, the temperatures recorded at the output of the ST and the UAF were 21.06 ± 1.97 • C and 20.91 ± 2.31 • C, respectively. Although the temperature of the anaerobic stages is within the optimal range for denitrification, slight temperature changes may cause variations in the composition of the communities related to nitrogen removal. Temperature is also essential for anaerobic digestion, especially for the hydrolysis of complex organic compounds, whose breakdown is highly sensitive to temperature [87]. Methanogenic activity is reduced 10-20 times at low temperature (<15 • C) in comparison to activity at 35 • C [94]. Advantageously, wastewater treatment systems in tropical or subtropical regions are less affected by this issue, since wastewater temperature remains stable at above 25 • C throughout the year [95].
EC is used to measure the number of ions relative to salinity, a characteristic of wastewater that significantly affects bacterial communities in treatment systems [96]. Figure 9b shows how EC affects the bacterial families related to organic matter degradation in specific samples of the UAF. The most abundant families in the UAF (of those related to organic matter degradation) were the Syntrophaceae, Clostridiaceae, Sinobacteraceae, Geobacteraceae, and Ruminococcaceae (Figure 7a), suggesting that these bacteria are more sensitive to varia-tion in the salinity of wastewater. All these families have been identified in high salinity environments. The families Syntrophaceae and Clostridiaceae have been found in an acclimated marine sediment-derived culture used for biomethane production under high salinity conditions [97]. The presence and overgrowth of Sinobacteraceae has been related to high salinity conditions in shrimp culture enclosure ecosystems [98]. Further, the member Glk. subterraneus of the family Geobacteraceae, is a halotolerant bacterium that has been associated with high current generation (>1 A m −2 ) in electroactive biofilms used for microbial fuel cell technology [99]. Finally, the relative abundance of Ruminococcaceae is known to increase in salt stress conditions in a UAF used to produce methane from molasses wastewater [100]. The results found here suggest that our experimental system is tolerant of high salinity conditions, increasing its viability and applicability for treating different wastewater types.
Buffer capacity has been reported in anaerobic digestion processes, which explains why pH does not significantly affect bacterial communities in the anaerobic phases [101]. Higher nitrification rates have been reported in CW with pH ranges of 7.0-7.5 [101,102], while denitrification processes are optimal at a pH of 7.5 [103]. The pH measured at the outflow of the HFCW was 7.4 ± 0.4, which is within the optimal range for nitrification and methanogenic bacteria. The results also suggest that pH has a higher influence on communities in the HFCW than on those in the anaerobic stages ( Figure 9; Table 2; Tables S15 and S16). Variations in the pH of wastewater can cause stress in bacterial communities, as the intracellular pH of most microorganisms is close to neutral [104]. All the pH measurements recorded in this study (Table 1) were close to 7. However, variations in pH may occur inside the HFCW in specific regions (further pH measurements across the HFCW are required to prove this statement). Interestingly, bacterial communities at the inlet of the HFCW were the most sensitive to pH variation. The most abundant nitrogen degrading families were the Rhodocyclaceae, Xanthomonadaceae, Chromatiaceae and Comamonadaceae and the most abundant organic matter degrading families were Syntrophaceae, Methylococcaceae, Bacteroidaceae, and Sinobacteraceae. The Rhodocyclaceae and Comamonadaceae are denitrifying bacteria that use nitrate or oxygen as electron acceptors and short-fatty acids as electron donors [105], and these families have been reported to be the most abundant in other wetland systems [64,106]. Additionally, the Xanthomonadaceae is involved in forming microbial biofilm and granules as they participate in the production of extracellular polymeric substances [107]. Consequently, as these families are affected by pH variation, it is plausible that this parameter is a determining factor for nitrogen removal and the formation and stability of the wetland microenvironment.
It is reported that the Syntrophaceae family in syntrophic partnership with methanogens (Methylococcaceae) can degrade organics to H 2 and CO 2 to methane in methanogenic environments [108,109]. Therefore, manipulation of these families may prove useful in reducing the amount of methane produced by the system. However, further investigation is needed to determine the specific species involved in methane oxidation and the specific conditions required to favor their prevalence.
It can be observed in Figure 9 that vectors representing the concentrations of pollutants (BOD 5 , COD, TKN, NH 3 -N, ON, TSS, NO 3 − , NO 2 − ) are grouped in one single quadrant (quadrant II for nitrogen degrading families and quadrant I for organic matter degrading families). The observations made in the ST can be projected perpendicularly onto the vectors of the physicochemical parameters and give an indication of their corresponding values in such observations. In this case, greater concentrations were always observed at the initial stage (ST) and lower concentrations are found the next two treatment stages (UAF and HFCW).
Temperature, closely linked to climatic variation, is located in the same quadrant as the vectors of the physicochemical parameters. Thus, the vectors of the bacterial families that point in the same direction or in the opposite direction to the vectors of the water quality parameters and temperature (forming small angles or close to plane angles) are highly influenced by the pollutant inlet concentrations and variations in climate. Conversely, the vectors of bacterial families that form perpendicular angles with the vectors of the water quality parameters indicate that these families are barely influenced by variations in the water quality parameters. The families related to nitrogen degradation that were barely influenced by the physicochemical parameters comprised the Pseudomonadaceae, Bacillaceae, Propionibacteriaceae, Xanthomonadaceae, Caulobacteraceae, and Chromatiaceae. Similarly, the Syntrophaceae, Geobacteraceae, Desulfobacteraceae, Acidaminobacteraceae, Syntrophobacteraceae, Hydrogenophilaceae, Methylococcaceae, Xanthobacteraceae and Desulfomicrobiaceae were the organic matter degrading families that were mostly unaffected by water quality. It is important to note this group of bacteria can be considered to provide robustness to the treatment system.
Finally, correlations between specific bacterial families can be analyzed using the data presented in Figure 9. For example, a correlation between the Sinobacteraceae and the Saprospiraceae can be observed in Figure 9b, as the angle between them is close to zero. Similarly, these families are negatively correlated to nitrogen concentration (TKN, ON, and NH 3 -N), suggesting the existence of a mutualistic interaction between them by which nitrogen degradation is enhanced, however, this type of interaction must be further investigated. Furthermore, the families Sinobacteraceae and Saprospiraceae have been reported to be consumers of organic matter.
In addition, it can be deduced from Figure 9a that the occurrence of the Hyphomicrobiaceae and the Bradyrhizobiaceae is positively correlated, and they are negatively correlated with COD and BOD 5 concentrations. These correlations prove the interconnection between biological pathways for pollutant degradation. For instance, nitrogen and organic matter degradation are intimately related in the denitrification pathway, through which an electron is transferred from carbonaceous compounds to gaseous nitrogen [110]. Additionally, proteins contained in organic matter represent a significant amount of organic nitrogen. Protein degradation causes the accumulation of different nitrogenous compounds in wastewater [111]. Future studies should focus on the complex interconnections of pollutant removal pathways and microbial interaction networks.

Bacterial Communities in Multi-Stage Wwtps Located in Subtropical Regions
To the best of our knowledge, very few studies have focused on analyzing the structure and diversity of microbial communities in multi-stage decentralized WWTPs treating domestic wastewater in tropical and subtropical regions. Bedoya et al. [17] reported the microbial communities within the biosolids of a centralized WWTP treating municipal wastewater in Colombia. This plant uses an activated sludge system to treat wastewater generated by approximately 500,000 people (influent~1.8 m 3 /s). The authors found that Proteobacteria (66%) was the major phylum, followed by Actinobacteria, Firmicutes and Bacteroidetes. Similar results were obtained in the present study, in which Proteobacteria were the most abundant in all three treatment units of the system (34% in the ST; 38% in the UAF; 60% in the HFCW). Bacteroidetes and Firmicutes were also some of the most abundant phyla in anaerobic reactors (ST and UAF); additionally, Bacteroidetes and Actinobacteria were highly abundant in the HFCW.
The study of Desta et al. [28] reported the microbial communities within a multi-stage system treating tannery wastewater located in Modjo, Ethiopia. The system was integrated with two anaerobic reactors, followed by one aerobic reactor and a vertical-flow constructed wetland (VFCW) planted with Phragmites australis as a final step. These authors reported a relative abundance of 53% for Firmicutes and 24% for Proteobacteria in the aerobic reactor, a relative abundance of 52% for Firmicutes and 14% for Proteobacteria in the anaerobic reactor, and 44% for both Firmicutes and Proteobacteria in the VFCW [28]. In contrast, in this study, Firmicutes presented a lower relative abundance in the ST (26%) and the UAF (21%), while Proteobacteria showed a higher relative abundance in the ST (34%) and the UAF (28%). Although the system reported by Desta et al. [28] and that described in the present study are both in a subtropical region, the differences in abundance of these important phyla could be attributed to the wastewater characteristics and differences in system configuration, since in Desta et al. [28] the system described included an aerobic reactor that functioned as a pretreatment for the VFCW.
Song et al. [29] studied six activated-sludge WWTPs located in different climatic regions (tropical, subtropical, and temperate) and reported that Proteobacteria, Bacteroidetes, Chloroflexi, Acidobacteria and Nitrospirae were the major phyla in the WWTPs analyzed. However, Actinobacteria, Bacteroidetes, and Chloroflexi were more abundant in subtropical and temperate WWTPs compared to those in tropical regions. In the present study, Proteobacteria, Firmicutes, Bacteroidetes and Caldiserica were the most abundant phyla in the ST and the UAF, while Proteobacteria, Bacteroidetes, Actinobacteria and Chloroflexi were the most abundant in the HFWC. Based on these results, the microbial composition of the HFCW studied here is the most similar to that of the subtropical WWTPs studied by Song et al. [29]. Activated sludge systems and HFCWs are open processes highly affected by climate variations, while ST and UAF are confined spaces where more stable conditions can be maintained even with climate variability. In general, closed systems can present better conditions for bacteria susceptible to climatic variability. In the same study [29], the Nitrospirae phylum presented a higher abundance in moderately high temperatures, reaching higher nitrogen removal rates in tropical regions than in subtropical and temperate regions. This fact can be considered when improving the nitrogen removal rates in the decentralized system studied here. For example, a greenhouse could be installed to house the HFCW to provide higher temperatures inside the system, which could also enhance plant growth and allow for the cultivation of a wider variety of plants.
According to the RDA (Figure 9), family members of Mycobacteriaceae, Microbacteriaceae, Moraxellaceae, and Porphyromonadaceae are positively influenced by temperature. Conversely, temperature showed no significant effect on the families Rhizobiaceae, Nitrospiraceae, Methylococcaceae and Desulfomicrobiaceae. Previous studies have reported that the members of the family Nitrospiraceae are involved in nitrification, a process that may be affected by temperature changes in WWTPs located in regions with temperate climates. However, in tropical and subtropical regions, as in the case of the present study, Nitrospiraceae are not affected significantly by temperature, as the weather is warmer and more stable throughout the year. Accordingly, open systems, such as HFCWs, located in tropical/subtropical regions are adequate in maintaining more stable conditions for microbial communities throughout the year to enhance the wastewater bioremediation processes.
In general, a dominance of Proteobacteria, Actinobacteria, Bacteroidetes, Chloroflexi and Firmicutes can be found in tropical and subtropical regions, based on the studies of WWTPs located in these regions. However, the dominance of a specific phylum in these systems cannot be generalized, as their abundances vary widely from one study to another. Further studies in subtropical regions are required to generalize the microbial communities that may be expected in a WWTP, but these must also consider other factors that affect bacterial distribution, such as the type of wastewater, process features and design, and seasonality, among others. Detailed studies can also be useful in finding bacterial species that are favored in these regions and in developing strategies to use them to enhance the performance of decentralized technologies.

Conclusions
Passive wastewater treatment technologies are especially appropriate for tropical and sub-tropical regions as the climatic conditions facilitate adequate and stable performance throughout the year. Furthermore, the cost and infrastructure required to implement these systems in developing countries are acceptable. This study's contribution lies in the detailed description of the presence and distribution of the bacterial communities throughout a complete wastewater treatment system composed of ST, UAF, and HFCW. The pollutant removal pathways implemented by specific bacterial families within the different treatment stages were elucidated. Additionally, we described possible microbial interactions that enhance the removal of specific pollutants, as well as the influence of physicochemical parameters on the composition of the bacterial communities.
The characterization of the spatial variation in the microbial communities in our experimental system provides an in-depth understanding of the complexity of bacterial communities in a subtropical WWTP that combines an ST, a UAF, and an HFCW, and lays the foundation for future studies where manipulation of these microbial communities can result in better WWTP performance. Future studies should also focus on the complex interconnections of pollutant removal pathways and microbial interaction networks.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/w13091157/s1, Table S1. Removal efficiencies of the overall system and for stage; Figure S1. Rarefaction curve of 16S rRNA gene samples; Table S2. Percentage of reads classified as bacterial phyla > 1% using Greengenes 13_8. Comparison between the three treatment stages (ST, UAF and HFCW) for nitrogen degrading families; Table S3. Percentage of reads classified as bacteria phyla > 1% using Greengenes 13_8. Vertical variations within the UAF; Table S4. Percentage of reads classified as bacteria phyla > 1% using Greengenes 13_8. Vertical variations within the HFCW; Table  S5. Percentage of reads classified as bacteria phyla > 1% using Greengenes 13_8. Longitudinal variations within the HFCW (Inlet, Middle and Outlet); Table S6. Percentage of reads classified as bacteria family > 1% using Greengenes 13_8. Comparison between the three treatment stages (ST, UAF and HFCW) regarding the nitrogen degrading families; Table S7. Percentage of reads classified as bacteria family > 1% using Greengenes 13_8. Vertical variations of nitrogen degrading families within the UAF; Table S8. Percentage of reads classified as bacteria family > 1% using Greengenes 13_8. Vertical variations of nitrogen degrading families within the HFCW; Table S9. Percentage of reads classified as bacteria family removal > 1% using Greengenes 13_8. Longitudinal variations of nitrogen degrading families within the HFCW (Inlet, Middle and Outlet); Table S10. Percentage of reads classified as bacteria family > 1% using Greengenes 13_8. Comparison between the three treatment stages (ST, UAF and HFCW) for organic matter degrading families; Table S11. Percentage of reads classified as bacteria family > 1% using Greengenes 13_8 involved in organic degradation at vertical variations within the UAF; Table S12. Percentage of reads classified as bacteria > 1% using Greengenes 13_8 involved in organic degradation of vertical variations within the HFCW depth; Table S13. Percentage of reads classified as bacteria family > 1% using Greengenes 13_8. Longitudinal variations of organic matter degrading families within the HFCW (Inlet, Middle and Outlet); Table S14. PerMANOVA and ANOSIM analyses for the vertical (UAF and HFCW) and longitudinal variations (HFCW) for nitrogen and organic matter degrading families; Figure S2 Table S15. RDA for nitrogen degrading families (first four significant redundancy components; Table S16. RDA for organic matter degrading families (first four significant redundancy components).

Data Availability Statement:
The sequencing run data presented in this study are deposited and public available at the Sequence Read Archive on the NCBI with the following link: https://www. ncbi.nlm.nih.gov/bioproject/PRJNA700667 (accessed on 12 April 2021).

Acknowledgments:
The authors appreciate the effort and excellent work carried out by CIATEJ's analytical services unit, who at all times supported the timely registration and processing of samples for analysis. Also, the authors are grateful to Tecnológico de Monterrey, Campus Guadalajara for allowing access to the laboratories and facilities of the Department of Bioengineering.