Linking Bacterial Communities to Optical-Derived Properties of Porewater DOM in Sediments in the Coastal East China Sea

Bacterial communities and porewater dissolved organic matter (DOM) pool are intimately interactive in sedimentary environments. Estuarine coastal regions are an interactive area between terrestrial and marine influences in terms of DOM origins and freshness. Yet, we know little about the relationships between the bacterial communities and DOM in those regions. In this study, porewater DOM samples were collected from 42 sites in the coastal East China Sea. The porewater DOM optical properties were determined by fluorescence and absorption spectra, while the corresponding bacterial community compositions of those sediments were examined by 16S rRNA gene high-throughput sequencing. The results showed that bacterial species richness was positively correlated with multiple terrestrial indicators based on the optical properties of DOM, which implied that heterogeneous DOM from terrestrial origins might harbor a wider spectrum of bacterial taxa in marine sediments. The analysis of the co-occurrence network of the bacterial communities showed that the edges and density for samples with low DOM freshness were 3.4 times and 3 times those for the samples with high DOM freshness, respectively. This suggested that the connection among the bacterial taxa under the lower DOM freshness condition were enhanced and that reduced freshness of DOM may encourage more complimentary utilization of resources. The findings provide a new insight into such interactive processes of heterogeneous organic matter utilization meditated by microorganisms in coastal sediments.


INTRODUCTION
Marine dissolved organic matter (DOM) represents one of the largest reservoirs of reduced carbon on the surface of the earth, matching the size of the atmospheric reservoirs (Hansell et al., 2009). In marine sediments, DOM undergoes a series of transformations due to complex interaction with other components in oceans and plays an important role in global carbon cycling (Chen and Hur, 2015). DOM in sediments is traditionally considered to be recalcitrant (Weston and Joye, 2005;Chipman et al., 2010;Robador et al., 2010). However, some recent studies revealed that some parts of DOM in marine sediments have the potential to be bio-reactive (Rossel et al., 2016;Cai et al., 2019;Tobias-Hünefeldt et al., 2021). Since prokaryotes (especially for bacteria) are the pivotal component responsible for transformation of DOM in marine sediments (Wang et al., 2021), it is important to investigate the interactions, between DOM and bacteria in marine sediments. For all heterotrophic and some autotrophic bacteria in sediments, stepwise processes during the decomposition of sedimentary DOM sever as their main carbon and electron donors source (Burdige and Komada, 2015;Mahmoudi et al., 2017), and porewater DOM was believed as the most active delegate in sedimentary DOM pools (Derrien et al., 2019). Two tangled groups of bacteria processes could be summarized from conceptual frameworks: consumption and production. On the one hand, bacteria could utilize DOM as both carbon sources and electron donors (Coates et al., 2002;Ijiri et al., 2012). DOM could be either remineralized or incorporated into bacterial biomass through bacterial metabolism (Azam et al., 1983;Ducklow et al., 1986;Goldman and Dennett, 2000). On the other hand, chemoautotrophs could transform inorganic matter into DOM directly, while heterotrophs also contributed to DOM in the form of extracellular secretion like hydrolytic enzymes and siderophores (Vraspir and Butler, 2009;Arnosti, 2011). In most cases, consumption and production were coupled. For instance, microbial carbon pump proposed the important role of microbes in transforming the labile DOM into refractory DOM (Jiao et al., 2010).
Many studies have found that specific DOM compounds or compound classes may have the effect to shape specific bacterial groups, while different bacteria have preferences for specific DOM portions (Cottrell and Kirchman, 2000;Dyhrman et al., 2009;Poretsky et al., 2010;Amaral et al., 2016). However, specific DOM substrates cannot represent the chemical complexity of natural DOM for lack of environmental relevance. Similarly, natural bacterial communities generally have high biodiversity and hence possess diverse ecosystem functioning (Smith et al., 2018). Therefore, it is necessary to give both DOM and bacteria more thorough descriptions in order to have a holistic and integrative understanding of their relationship. Although it is difficult to couple complexities of DOM and bacterial community considering chemodiversity and bacterial diversity, recent advances in DOM characterization and high-throughput sequencing make it possible to do comprehensive analyses about the relationship between DOM and bacteria. With the help of optical approaches, chromophoric DOM (CDOM) could be identified from the bulk DOM pool and used to indicate properties such as molecular weight (MW), aromaticity, etc. (Li and Hur, 2017). A part of CDOM that emits fluorescence exposed to ultraviolet light is named fluorescent DOM (FDOM) and could be further characterized as different components with fluorescence spectroscopy and parallel factor analysis (PARAFAC) (Stedmon and Bro, 2008). Meanwhile, the popularity of 16S rRNA gene sequencing and associated marker-based prediction tools make exploring bacterial community compositions and structure a handy and efficient work.
For bacterial community compositions, LaBrie et al. (2021) employed optical-derived properties of DOM to reveal that bacterial spatial clusters were better related with more complex DOM. Moreover, bacterial community compositions were reported to be associated with the aromaticity of DOM indicated by the spectral slope and showed depth-specific patterns in the Atlantic Ocean (Guerrero-Feijoó et al., 2017). A laboratory study confirmed the validity of optical properties in inferring DOM sources even after intense biological modification, and it has been reported that DOM sources had impacts on extracellular enzyme activities of bacterial communities (Freixa et al., 2016;Hansen et al., 2016). Besides, the lability or freshness of DOM was believed to be a key driver of bacterial metabolism (Loṕez-Urrutia and Morań, 2007;Nguyen et al., 2012). After high-labile DOM was taken up by microbes in the epipelagic zone, a small portion of DOM may survive the thorough remineralization with a certain degree of lability, which would be transported to the deep sea and consequently further utilized by deep-sea microbial communities (Carlson et al., 2004). The destiny of fresh DOM is of great importance for the global carbon cycle and it is basically decided by bacterial uptake (Puddu et al., 2003). The uptake rates depend not only on the state of the local bacterial community, but also on the chemical composition of organic substrates (Moriarty and Bell, 1993), which reflect linkages between DOM and bacteria. Apart from modifying the bacteria community with chemical compounds traded between DOM and bacteria, interspecies relationships and co-occurrence patterns within the bacterial community may also change due to the metabolic network variance stimulated by DOM (Bai et al., 2020). Zhou et al. (2021c) found that elevated bioavailable DOM induced by eutrophication enhanced the modularity of bacterial interspecies networks. Based on an incubation experiment, there were different dominant taxa in co-occurrence networks associated with different FDOM components (Li et al., 2020b). Also, the quantity of phytoplankton-derived DOM was reported to affect the bacterial niche partitioning (Sarmento et al., 2016).
Compared with water column DOM pools, biogeochemical processes of DOM pool in sediments have been rarely revealed. As the most reactive regions within sediments (Fox et al., 2018), pore waters were believed to be an indispensable intermediate of sedimentary carbon preservation and an important part for estimating the global carbon budget (Burdige, 2001;Chen and Hur, 2015). Some previous studies have reported porewater DOM as a net upward benthic efflux to overlying water (Wang et al., 2014;Burdige and Komada, 2015). Compared to the water column, the DOM pool in sediment pore waters showed a higher average molecular weight, and a higher humidification index (HIX) (Ziegelgruber et al., 2013;Chen and Hur, 2015). There were significant differences in fluorescent components and their distribution between sediment pore water and water column (Chen and Hur, 2015;Li et al., 2015;Li et al., 2020a). While abovementioned studies have reported remarkable findings, DOM-bacteria processes in sediments remain poorly resolved and studied in this field.
Among the East Asian margin seas, the East China Sea (ECS) receives a great amount of terrestrial organic matter from rivers (Yang et al., 2011). Coupled with the complex current system and strong ocean-land interactions, ECS provides a broad gradient for our study. Considering the fact that porewater DOM is the prevailing substrate for bacteria in marine sediments, we hypothesized that the bacterial communities might be significantly controlled by the varied freshness degrees or terrestrial contributions quality of the porewater DOM in the ECS. The objectives of our study thus were to (1) clarify variations in optical-derived properties of DOM and the bacterial community and (2) elucidate the linkage between these variations in coastal sediments.

Study Area and Sample Collection
A total of 42 samples from 10 section were collected from the coastal area of ECS with a grab sampler during a cruise from August to September 2015 ( Figure 1). We collected surface sediment samples (ca. 10cm) from 8-66 m depth using a CTD (Sea-bird, USA). Salinity and dissolved oxygen of overlying waters varied between 17.65‰ -34.42‰ and 1.063 -6.098 mg/L, respectively. Samples were stored in the refrigerator of -20°C and protected from light until further processing. Total organic carbon (TOC), carbon isotope, chlorin content and glycerol dialkyl glycerol tetraethers (GDGTs) content measurement were performed as described in our previous studies (Li et al., 2020a).

DNA Extraction, Sequencing and Sequence Analysis
Following the manufacturer's protocol, 0.5 g of sediment samples from each site were used to extract DNA with FastDNA Spin Kit for Soil (MP Biomedicals, USA). The V3-V4 regions of bacterial 16S rRNA genes were amplified with primers 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTAC NVGGGTWTCTAAT-3′). The pair-end sequencing was performed on the Illlumina MiSeq platform by Shanghai Majorbio Bio-pharm Technology Co., Ltd (Shanghai, China). 'Quantitative Insights into Microbial Ecology' pipeline (QIIME2 2019.7.0 release) (Bolyen et al., 2019) was employed for sequence analysis together with R. DADA2 (Callahan et al., 2016) plugin in QIIME2 was applied to denoise sequencing reads and produce high-resolution amplicon sequence variants (ASVs) for subsequent analysis. Taxonomic analysis was performed with a classify-sklearn classifier trained against SILVA 16S rRNA gene reference database (Quast et al., 2012). Following quality filtering, a total of 650,936 reads from 16,741 ASVs were collected with an average of 15,498.48 reads per sample. To improve downstream diversity and composition analyses, sequences were rarefied to the lowest number of sequences per sample (n = 11,455 sequences). Thus, there were 16,703 ASVs removed (0.23% of the total number of ASVs). Additionally, richness estimators in alpha-diversity indices including Chao1, the observed richness (sobs), ace and phylogenetic diversity (pd) were also calculated. Slopes of rarefaction curves based on sobs gradually decreased increasing sequencing depth, and finally tended flatten out ( Figure S1), indicating that the sequencing depth of samples was adequate. Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2, v2.2.0-b) were used to predict the functional potentials of the bacterial communities. Raw sequence data are available from the NCBI Sequence Read Archive under accession number PRJNA761364.

DOM Optical Characterization and Analysis
About 20 g of sediments were used to extract pore water samples by centrifugation at 4000 rpm for 20 min. After centrifugation, pore waters collected were filtered with 0.2-mm polyethersulfone membranes (Millipore Millex-GP, USA). CDOM absorbance spectra were obtained by UV2600 UV-Vis spectrophotometer (Shimadzu, Japan). The scanning range was 200-800 nm with 0 nm increments in a 1-cm cuvette and Milli-Q water (Milli-Q Integral 10, Millipore, France) was used as blank control. The CDOM absorption coefficient at the wavelength of 254 nm and 340 nm (a 254 and a 340 ) were used to estimate the CDOM concentration, with a 254 refers to information of relatively simple compounds and a 340 refer to that of complex molecules (Stedmon and Nelson, 2015). CDOM spectral slope from 275 nm to 295 nm (S 275-295 ) and spectral slope ratio (SR) were calculated as suggested by (Helms et al., 2008) and both showed negative correlations with CDOM MW. We employed Excitationemission matrices (EEMs, Figure S2) together with PARAFAC FIGURE 1 | Map of the study area with 42 sampling sites and the spatial distribution of salinity of bottom water in the East China Sea (Ocean Data View v5.2.1). (Stedmon and Bro, 2008) to characterize FDOM as five components (C1-C5, Figure S3). Detailed information about FDOM characterization, HIX, fluorescence index (FI) are available from our previous study (Li et al., 2020a). Biological index (BIX) was calculated at the excitation wavelength of 310nm, by dividing the fluorescence intensity at the emission wavelength of 380nm by that at 430nm. BIX was a freshness index, with higher values indicating a greater proportion of recently autochthonous contribution fresh DOM and samples with BIX > 1 being considered to be marine autochthonous (Huguet et al., 2009). We categorized 42 study sites into two groups: Group Fresh (n = 26) with BIX > 1 and Group Non-Fresh (n = 16) with BIX < 1.

Statistical Analysis
To investigate the distance among samples or between different groups based on their bacterial community, non-metric multidimensional scaling ordination analysis (NMDS; Bray-Curtis distance) was conducted using the function 'metaMDS' within the 'vegan' package in R. Also, the significance of differences between Group N and Group F based on NMDS was tested with analysis of similarities (ANOSIM) with the function 'anosim' within the 'vegan' package in R (Oksanen et al., 2013). Levins' niche breadth index for each bacterial taxon (Levins, 2020) and niche breadth at the community level (Bcom) were calculated as proposed by Jiao et al. (2020) Before analyzing the relationship between environmental variables and bacterial communities, variables with variance inflation factors larger than 10 were removed to avoid collinearity. Then, detrended correspondence analysis (DCA) was carried out to decide whether redundancy analysis (RDA) or canonical correlation analysis (CCA) should be adopted for downstream analysis. In this study, RDA was selected in all analyses since the longest lengths of the gradient in DCA were less than 3. Spearman correlation coefficients and their corresponding p values were calculated using the 'Hmisc' package in R and visualized using the 'corrplot' package in R (Wei et al., 2017;Harrell Jr and Harrell Jr, 2019). For relationships between environmental factors and bacterial assemblages, Mantel and partial Mantel tests were performed with 'vegan' (Oksanen et al., 2013). Due to the data requirements during calculations, some analysis involved with d 13 C employed the negative value of d 13 C (labeled as -d 13 C). To explore the cooccurrence network of bacteria, pairwise Spearman correlations were calculated and networks were visualized with Gephi v0.9.2. Moreover, topological roles of nodes within the networks were decided by within-module degree (Zi) and Participation Coefficient (Pi) calculated with 'GAINT' plugin in Cytoscape v3.8.2. A linear discriminant analysis Effect Size (LEfSe) analysis was used to identify taxa that were significantly different within two groups.

General Pattern of Optical Parameters
A suite of optical parameters of the pore water from the sediments were derived by fluorescence and absorption spectra.
The absorbance-derived quantitative parameters a 254 and a 340 showed increasing trends (r = 0.498 and r = 0.431, respectively; p < 0.01; Figure S4) from coastal area to open ocean, proclaiming higher concentration of CDOM in more off-shore sites. In brief, five components were identified with EEMs-PARAFAC, two of which were humic-like and three are protein-like. For the remainder of this article, we will refer to C1-C5 as the component itself or its fluorescence intensity of the five components: C1_r-C5_r, Humic_r and Protein_r as the proportions of total intensity (i.e., relative intensity) of the five components, humic-like components and protein-like components, respectively. The fluorescencederived parameters HIX (0.32 -3.61) and FI (1.07 -2.08) indicated that this FDOM pool consisted of both terrestrial and marine end-members. Freshness index BIX ranged from 0.42 to 1.37, with an average value of 1.05. In total, 26 pore water samples had BIX values greater than 1, indicating that they were mainly of autochthonous origin and contained fresher DOM (Huguet et al., 2009). The other 16 samples, however, had BIX values lower than 1 and may contain more non-fresh DOM ( Figure S5).
Heatmaps at the phylum and genus level ( Figures S6A, B) were constructed to illustrate the similarity between samples or taxa. At the phylum level, Proteobacteria was the most distant one and it occupied a single cluster. Dominant taxa Chloroflexi was clustered together with Desulfobacterota, Actinobacteriota and Acidobacteriota. For genera, four clusters were identified without any overt patterns. On the basis of both levels, the cluster analysis indicated no tendencies for different samples to be potentially grouped according to the relative distance from the land. According to qPCR analysis, samples from ECS contained (1.3 ± 1.5) × 10 10 bacterial copies per gram sediment, ranging from 7.1× 10 10 (copies/g sediment) at site A8_2 to 4.9× 10 8 (copies/g sediment) at site A3_2. No clear spatial variation of qPCR results was reported. As suggested by Lou et al. (2018), the abundance of each taxon in the soil bacterial community could be estimated by multiplying total bacterial quantities (the copy number of total bacteria from qPCR) and the corresponding relative abundance from high-throughput sequencing. Among 42 sites, Proteobacteria was the most dominant phylum with 3.8× 10 8 (copies/g sediment) on average and a yet no-rank genus in Order Actinomarinales the most dominant genus with 8.6× 10 8 (copies/g sediment) on average. NMDS analysis showed that 42 sites could not be grouped as a separated community with many sharing sites between two clusters ( Figures 3A, B).

Environmental Variables and Their Relationships With the Bacterial Community
The VIF analysis retained 18 variables for further analysis. RDA was carried out on these 18 environmental variables and bacteria relative abundances from the phylum to genus level, and then factors with p < 0.05 in RDA were applied to the Mantel test ( Table 1), suggesting that relative abundances of bacterial community compositions at all levels were positively correlated with their corresponding factor matrix determined after processes above. Unlike relative abundances, the bacterial abundance did not show any significant correlation with the corresponding factor matrix at all levels. It seemed that relative abundances of bacterial community could better reflect the effects of environmental variables.
To further investigate the relationship between bacteria and environmental variables, a correlation heatmap (Figure 4) was potted between 11 factors and the bacterial relative abundance of the 30 most abundant orders, based on the number of significant factors. Obvious correlations were noticed universally between bacterial relative abundances and factors that survived the elimination process. Among those factors, longitude and latitude reflect the location of the sampling site and therein the tradeoff between terrestrial and marine contributions contribution to sedimentary organic matter. Crenarchaeol and branched GDGTs (bGDGTs) together determined the proxy of the Branched and Isoprenoid Tetraether index (BIT) that correlated with the relative proportion of terrestrial input, while chlorin served as a measure of sedimentary primary productivity in the study area. The range of d 13 C (from -25.95‰ to -23.17‰) suggested a mixing source for porewater DOM while TOC decreased from nearshore to offshore. As a result, both GDGTs and d 13 C illustrated the effect of DOM sources on bacterial communities. These seven variables could refer to DOM sources and are therefore able to affect the properties of DOM. The remaining variables, C5 and C4_r, represented respectively the intensity and relative intensity of  FDOM protein-like, while BIX and FI measured freshness and differential sources of FDOM respectively. Alpha-diversity indices including chao1, sobs, ace and pd were calculated for bacterial richness at the genus level. They got range of 112-26.67 for chao1, 112-262 for sobs, 112-262.93 for ace and 0.9994-1 for pd. After the process as mentioned above, factors matrix corresponding to richness matrix (i.e., chao1, sobs, ace and pd) for Mantel test included C5, bGDGTs and TOC. Mantel_r value between two matrixes was 0.2214 (p = 0.0067). This result suggested the relevance existing between bacterial richness and environmental factors. Given the nature of statistics that parts of data should be discarded to gain more important information (Stigler, 2016), details in the data will be inevitably lost more or less in the statistical process. To gain a deeper understanding of the relevance, we developed a correlation matrix for diversity indices and all 34 variations out of the consideration of too few factors retained after processes ( Figure 5). Still, quite a lot of correlations showed between richness and other variables. Notably, all four richness indices were positively correlated with BIT (r = 0.583 -0.531, p < 0.001) and Humic_r (r = 0.388 -0.42, p < 0.05), and negatively correlated with Protein_r (r = -0.42 --0.388, p < 0.05), suggesting higher bacterial richness in the sites with higher terrestrial input and higher humic-like proportions. According to the rank abundance distributions obtained at the genus level ( Figure S7), samples with BIT > 0.1 generally had greater ranks than those with BIT < 0.1, further indicating that sites with higher terrestrial contributions harbored higher bacterial diversity. It was also notable that the community-level (or sample-level) niche breadth B com showed a significant negative correlation with humic-like proportion (r = -0.4252, p < 0.01; Figure S8).

Bacterial Co-Occurrence Networks and Comparative Analysis Between Two Groups
As indicated above, DOM properties had profound effects on composition and diversity of the bacterial community. All 42 samples were divided into two groups (Group F with BIX ≥ 1, otherwise Group N) to explore the relationship between DOM bioavailability and bacterial communities. Ordination analysis (MNDS) was performed from the level of the phylum to genus ( Figures 6A, C) to validate of the grouping. Both relative abundances and abundances at all five levels varied significantly between two groups, with higher r values and lower p values for relative abundances at all five levels than that for abundances. In concordance with the dominant marine source of DOM in this area, most of the BIX values were in the highly fresh range of 0.8 to 1.2 ( Figure S5). To avoid pseudopositive results from the grouping procedure (i.e., setting the breakpoint as BIX = 1), ordination analysis was applied again with samples with BIX values between 0.8 to 1.2 ignored. From Figures 6B, D, bacterial community compositions in terms of relative abundances were still significantly different between two groups with higher r values, while the difference in abundance was no longer significant. This result confirmed the response of bacterial relative abundance to DOM with different freshness as opposed to abundances. Nevertheless, total bacterial biomass estimated by qPCR in samples of Group F were significantly higher than that of Group N (p = 0.02, Wilcoxon rank-sum test), in part due to more bioavailable substrates, although no significant differences were found in the richness indices between Group F and N (p > 0.05 for four richness indices, Wilcoxon rank-sum test).
Then we performed the subsequent network analysis with relative abundances at the genus level to identify the cooccurrence patterns in two groups. It was evident Group N possessed a more complex and denser network than Group F did ( Figures 7A, B). There were 583 and 544 nodes within the networks of Group N and Group F, respectively. The edges and density of Group N were 3.4 times and 3 times those of Group F, respectively. Significantly higher degrees or edges (p < 0.001, Wilcoxon rank-sum test) connected to each node appeared within the network of Group N than Group F. Based on the value of within-module degree Zi and participation coefficient Pi ( Figures 7C, D), all taxa within two networks belonged to either peripherals or connectors and these two categories were believed to represent specialists and generalists, respectively (Deng et al., 2012). There were 83.6% of nodes belonged to generalists and 15.3% to specialists for Group F, while numbers for Group N were 76.5% and 23% belonged to generalists and specialists, respectively. Results from LEfSe analysis revealed 213 Bacterial community compositions shown as "taxonomic level + r or a" in the first column means relative abundance (r) or abundance (a) at the specific level. BIX, biological index; FI, fluorescence index; bGDGTs, branched glycerol dialkyl glycerol tetraethers; TOC, total organic carbon; a 254 , absorption coefficient at the wavelength of 254 nm; C1-C5, fluorescence intensity in the Raman Unit of C1-C5; C1_r-C5_r, proportions of total intensity (i.e., relative intensity) of C1-C5. Significant factors are shown within brackets in the second column. Significance labels: **p < 0.01, ns (not significant) p > 0.05.
significantly differential bacterial features (LDA score > 2, p < 0.05) from phylum to genus levels between two groups. According to the ranking of the LDA score from large to small, the top 20 taxa were chosen to represent the impact on the difference between groups (Figure 8). Compared with Group F, Group N possessed lower relative abundances of Actinobacteriota, Bacteroidota, Alphaproteobacteria, Desulfobulbia, and Gammaproteobacteria and higher of Anaerolineae and Desulfobacteria.

Effects of DOM From Different Sources on Variation in Bacterial Richness and Abundance
DOM sources have been considered important in affecting the bacterial community composition and species diversity (Nagata, 2008;Landa et al., 2014). As our results showed, the relative abundances of several bacterial orders were strongly correlated with FDOM components including the marine autochthonous components C5 (Jamieson et al., 2014), as well as some FDOMderived parameters such as BIX and FI. Also, both the relative abundances of several bacterial orders and the richness of the bacterial communities were closely correlated with variables indirectly depicting DOM sources, for example, the location of sites (relative distance from land), GDGTs-derived terrestrial proxy (BIT) and chlorin (Harris et al., 1996;Pitcher et al., 2011) representing sedimentary productivity in the sediments of the East China Sea. In line with previous studies (Landa et al., 2014;Li et al., 2020b), our results suggestthat the variations in DOM sources had significantly impact on the bacterial communities. Species richness indicators, such as chao1, ace, sobs and pd, were positively bound up with C2_r ( Figure 5), a reported terrestrial proxy (Li et al., 2020a). As a result of terrestrial vegetation diversity and regional anthropogenic activities, DOM from terrestrial environments is more chemically heterogeneous, which may further induce changes in aquatic bacterial communities (Lobbes et al., 2000;Grieve and Marsden, 2001;Lennon and Pfaff, 2005). In this study, a high terrestrial contribution might result in a high DOM heterogeneity and consequently supported more diverse bacterial taxa. Since samples with higher C2_r had lower values of B com (r = -0.4436, p < 0.01; Figure S8), the amendment of terrestrial DOM pool may support a narrow niche breadth at the community level. In the case of less available resources, a wide community-level niche breadth would benefit the bacterial taxa with high metabolic flexibility to acquire enough resources (Qin, 2009;Jiao et al., 2020). In this study, samples with lower FIGURE 4 | Correlation analysis between significant factors and bacterial taxa (top 30) at the order level. Spearman's rank correlations are indicated by colors, with red for positive correlations and blue for negative ones. BIX, biological index; FI, fluorescence index; bGDGTs, branched glycerol dialkyl glycerol tetraethers; TOC, total organic carbon; C1-C5, fluorescence intensity in the Raman Unit of C1-C5; C1_r-C5_r, proportions of total intensity (i.e., relative intensity) of C1-C5. Significant factors are shown within brackets in the second column. Significance labels: ***p < 0.001, **p < 0.01, *p < 0.05. terrestrial input might have lower DOM heterogeneity and nutrition. Hence, bacteria had to adopt flexible strategies and shift to the utilizing what were available to them, while the niche breadth at the community level becomes wide with insufficient resources for bacterial utilization. It was reported that DOM heterogeneity could provide diverse compounds as resources to heterotrophs and may lead to higher diversity (Smith et al., 2018). Zhao et al. (2019) showed an increase in chemodiversity after adding DOM from new sources (virus-derived DOM), although other studies reported that DOM sources had minimal effects on bacterial communities Rink et al., 2007).

Effect of DOM Freshness on Shaping Bacterial Community Co-Occurrence Patterns in Marine Sediments
In general, a higher BIX value meant more freshly autochthonous DOM that is readily bioavailable (Sundh, 1992;Huguet et al., 2009;Romera-Castillo et al., 2011). In this study, the bacterial community compositions of those samples with lower BIX values (defined as non-freshness group, Group N) were distinctly different from those samples with higher BIX values (defined as freshness group, Group F) ( Figures 6A, C), indicating the preference of bacterial communities for substrates with varying levels of freshness. Except for the community composition, the biological interactions among different bacterial taxa were also considered to be affected by DOM in natural environments (Montoya et al., 2006;Osterholz et al., 2016). A network analysis was employed to detect co-occurrence patterns of sedimentary bacterial community in this study ( Figures 7A, B). A more complex and denser network was found in those non-freshness samples than in freshness samples, which suggests that the non-freshness DOM could enhance the bacterial species relationships in the sediments in the East China Sea. As DOM becomes less fresh and more refractory, a bacterial species that can adopt appropriate survival strategies of reducing energy loss would survive more easily One of the possible adaptive strategies may be to increase the exploitation of specialized resource due to genome streaming (Swan et al., 2013). This might be one possible reason for higher Red circles indicate sites of the Group F while blue triangles indicate sites of the Group N. Significance of difference between two groups based on NMDS was tested via ANOSIM. Significance labels: ***p < 0.001, **p < 0.01, *p < 0.05.
proportion of specialists in Group N than that in Group F. However, microbial transformation of DOM have been reported as a synergistic process involving different functional enzymes from different bacterial taxa and sometimes even utilization of simple DOM compounds required the participation of a wide phylogenetic spectrum of microbes (Cottrell and Kirchman, 2000;Zimmerman et al., 2013;Logue et al., 2016;Yang et al., 2020). To meet the DOM exploitation demand, bacteria in Group N should be more complementary to each other because they were more specialized than Group F, resulting in more interspecies connections within the network. Our results were agreeable with earlier studies indicating that fresh phytoplankton-derived DOM could be adverse to bacterial specialization (Sarmento et al., 2016). Another possible explanation for the difference between two groups may come from bacterial functional prediction. A metabolic pathway ko02024 responding to quorum sensing function in bacteria were significantly up-regulated (p < 0.01) in Group N ( Figure  S9). Under different physicochemical conditions, bacterial composition of functional groups would differ due to the selection of suitable metabolic pathways (Louca et al., 2016). Previous study had reported that some specialists of bacteria may produce specific DOM portions for quorum sensing (Hmelo et al., 2011). In agreement with previous studies (Zhou et al., 2021a;Zhou et al., 2021b), DOM freshness may influence bacterial metabolisms and thus be a driving force to shape the bacterial community co-occurrence. In turn, the complimentary utilization of N group DOM among different bacterial specialists would play a role in the marine cycle of DOM elements (McCarren et al., 2010).
There were two similar taxa, Desulfobulbia and Desulfobacteria, which are both enriched in Group F and Group N, respectively. They previously belonged to the phylum Gammaproteobacteria but now have been proposed as two new phyla (Young, 2007;Waite et al., 2020). Although both taxa were capable of reducing sulfate, Desulfobulbia had an additional assimilation pathway reducing nitrate to ammonia (i.e., denitrification) (Waite et al., 2020). It was reported that bacterial denitrification degree was correlated positively with the consumption of the B peak of FDOM (Barnes et al., 2012). In this study, C5 corresponding to peak B was identified as a tyrosine component and was considered to be biodegradable DOM (Coble, 1996;Jamieson et al., 2014). Thus, in Group F, C5 was expected to be much higher than Group N in their initial stage. However, Wilcox test results showed that A B D C FIGURE 7 | Bacterial co-occurrence networks of Group Fresh (A) and Group Non-fresh (B). Each node represents a genus, with larger size indicating higher degree. Edges were colored as their source nodes. Only Spearman's rank correlations |r| > 0.6 and corresponding p < 0.05 are shown as edges. In each network, nodes with same color belong to same module. Between two networks, modules with the same color in two networks did not necessarily mean that they are the same module. Based on participation coefficient and within-module degree, topographic roles of nodes are organized into two categories: peripherals (red) and connectors (blue) for Group Fresh (C) and Group Non-fresh (D). C5 in Group F was significantly lower than in Group N, which indicated a higher net consumption of C5 in the fresh F and perphaps indicated a high expression of denitrification by Desulfobulbia. Another differential taxon enriched in Group N was Anaerolineae, within which there were many taxa of filamentous bacteria in Group N. Filamentous bacteria were reported to show the high stability and low population change with varied resources load (Mielczarek et al., 2012), which may explain their enrichment in Group N owing their powerful competitiveness for varying DOM with different freshness degrees. As abundant taxa in Group F, Actinobacteriota were affected by both aromatic and protein-like substrates binding to humic compounds (Liu et al., 2019). In this study, DOM mainly originated from marine autochthonous production with a high BIX value. Thus, there were more humic-like carriers as well as more 'passenger' compounds like protein-like substrates, which may have supported the growth of Actinobacteriota. Alphaproteobacteria have been considered as a strong competitor for amino acids in marine ecosystems and they also subsist on almost all kinds of DOM compounds (Malmstrom et al., 2004;Kujawinski, 2011). Among the taxa, some clades like SAR11 became more important and produced more fresh DOM for other bacterial taxa (Poretsky et al., 2010). Upon exposure to marine autochthonous DOM, the abundance of both Bacteroidota and Gammaproteobacteria were elevated (Amaral et al., 2016). Both taxa were viewed as opportunistic organisms (Alonso and Pernthaler, 2006;Alonso-Saéz and Gasol, 2007) and this may explain why they were more likely to thrive in Group F with more bioavailable DOM. In this study, bacterial communities, either their co-occurrence patterns or dominant taxa, were associated with DOM freshness. These results implied that bacterial communities would probably change in their structure or composition as DOM converting from fresh to recalcitrant.

CONCLUSIONS
In this study, the bacterial communities were strongly linked to DOM optical-derived properties such as sources and freshness of DOM. Variations in the DOM properties and the bacterial communities indicated mutual responses to each other. Bacterial community compositions either as bacterial assemblages or individual bacteria taxa had significant correlations with DOM properties and the relative bacterial abundance responded to those properties more intensely than their abundance. Higher values of alpha-diversity were found in sediments with higher terrestrial contributions, which could be explain by DOM heterogeneity or nutrition within the allochthonous matter derived from variable ecosystems on land. It is possible that the high richness may enhance the complexity of the interspecies interaction since more taxa participated in the co-occurrence network. Although there were no significant differences in bacterial richness between Group F and Group N, bacterial biomass was higher in Group F due to more bioavailable substrates provided by higher freshness DOM. Network analysis identified distinct cooccurrence patterns between two groups. By comparison, bacterial taxa in Group N were connected much more frequently and had more specialists, which might be an adaptive strategy under low-resource conditions to reduce energy loss by elevating cooperation and specialization during the utilization of DOM. Differential taxa in two groups showed physiological characteristics to adopt their local environmental conditions. Our findings demonstrate that changes in DOM composition would not only modify bacterial community compositions but also enforce interspecies interactions within the bacteria community to respond to the increasing recalcitrance of DOM. Furthermore, this study also provides a FIGURE 8 | LEfSe analysis revealed the taxa that are significantly different between two groups base on the relative abundance. According to the ranking of LDA score from large to small, the top 20 taxa are shown in the figure with different colors and styles. framework to clarify the relationships between DOM and bacterial communities in marine sediments.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the NCBI repository (https://www.ncbi.nlm.nih.gov/), accession number PRJNA761364.