The Aojiang River (58 m³/s discharge) is located in Lianjiang County, Fujian, China. Tidal flats are in the south of the Aojiang Estuary (Figure 1A). The nearshore water is very turbid (yellowish), and the semi-diurnal tides show a range of 3.8 m (China Oceanic Information Network, http://www.nmdis.org.cn). During low tide, vast flats of up to a few kilometers long emerge from the water.

Sinonovacula constricta and R. philippinarum are the two species that are cultivated in the flats south of the Aojiang River mouth (Figures 1B–D). Sinonovacula constricta prefers finer sediments, digging 30–40 cm into the sediment, whereas R. philippinarum prefers coarser sediments, digging shallower (~10 cm). The S. constricta cultivation area is 705 hectares, with an annual production of 20,667 tons, while the R. philippinarum cultivation area is 525 hectares, with an annual production of 19,200 tons (Committee, O.O.L.L.C.C., 2017). Due to the bivalves vertical migration and harvesting, top 0−10 cm sediment were repeatedly under strong mixing or disturbance. Each year, larvae of both species are sowed in spring, and harvest occurs in late summer. For S. constricta, the cultivation area is manually maintained like a paddy field but without paddy plants. Sinonovacula constricta are cultivated in a field surrounded by a ridge (Figure 1C). The ridges are maintained by inserting Miscanthus plants into the surface sediment. Miscanthus are manually carried to the fields via tricycles and regularly re-supplied during low tide to maintain the ridge shape.

The bivalves feed themselves via burrows, which enhanced the vertical mixing of sedimentary materials. Due to the bivalves vertical migration and harvesting, top sediment (e.g., 0-10 cm) were repeatedly under strong mixing or disturbance. In the field, we saw “black soils” randomly distributed from very surface to the deep (~40cm deep), without any lamination structure. We repeatedly found this profile distribution feature in cultivating and non-cultivating seasons.

The external carbon added to the 0-10 cm top sediment is composed of sedimentation from overlaying water and from benthic algae production in the very surface of sediment. For the two months after harvest, these added carbon is estimated to account for a very trace proportion (2.7%) of the total carbon pool. Further given that the added carbon (and microbe) took place in the very surface of the sediment, the top 0-10 cm sediment after harvesting can be regarded as a close system. The detailed calculation process can be found in Supplemental Materials.

The S. constricta and the R. philippinarum are harvested in summer (late August early September) every year. Field sampling in this work was conducted in late October, 2020 (Figure 1B). For the S. constricta cultivation area (station YC), three subsamples were collected in the field where S. constricta was cultivated and another three subsamples were collected on the ridge (S. constricta-free). For the R. philippinarum area (station TT), two sites were chosen. One site was more seawards and the other more landwards. Again, three subsamples were collected for each site. Mixed surface sediment samples (0–10 cm deep) of the flats were collected using a clean spoon during low tide. Sediment samples were stored frozen until laboratory determination.

The DI was calculated by the following equation (Dauwe et al., 1999):

where var_i, AVG var_i, STD var_i, and fac.coef._i are the mol%, mean, standard deviation, and factor score coefficient of AA_i, respectively. Factor score coefficients were calculated using principal component analysis.

For the Shannon index, 2 was the log base, and for the Simpson index, the equation is 1 − d o m i n a n c e = 1 − ∑ p i 2 where p_i is the proportion of the community represented by OTU_i. For the biological and environmental data, a redundancy analysis (db-RDA) was conducted via PRIMER 7. In the PRIMER settings, the distance method is Bray-Curtis similarity, selection criterion is adjusted R² and selection procedure is step-wise.

The network analysis was conducted via molecular ecological network analysis pipeline (Deng et al., 2012). The network analysis followed the default settings and the RMT threshold was set as 0.95 (p = 0.001, r² = 0.5). The visualization of the network analysis results was further performed by Cytoscape (Shannon et al., 2003).

Mean grain size in S. constricta sites (field + ridge) and R. philippinarum sites (landwards + seawards) was 46.5 ± 15 μm and 146 ± 37 μm, respectively. Within the S. constricta sites, field sediment was finer (YC1, 35 ± 1 μm) compared to that from the ridge (YC2, 58 ± 13 μm) (Table 2). Within the R. philippinarum sites, larger grain size sediment was found in landward samples (TT2, 174 ± 34 μm) relative to seaward samples (TT1, 119 ± 8 μm) (Table 2). This was consistent with the grain size preference of S. constricta (silt) and R. philippinarum (sand), respectively. For sedimentary Cl^- and S O 4 2 − , a seaward increase was recorded, namely elevated values were found in samples from seawards sites [TT1, (Cl^-) = 208 ± 12 mmol/kg, ( S O 4 2 − ) = 11 ± 0.4 mmol/kg] compared to landwards sites (TT2, (Cl^-) = 160 ± 11 mmol/kg, [ S O 4 2 − ] = 9 ± 0.7 mmol/kg); however, there were no such significant differences found between the S. constricta field and ridge (Table 2). For sedimentary N O 3 − , the highest values were observed in the S. constricta field (YC1, 0.48 to 1.0 mmol/kg), with the rest of the samples (YC2, TT1, and TT2) ranging between 0.33 and 0.45 mmol/kg. Statistically, sedimentary Cl^-, S O 4 2 − , and N O 3 − showed no clear difference between S. constricta and R. philippinarum sites, but there was a significant difference in grain size.

Regarding bulk carbon, S. constricta sites showed a higher carbon content relative to R. philippinarum sites (i.e., 0.46% ± 0.02% vs. 0.16% ± 0.02%). Within each site, there was no clear difference in carbon content between the field and ridge, or between landwards and seawards sites (Table 2). The lowest δ¹³C (-27.22‰) was recorded in the ridges (Table 2).

Sedimentary total AAs ranged from 1.8 mmol/g OC to 5.6 mmol/g OC (OC: organic carbon). Samples from R. philippinarum landward sites (TT2) showed the lowest AAs (1.9 ± 0.2 mmol/g OC) and samples from S. constricta field ridges (YC2) showed the highest values (3.4 ± 2 mmol/g OC) (Table 2). Glycine, the non-chiral and most simple AA compound, was the most abundant (mean 1.2 mmol/g OC), followed by L-glutamic acid (0.19 mmol/g OC) and L-aspartic acid (mean 0.18 mmol/g OC). Regarding D-AAs, which are usually related to microorganism activities, D-arginine and D-aspartic acid were the two most abundant, with mean concentration as high as 22 μmol/g OC and 20 μmol/g OC, respectively. Correspondingly, the D/L ratio of arginine (D/L arg; mean 0.33) and D/L ratio of aspartic acid (D/L asp; mean 0.12) was the highest among chiral AAs. D/L ratios of key AAs in samples from S. constricta sites were higher (arginine, serine, and aspartic acid) or similar (threonine and alanine) when compared to samples from R. philippinarum sites. Within S. constricta sites, the field samples showed higher D/L ratios compared to the ridge samples, and within the R. philippinarum sites, a seawards decrease of D/L ratio was observed (namely D/L landwards sites > D/L seawards sites). The D/L ratio patterns were similar to those of the quantitative PCR results for bacteria and archaea. Namely, bacteria and archaea showed the highest abundance in S. constricta field samples (1.8 × 10⁹ copies/g wet weight [ww] and 1.0 × 10⁸ copies/g ww, respectively), followed by S. constricta ridge samples and R. philippinarum landwards site; the R. philippinarum seawards sites had the lowest bacteria (9.5×10⁸ copies/g ww) and archaea (4.6 × 10⁷ copies/g ww) abundances.

The mean DI of all 12 sediment samples was -0.91. The DIs of samples from R. philippinarum sites were commonly lower compared to those of S. constricta sites, and the lowest DI was found in R. philippinarum landwards sites (-1.06, TT2a). The highest DI (-0.47) was recorded on the ridge (YC2b), but overall DI within the S. constricta field and ridges was comparable (Table 2).

Statistically, carbon parameters (DI, δ¹³C) were significantly different between R. philippinarum sites and S. constricta sites, but AA composition did not significantly differ.

The microbial (bacteria and archaea) abundance obtained by quantitative PCR (copies/g ww) was significantly different between R. philippinarum sites and S. constricta sites (p = 0.002) (Figure S1). The bacterial and archaeal biodiversity quantified by 16s rRNA gene high-throughput sequencing also had significant difference between R. philippinarum sites and S. constricta sites, respectively (Figure S2).

Taxonomic assignment suggested that family Thermoanaerobaculaceae was the major bacterial group in S. constricta sites, followed by Desulfobulbaceae and Desulfobacteraceae; Similarly, the bacterial community of R. philippinarum sites, was also dominated by Thermoanaerobaculaceae, followed by NB1-j, Wosesiaceae and Microtrichaceae (Figure 2A). There was a clear difference between R. philippinarum sites and S. constricta sites as revealed by hierarchical clustering analysis (Figure 3A).

For archaeal fraction, family Nitrososphaeria was the most abundant archaeal group in R. philippinarum sites, followed by Bathyarchaeia and Woesearchaeia. The top two most abundant communities belonged to Bathyarchaeia and Nitrososphaeria, followed by Woesearchaeia, Methanosarcinaceae and Lokiarchaeia (Figure 3B). Bathyarchaeia were very abundant in samples from S. constricta sites, but showed a much lower abundance in R. philippinarum sites (Figure 2B). Similarly, a clear assemblage difference was observed between the R. philippinarum and S. constricta sites based on hierarchical clustering analysis (Figure 3B).

The db-RDA was applied to archaea and bacteria, and yielded a good r² (both r² > 0.8 for archaea and bacteria; Figure 4). For bacteria, DI and N O 3 − explained 38% and 15% of the variation, respectively (Figure 4A). For archaea, DI explained 42%, and N O 3 − 22% of the variation (Figure 4B). Both archaea and bacteria are primarily constrained by DI and N O 3 − . Instead, the carbon source as quantified via δ¹³C was not the key influencing factor (Figure 4).

At family level, the relation between OTUs and DI is shown in Figure 5. For bacteria, 22 out of 350 families showed a negative relationship with DI (i.e., the OTU was negatively related to DI), with a subtotal OTU reading of 270,636 in comparison to a total bacterial OTU reading of 541,048. That is, 22% of bacteria families occupied 50% of the total OTU readings (Figure 5A). The remaining 78% families contributed to the rest 50% of OTUs (Figure 5A), and their OTUs were either positively related to DI or showed no relationship with DI (Figure 5A). With respect to archaea, the unbalanced composition was even more pronounced (Figure 5B). Namely, only two out of 85 families showed a logical OTU relationship with DI, which accounted for 57% of OTUs (Figure 5B). The top two families of archaea were Nitrosopumilaceae (OTU = 329001) and an unidentified family of Desulfurococcales (OTU = 644).

Network analysis result is shown as Figure 6. Edges in the network (lines that connect nodes) indicate the statistical relationship between two given nodes. They can be negative (np; brown in color) or positive (pp; cyan in color) (Figure 6). A first look at the network is complicated (Figure 6A). Given that carbon degradation is the result of microorganism utilization, and by this utilization microorganisms grow and reproduce, an increase in microorganism biomass should be logically accompanied with advanced carbon degradation. Therefore, microorganism that showed negative relation with DI is termed as winner in the network analysis. On the opposite, microorganism that showed positive relation with DI is termed as loser, and microorganism that showed no relation with DI is termed as centrist. Accordingly, the nodes in the network analysis were further grouped into red (winner), yellow (loser), and blue (centrist) according to their relation with DI (Table 3, Figure 6B). The negative edges (brown lines) were the dominant connections for winner external relations in the network (Figure 6B). But the highest positive versus negative edge ratio was observed for loser external relations (513 pp vs. 356 np; pp-positive edge, np-negative edge; Table 3). That is, the ratio (513/356) is higher when compared to that of centrist and winner (457/425 and 247/717, respectively; Table 3).

AAs belong to the labile pool in bulk carbon, and their individual synthetic and metabolic pathways differ. Physiologically, some AAs (e.g., the simple-structured glycine) are at the base level (being the substrate) in many synthesis pathways while other AAs (e.g., the complex-structured phenylalanine) are the target product of several synthesis steps. In addition, individual AAs have varying associations with the cell wall or cytoplasm, or play different roles in embedding into structural matrices (Cowie and Hedges, 1996). As a result, during microorganism utilization and carbon degradation, certain AAs tend to accumulate (increase in mol%), whereas others tend to be lost (decrease in mol%) (Cowie and Hedges, 1992). Based on such varying accumulation/removal patterns, the AA-based molecular indicator, DI, is a statistical term that quantifies the carbon degradation degree based on individual AA mol% changes (Dauwe and Middelburg, 1998). It is much more specific to carbon degradation status relative to OC%, and has been widely applied in carbon degradation studies (Amon et al., 2001; Arnarson and Keil, 2001; Dittmar and Kattner, 2003; Amelung et al., 2008; Davis et al., 2009; Peter et al., 2012). AAs show a similar mol% distribution pattern among various organisms ranging from bacteria/fungi to plankton and vascular plant tissues (Cowie and Hedges, 1992). Fresh live matter, including bacteria, shows a DI value close to +1.5, while that of heavily degraded sedimentary carbon can be as low as -2.2 (Dauwe and Middelburg, 1998). Such feature means that DI is a good biomarker for quantifying carbon degradation status in marine organic biogeochemistry studies.

Bacteria biomass interference to bulk carbon is assessed via D-alanine concentration, assuming that all D-alanine is dominantly contributed by bacteria due to its key role in peptidoglycan construction. Using a conversion factor of 108 nmol/mg bacterial carbon (Kaiser and Benner, 2008), the fraction of bacterial carbon to bulk carbon in the tidal flats ranged from 7% to 15%, with a mean value of 10%. Therefore, bacterial carbon is one order of magnitude lower relative to bulk carbon. Given the additional D-AA sources in higher plants (Strauch et al., 2015), and algae (Yokoyama et al., 2003), the above estimated bacterial carbon fraction in bulk carbon would be even smaller. Furthermore, D-AAs have been detected in some archaea (Nagata, 1999), but their interference with bacteria carbon estimations, as well as the determination of archaea carbon contribution to bulk carbon, has not been assessed due to lack of a conversion factor.

The db-RDA analysis suggests that carbon degradation and related nitrogen cycling (e.g., ammonia-oxidation) are the dominant two processes that microorganisms were responding (Figure 4). It is not surprise to see the bacteria (Figure 2A) are related to carbon cycling as many of the bacteria are heterotrophic or related to carbon degradation. For example, thermoanaerobaculaceae prefer sugar, organic acids, and protein-like materials (Dedysh et al., 2021). Woeseiaceae are widely involved in the cycling of detrital proteins in marine benthic environments (Hoffmann et al., 2020). Sandaracinaceae may be related to utilization of ethanol and fatty acids (Probandt et al., 2017) and actinomarinales are usually characterized by their utilization of AAs as a carbon source (López-Pérez et al., 2020).

As nitrogen cycling process, the most abundant archaea, Nitrosopumilaceae, is related to nitrogen cycling (Figure 2B). Candidatus Nitrosopumilus were the most abundant archaea at the species level; these archaea are chemolithoautotrophic and grow by aerobic ammonia oxidation (Qin et al., 2016). The strong negative relation between Nitrosopumilaceae OTUs and DI (Table 3) indicates that the ammonia oxidizing archaea benefits from the carbon degradation. The source of ammonia in the sediment remains unknown in the present study, but the degradation of organic matter is a credible source, in which organic nitrogen is decomposed as the corresponding inorganic form ammonia.

The second most abundant archaea, Bathyarchaeia, are directly related to the carbon cycle. They are capable of utilizing a variety of organic carbon from labile to refractory sources (Zhou et al., 2018a), and some Bathyarchaeia are well-known for their growth being stimulated by lignin, a series of refractory and aromatic polymer compounds in higher plants, and such stimulation for growth is even higher compared to protein source stimulation (Yu et al., 2018). This is consistent with the extensive use of higher plant Miscanthus in S. constricta sites (Figure 1C).

Both uneven distribution between family numbers and OTUs (Figure 5) and contrasting external connections between winner, loser and centrist (Figure 6B) suggest severe competition among microbes during carbon utilization. On one hand, it is not surprising to see negative edges between winner and loser groups (Figure 6B) since these two groups had an opposite relationship with DI, but the edges between winner and centrist were also negative, with very few exceptions (Figure 6B). This means that the winner microorganisms that contributed to and reproduced under carbon degradation were overall suppressing both the losers (yellow nodes) and the centrists (blue nodes). On the other hand, the highest pp/np ratio for losers (Table 3) is achieved by establishing positive connections with centrists, which suggests the loser’s survival strategy: uniting centrists.

Regarding how microorganisms respond to changing organic matter, there are two hypotheses. The first one is the taxon-level accessibility model. This suggests that different microorganism populations specialize in the metabolism of certain carbon pools, and hence changes in carbon degradation status would be accompanied with observable changes in microorganism community structure (e.g., Teeling et al., 2012). The second is the cellular-level accessibility model. This suggests that cellular level processes are more important in dealing with variable carbon and hence microorganisms can process multiple carbon pools as needed, thereby observable microbial community changes can hardly be found with changing carbon degradation status on a short time scale (e.g., Mou et al., 2008). Although the microbial composition remained unchanged in the beginning of sedimentary carbon incubations, at a later stage (incubation 4–8 days), as the pool of available organic matter became more recalcitrant, an increasing abundance of bacterial taxa that were able to produce the enzymatic machinery for processing complex biopolymers was observed, supporting the taxon-level accessibility model (Mahmoudi et al., 2017). In sediment-derived dissolved carbon incubations, the microbial community in the early stages of incubation was influenced by relatively labile tannin- and protein-like compounds; while in later stages the community composition evolved to be most correlated with less labile lipid- and lignin-like compounds (Wu et al., 2018).

These previous findings (Mahmoudi et al., 2017; Wu et al., 2018) support the notion of competition being the feedback to carbon degradation status (Figures 5, 6). Two months after harvesting, the microorganism and labile carbon degradation coupling was expected to be closer to the steady state, rather than still at a unsteady or initial stage. In our work, the microorganisms showing a negative relationship with DI (red nodes; Figure 6) are the key implementors and hence beneficiaries of carbon degradation (winner). The winner gained the most energy and substrates from the degradation of carbon, which thereby placed them at an advantage in the competition of the entire microbial community. In contrast, the negative or poor performers in carbon degradation process (yellow and blue nodes) failed to gain maximum resources from carbon degradation and hence were at a competitive disadvantage in the community (Figure 6; Table 3). Within the winners (red nodes), it is expected that the most abundant archaea, Nitrosopumilaceae, occupy a unique niche, namely oxidization of carbon degradation products, which enables them to benefit from the carbon degradation process (Figure 2).

Most bacteria in the flats were also identified in previous incubation studies, so a comparison with previous work is at work. Variovorax, Methylotenera, Sphingomonas, and Rhodopseudomonas (specie level) were all identified in tidal flats (present study) and previous sediment-derived carbon incubations (Wu et al., 2018). However, these bacteria showed an increased abundance with advanced carbon degradation in previous incubations (Wu et al., 2018), whereas in tidal flats they either showed no relationship or a decrease abundance with advanced carbon degradation. In other words, these bacteria were winner in previous studies but became centrist or loser in current work. A similar discrepancy also occurred regarding Bacillales and Paenibacillaceae when comparing the findings of tidal flats results (present study) and previous sedimentary carbon degradation incubations (Mahmoudi et al., 2017). One possible explanation is that the chemical and biological habitats differed and hence so did the given bacteria performance. However, the varying performance of the same family or specie during carbon degradation process in different environments also implies that the microbial community structure under a certain carbon degradation situation is the specific outcome of that given competition, and there is no specific family or genus that can always beat others in the competition. If this assumption is correct, it is expected that the opportunity for a structure and network shift lies in changing external conditions (e.g., addition of new carbon sources), which creates another round of competition, while it is interesting to check the performance of microorganisms that occupy a unique niche in the competition, like Nitrosopumilaceae in this work.

Tidal flats of S. constricta sites stored carbon in a higher efficiency than R. philippinarum sites due to the higher carbon content (Table 2). However, from a carbon degradation perspective, S. constricta sites showed higher carbon degradation potential (higher DI values), whereas carbon in R. philippinarum sites were likely to lose less due to a more advanced degradation status (smaller DI values) (Table 2).

In microbial ecology, the bacteria competition strategy has been previously addressed following macroecology theory (Andrews et al., 1986; Klappenbach et al., 2000), namely that microorganisms can be regarded as r-strategists (copiotrophs) and K-strategists (oligotrophs). Oligotrophs have low growth rates and high resource use efficiency, and hence are characterized by an ability to grow under low substrate concentrations, whereas copiotrophs are the opposite, with higher growth rates and lower resource use efficiency. For young/disturbed ecosystems, the K/r ratio tends to be lower compared to mature/undisturbed ecosystems (Odum, 1969; Odum, 1985). It is indicated that Acidobacteria and Actinobacteria behave as K-strategists, whereas Proteobacteria and Bacteoidetes are r-strategists, based on their contrasting relationship between carbon mineralization rate and bacteria biomass (Fierer et al., 2007). Previous research has shown that the tracking of microorganisms at finer phylogenetic and taxonomic resolutions (e.g., family level or lower) is needed (Ho et al., 2017), but such classification remains lacking at present. Regarding phylum level, the K/r ratio for S. constricta sites and R. philippinarum sites was 0.27 ± 0.04 and 0.48 ± 0.03, respectively. Such a K/r ratio pattern indicates that microbial ecosystem of the S. constricta sites was younger and more disturbed, whereas that of R. philippinarum sites was more mature and less disturbed. In a forest ecosystem, the microbial system changes from a K-dominated to a r-dominated community stimulates carbon degradation (Zhou et al., 2018b), and if this is extrapolated to the tidal flat system, this would mean that S. constricta sites experience higher carbon degradation rates (and hence result in greater CO₂ release to the atmosphere) compared to R. philippinarum sites. This is in agreement with the carbon degradation potential derived from DI (Table 2). It should be noted that the overall net CO₂ release assessment in the respective aquacultural tidal flat sites requires additional direct CO₂ gas measurements and sediment-water interface dissolved carbon flux observations.