Fungal-Bacterial Cooccurrence Patterns Differ between Arbuscular Mycorrhizal Fungi and Nonmycorrhizal Fungi across Soil Niches

Soils near living and decomposing roots form distinct niches that promote microorganisms with distinctive environmental preferences and interactions. Yet few studies have assessed the community-level cooccurrence of bacteria and fungi in these soil niches as plant roots grow and senesce.

T he biological communities and biogeochemical cycling in soil near roots versus distant from roots reflect distinct ecological niches. Rhizosphere soil (directly next to living roots) receives 10 to 40% of plant photosynthetic products as rhizodeposits (1) that include low-molecular-weight compounds easily consumed by microorganisms, metabolites related to plant defense and stress tolerance, and root debris resulting from death and sloughing off of root structures. These rhizodeposits frame the formation of a rhizomicrobiome occupying the rhizosphere niche (2). In soil distant from and less influenced by roots, litter-derived older organic materials become major resources. The decomposition of macromolecular organic compounds and mineralization of nutrients by saprotrophic organisms in bulk soil containing root litter (the detritusphere) are critical processes in terrestrial nutrient cycling (3). The distinct physicochemical environments within the rhizosphere and the detritusphere influence how associated soil microbial communities assemble and function. Despite rich research on the microbial composition of these soil niches, we know little about the covariation patterns of bacteria and fungi in these soil environments. It has been suggested that many bacteria depend on fungal hydrolysis of recalcitrant plant materials such as lignin and cellulose, forming intermediate products that are primary bacterial substrates in soil (4). In return, bacteria contribute to the resources needed by fungal decomposition processes by providing nutritional benefits, for example, through in situ N fixation (5). Exploring fungal-bacterial cooccurrence patterns in rhizosphere and detritusphere niches can identify evidence of niche sharing and potential interactions that collectively affect carbon (C) cycling and nutrient availability to plants.
Bacteria and fungi are the most diverse groups of soil inhabitants, and their versatile physiologies help regulate many soil processes (6). Cooperatively, saprotrophic fungi and bacteria can increase the efficiency of litter decomposition (4,5,7). Fungal hyphae can connect patches in unsaturated porous soil, acting as transport vectors or "fungal highways" for bacterial dispersal and access to substrates (8)(9)(10)(11). Yet bacteria also consume fungal exudates and biomass released after hyphal damage or death and actively attack living hyphae through mycophagy, promoting the transfer of C from biotic forms to CO 2 , dissolved organic C, and mineral-associated C (12)(13)(14). In addition, antagonistic interactions between certain bacterial and fungal pathogens can promote plant pathogen resistance (15). As individual groups of organisms, fungi and bacteria are known to dynamically compete for space, nutrients, and essential minerals (16)(17)(18), especially under nutrient limitation (19). In aggregate, the location of these direct fungal-bacterial interactions is termed the hyphosphere (20). The hyphosphere is a hot spot for bacterial-fungal interactions with strong effects on bacterial niche development (13), the consequences of which are only beginning to be understood. Studies suggest that mycorrhizal and saprotrophic fungi modulate bacterial community composition (21)(22)(23)(24)(25), thereby having differential effects on C in rhizosphere, hyphosphere, and detritusphere soil habitats (26). Yet most of the existing bacterialfungal interaction research has been carried out under highly controlled laboratory conditions (i.e., on artificial media), and much remains unclear about the distribution and dynamics of bacterial and fungal communities in natural root-soil-microbe systems.
The associations between arbuscular mycorrhizal fungi (AMF) and bacteria are particularly important to plant nutrient acquisition and the redistribution of the plant photosynthate into belowground C pools. AMF are strict biotrophs that depend on plantderived C via the exchange of plant-fixed C for nutrients that AMF absorb from soil (23) and are estimated to colonize 71% of land plants (27). AMF-plant interactions occur in the rhizosphere during the establishment of symbiosis, but once established, they become conduits that redistribute C from the roots to areas beyond the rhizosphere (i.e., the "bulk" soil). Bacteria have been shown to inhabit the hyphosphere, the area around mycorrhizal hyphae (28,29). Classical examples include helper bacteria that facilitate the establishment of plant-mycorrhizal fungus symbiosis (30), plantgrowth-promoting rhizobacteria (PGPR) that work synergistically with AMF to enhance plant growth (31), and mycophagous bacteria that feed on live hyphae (32). While insightful research has been conducted on the influence of AMF-bacterial associations on soil and plants, little is known about the location of these hyphal-bacterial associations in soil and when they occur during the growth phases of plants.
Ecological networks have been increasingly reported as approximations of the complex relationships among members of diverse microbiomes (33)(34)(35), where nodes represent members and links depict characteristics of their relationships (directionality and weight, etc.). In plant and animal networks, links are based on prior knowledge of biotic relationships (e.g., food webs, pollination networks, and pathogenic networks). In macroorganisms, network behavioral studies have greatly expanded our understanding of community dynamics, stability, and evolution (36)(37)(38). In microbial ecology, however, networks are usually inferred based on abundance correlations between molecular markers (34), where the cooccurrence of microbial taxa can result from environmental filtering, resulting in niche partitioning or sharing (39) that is exceedingly difficult to document in the natural microbial world. Although these correlation-based networks have inevitable limitations in their capacity to distinguish different mechanisms underlying the links (39,40), they are thought to capture predicted and novel patterns in microbial community organization, including bacterial competition following substrate injection into groundwater (41), fungal functional guild partitioning in agriculture soils (42), plant-microbe associations (43), disease suppression against plant pathogens (44), increases in soil bacterium cooccurrence in response to climate change (45), and microbial metabolic dependence (46). In a previous related study (47), network analysis revealed highly complex cooccurrence patterns among bacteria living in the rhizosphere during the progression of plant growth. These patterns contrasted with the more simple and static networks in nonrhizosphere soil and suggest that niche sharing is critical to rhizosphere community assembly. In that same study, a keystone species identified based on network topology was found to harbor quorum sensing ability, suggesting that these networks captured potentially interacting bacteria. Networks based on cooccurrence patterns can serve as an important starting point for exploring potential interactions in complex soil microbiomes by identifying possible biotic interactions that are worthy of further experimental validation (39). Here, we explore cross-kingdom networks among bacteria and fungi to provide a hypothesisgenerating procedure to understand the interactions between these two consequential groups of soil organisms over space and time; the identified fungal-bacterial links can be further examined to assess possible interactions or niche sharing.
In this study, we investigated the cooccurrence patterns of bacterial and fungal communities over space and time, using amplicon sequences of bacterial 16S and fungal internal transcribed spacer (ITS) rRNA genes, in soil samples from a greenhouse experiment that has previously yielded information on bacterial succession and cooccurrence patterns (47,48) as well as potential trophic interactions revealed by stableisotope-informed genome-resolved metagenomics (49). This experiment used soil from a California annual grassland and the annual grass Avena fatua, which is common and naturalized to the area from which the soil was collected. Growth conditions approximated those of two consecutive growing seasons separated by a summer drydown period characteristic of the Mediterranean climate. Both rhizosphere and bulk soils were harvested in a multipoint time series and analyzed to compare soil microbial communities adjacent to and away from roots, at multiple plant growth stages, and in response to the effect of summer dry-down. We constructed cooccurrence networks using a random matrix theory (RMT)-based approach (45,50) and separately explored the cooccurrences of soil bacteria with both mycorrhizal and nonmycorrhizal fungi at the levels of communities, populations, and taxa in order to understand how soil niche and different life histories of organisms influence the dynamics of fungal-bacterial cooccurrence over space and time. We then used stable-isotope-enabled genome-centric metagenomics to investigate the links between bacterial and fungal taxa at different plant growth stages in soils from rhizosphere, bulk, and detritusphere zones. We sought to address the following: (i) What are the patterns of bacterial-fungal cooccurrence in rhizosphere versus bulk/detritusphere soil and over time through a plant growing season? (ii) How do cooccurrence patterns differ between bacteria and mycorrhizal versus nonmycorrhizal fungi? (iii) What taxa of bacteria and fungi cooccur, and what are the potential mechanisms behind these patterns? Compared to previous studies of fungal-bacterial interactions (primarily conducted under controlled conditions), our analysis is unique in its use of live soil and native fungal and bacterial communities.

RESULTS
Overall network topology. Using 16S rRNA and ITS amplicons (see Materials and Methods for basic sequencing results), we created separate nonmycorrhizal (including saprotrophic, pathogenic, endophytic, and unclassifiable) and mycorrhizal (AMF) fungal-bacterial networks for five time points in each of two growing seasons, for both rhizosphere and bulk (or bulk and decaying root litter) soils (16 biological replicates were used to construct each network).
The nonmycorrhizal fungal-bacterial networks ( Fig. 1; see also Table S1 in the supplemental material) have 310 to 783 nodes (operational taxonomic units [OTUs] whose relative abundance covaries with at least one other OTU) and 368 to 790 links (pairwise correlations of nodes), with an average degree (average number of links per node) ranging from 1.5 to 2.9. On average, 35% of the nodes are fungal, representing 12.5% of all detected nonmycorrhizal ITS OTUs at each sampling time point and soil location. On average, 85% of the links in the networks are positive (positive correlation between node relative abundances). There is no discernible pattern of node taxonomy changes over time at the phylum level (Fig. S2f).
A total of 156 ITS OTUs were identified as Glomeromycota, covering a wide range of the AMF phylogeny. Among the 26 AMF OTUs included in network construction (those with frequencies of $10 in 16 replicates), all of which were detected in both rhizosphere and bulk soils, 22 had positive or negative links with bacterial OTUs. Notably, AMF-bacterial links were detected in bulk soils in both seasons but were detected in the rhizosphere only in season 2 ( Fig. 2; Table S1). These AMF-based networks are relatively small, with 2 to 51 nodes and 1 to 41 links. Averaged across sampling times and locations, only 4% of the detected AMF OTUs were present in these networks, with averages of 6.3% for the bulk soil and only 1.1% for the rhizosphere. A total of 87% of the AMF network links are positive, and fungi make up 36% of all the nodes, similar to the nonmycorrhizal fungal-bacterial networks.
Sampling location shapes network complexity and succession. Sampling location had a strong impact on network topology, resulting in distinctly structured networks in rhizosphere and bulk soils. Although nonmycorrhizal fungal-bacterial networks ( Fig. 1; Tables S1 and S2) have similar numbers of links between rhizosphere and bulk soils, the rhizosphere networks contained 24% fewer (F 1,12 = 24.0; P , 0.001) nodes on average. Accordingly, the rhizosphere networks have a 30% higher average degree (F 1,12 = 34.0; P , 0.001) and 85% higher connectance (F 1,12 = 104.4; P , 0.001) than the bulk soil networks. Rhizosphere network modules are on average 47% larger (F 1,12 = 16.72; P = 0.002) than the bulk soil networks. In addition, the proportion of positive links is 13% higher (F 1,12 = 35.6; P , 0.001) in the rhizosphere, suggesting more cooccurrence than coexclusion of bacteria and nonmycorrhizal fungi.
For the nonmycorrhizal fungal-bacterial networks, the rhizosphere networks illustrate a far stronger successional pattern with plant growth than the bulk soils. Sampling week had a significant effect on the network size (F 1,12 = 6.2; P = 0.028), average degree (F 1,12 = 10.9; P = 0.006), average module size (F 1,12 = 5.99; P = 0.031), connectance (F 1,12 = 46.9; P , 0.001), and proportion of positive links (F 1,12 = 8.5; P = 0.013) in FIG 2 The succession of AMF-bacterial networks over different plant growth stages in two seasons. Node shapes represent AMF (circles), bacteria (triangles), or archaea (squares), and node colors differentiate phyla for bacteria/archaea and classes for AMF. The OTU identifier for each AMF node is marked and corresponds to those in Table S5 in the supplemental material. The bacterial nodes that are mentioned in the text are also marked with their OTU identifiers.
Fungal-Bacterial Cooccurrence across Soil Niches ® these networks (Table S2). General linear models (Fig. S2) show that the average degree (adjusted r 2 = 0.81; P = 0.024), connectance (r 2 = 0.78; P = 0.030), module size (r 2 = 0.77; P = 0.031), and proportion of positive links (r 2 = 0.72; P = 0.045) increased strongly and significantly over time for the rhizosphere networks in season 1. For the season 2 rhizosphere networks, only connectance showed a significant linear increase with time (r 2 = 0.90; P = 0.008). Network size (r 2 = 0.67; P = 0.056) decreases with time were marginally significant but with a meaningful effect size. In the rhizosphere networks, the average degree and average module size increased until week 9 but reversed in week 12 (plant senescence) in both seasons, suggesting that plant senescence decreased rhizosphere microbial network complexity. For both seasons, bulk soil networks remained static, with no topological parameter showing any significant change over time (Fig. S2).
Accounting for time points with no significant cooccurrence, the rhizosphere AMFbacterial networks ( Fig. 2; Tables S1 and S2) have on average 71% fewer nodes (F 1,12 = 7.7; P = 0.017) and 75% fewer links (F 1,12 = 8.1; P = 0.015) than the corresponding bulk soil networks. The differences in AMF-bacterial network size and complexity between rhizosphere and bulk soils are evident in Fig. 2.
For AMF-bacterial networks, due to a lack of data points in the rhizosphere networks, we did not detect any statistically significant change in network topological parameters over time with plant growth (Table S2). However, the bulk soil networks from season 2 (which contained decaying roots from the previous season) were significantly more complex than the season 1 bulk soils, where there was no root detritus. The numbers of nodes and links are also significantly different between these rhizosphere and bulk soil networks, with greater complexity in the bulk soils. In the season 2 bulk soils (which included root litter), we measured AMF-bacterial cooccurrence in the soil prior to A. fatua being planted (23 nodes and 15 links) ( Fig. 2; Table S1).
Contrasting patterns of bacterial cooccurrence with AMF versus nonmycorrhizal fungi. Bacterial cooccurrence with both AMF and nonmycorrhizal fungi was influenced by the presence of roots and plant growth timing but manifested quite differently between the AMF and nonmycorrhizal fungi. First, unlike nonmycorrhizal fungal-bacterial networks, AMF-bacterial cooccurrence was observed frequently and with increased complexity over time in bulk soil relative to rhizosphere soil, suggesting that although AMF are strict biotrophs, they tend to cooccur with bacteria away from roots. Second, growth season had a stronger impact on the cooccurrence of bacteria with AMF than with nonmycorrhizal fungi. There were over 5 times as many nodes (F 1,12 = 10.8; P = 0.007) and 6 times as many links (F 1,12 = 10.6; P = 0.007) in the AMF-bacterial networks for season 2 compared to season 1 ( Fig. 2; Tables S1 and S2). For nonmycorrhizal fungal-bacterial networks, although their sizes were similar in the two seasons (F 1,12 = 0.47; P = 0.507) (Table S2), the season 1 networks showed a slightly but significantly higher average degree (F 1,12 = 6.6; P = 0.024), connectance (F 1,12 = 12.8; P = 0.004), and proportion of positive links (F 1,12 = 5.9; P = 0.032) than those in season 2 (Table S2 and Fig. S2).
The complexity of AMF-bacterial networks over time tracked with the increasing richness of the detected AMF taxa, while the complexity of nonmycorrhizal fungal-bacterial networks followed an opposing trend to the declining pattern of nonmycorrhizal fungal richness (Fig. S3). More AMF OTUs were detected in bulk soil than in the rhizosphere; bulk soil AMF OTUs increased throughout the two seasons, regardless of summer dry-down. In the rhizosphere, the richness of nonmycorrhizal fungi decreased with plant growth, yet the network complexity indicated by average degree and connectance increased (Fig. S2). In addition, the summer dry-down brought the richness and level of complexity of the nonmycorrhizal fungal-bacterial networks to levels similar to those at week 0 in season 1.
Fungi comprise the majority of keystone species, and fungal nodes tend to occur in multiple networks. We identified keystone nodes (see Materials and Methods for criteria) in the fungal-bacterial networks (containing both nonmycorrhizal fungal-bacterial and AMF-bacterial links). Two to 18 module hubs were identified at each date and location, with a total of 172 module hubs belonging to 128 OTUs (Table S3). Among all the module hubs, 95% (163 out of 172) of the keystone nodes were fungi, mostly Ascomycota, while only 8 were bacteria. In contrast, for the composition of all the network nodes (including keystone nodes and those that were not), only 35% belonged to fungi. No connectors were identified in these networks.
A higher percentage of fungal OTUs than bacterial OTUs occurred in multiple networks. In the 18 nonmycorrhizal fungal-bacterial networks, over 36% of the 541 unique fungal OTUs occurred in more than 9 of these networks (Fig. S4b), while only 0.2% of the 2,825 unique bacterial OTUs occurred in more than half of the networks (Fig. S4a). Similarly, all except one bacterial node in AMF-bacterial networks were present in only a single network (Fig. S4c), while one AMF node (OTU 5348) occurred in 8 of the 13 AMF-bacterial networks (Fig. S4d).
Potential mechanisms underlying the observed fungal-bacterial cooccurrence as informed by genome-resolved metagenomics. In a separate analysis of a subset of these soil samples (49), we resolved multiple bacterial metagenome-assembled genomes (MAGs) from samples collected from weeks 0, 6, and 9 of this study. In five of these MAGs, we discovered the corresponding V4 region of the 16S rRNA gene that matched with high confidence (.99% similarity for at least 200 bp) to the representative sequences of the 16S rRNA OTUs that we used to construct our networks. These bacterial genomes include two Acidobacteria (Acidobacteria_68_21 and Acidobacteria_60_12), one Actinobacteria (Actinobacteria_70_13), and two Alphaproteobacteria (Rhizobium_leguminosarum_59_9 and Mesorhizobium_63_15) genomes. A total of 38 links, including 3 negative and 35 positive ones, connected bacterial nodes matching these genomes to fungal nodes. The linked bacterial and fungal nodes belonged to 8 bacterial OTUs and 31 fungal OTUs. Half of these links were in rhizosphere networks, and the other half were in bulk soil networks ( Fig. 3; Table S4).
The five MAGs corresponding to nodes in our networks harbored a number of genes that may have been involved in interorganismal interaction pathways (Table 1).
Based on their phylogeny and predicted ability to modulate plant hormones, Rhizobium_leguminosarum_59_9 and Mesorhizobium_63_15 may be plant symbionts. These two taxa were far more abundant in the rhizosphere than in bulk soils ( Fig. S5a and b) and were also highly 13 C enriched in the stable-isotope probing (SIP) analysis of the same samples (49). The 16S OTUs matching Rhizobium_leguminosarum_59_9 formed 21 positive and 2 negative links with nodes of Ascomycota, Chytridiomycota, and a few unclassified species in networks constructed for different plant growth stages. The OTU matching Mesorhizobium_63_15 formed five positive links with fungi (Table S4). Acidobacteria_68_21 could be a bulk soil specialist based on its high relative abundance in bulk soil samples (Fig. S5c) and the fact that it was not labeled by plant-derived 13 C. It harbors many defensive capabilities by carrying a large arsenal of biosynthetic gene clusters, a type IV secretion system, CRISPR-Cas systems, and, notably, insecticidal proteins, which could be used in defense or offense against other bacteria, nematodes, bacteriophages, and fungi. Its central metabolism is highly flexible, and its genome harbors a wide array of carbohydrate-active enzyme (CAZy) genes targeting different substrates. This MAG was the most abundant organism identified in the bulk soil in that related study (49). The OTU matching this genome had four positive links with fungal nodes (Table S4). The bacterial node matching Acidobacteria_60_12 formed only one negative link, with a Myxotrichaceae node (Table S4). Acidobacteria_60_12 exhibits antifungal potential by carrying 18 CAZy genes targeting fungal biomass compounds (49,51); a number of biosynthetic gene clusters, including one predicted to biosynthesize a phenazine-like compound, a documented antifungal compound (52); and antibiotic resistance modules and reactive oxygen species (ROS) protectants that could provide defense against fungi (49).
Fungal-Bacterial Cooccurrence across Soil Niches ® Finally, the Actinobacteria_70_13 genome carried a large number of CAZy CE14 genes, a chitin deacetylase, and some large genes likely used in adhesion to biological material (Table S4). The nodes matching this genome had five positive links with fungal nodes (Fig. 3).

DISCUSSION
By investigating fungal-bacterial networks in a native grassland soil, we found that different soil niches had distinctly different cooccurrence patterns for bacteria and AMF versus bacteria with nonmycorrhizal fungi. AMF were more frequently linked to bacteria in soils .2 to 4 mm from a root and in the presence of the previous season's root litter, while nonmycorrhizal fungi formed more complex networks with bacteria in the rhizosphere. Cooccurrence of species can result from combined effects of environmental filtering related to niche sharing, biotic interactions, as well as dispersal limitation (53). However, it is likely that dispersal limitation is a less important driver for cooccurrence in our experiment because the soil was sieved and homogenized before planting. Our results indicate that distinct niches partitioned by the availability and forms of C substrates are likely fundamental drivers of how microorganisms of different life histories and ecologies cooccur and potentially interact in complex soil environments.
Cooccurrence between AMF and bacteria. Known to be strictly biotrophic organisms that live on carbon derived from living roots, AMF had both higher richness and  relative abundance in soils outside the rhizosphere, especially in the presence of root litter in our second growing season. This apparent contradiction may be explained by AMF's life form in soils. First, as roots take up nutrients, a depletion zone is formed in their immediate vicinity (1 to 3 mm) (54)(55)(56). After colonizing roots, AMF are driven to extend the majority of their hyphae beyond this root-adjacent depletion zone to reach for nutrients that roots cannot access. In fact, it has been shown that AMF greatly extend the phosphorus depletion zone around clover roots (57,58), demonstrating their search for nutrients away from roots; this could lead to a higher detectable AMF relative abundance in nonrhizosphere soils. Second, AMF spores are formed on hyphae outside roots, and previous evidence shows that the nucleic acid concentration is much higher in fungal spores than in vegetative hyphae (59)(60)(61)(62)(63)(64). Thus, a combination of spores and hyphae may lead to the detection of more AMF using DNA-based molecular markers in nonrhizosphere soils than near roots. Finally, the second season's root litter (derived from the first season's root growth) likely included some AMF-colonized roots. These may have served as inocula, leading to AMF's quick propagation in the second season and increased our ability to detect their presence. This explains the progressive increase of AMF richness and relative abundance over the two growth seasons despite the intervening summer dry-down as well as the high AMF diversity that we measured at week 0 in season 2, before seedlings had been planted. Our observations are consistent with another study that reports slightly higher richness of AMF in bulk soil than in the rhizosphere (65) in a sorghum cultivar, indicating that this pattern may persist across different plant species and soil types. Since only relative abundance is available for AMF versus nonmycorrhizal fungi in molecular marker-based studies, further tests such as hyphal and spore counts could be helpful to confirm the absolute abundance of AMF in rhizosphere versus bulk soil.
In our rhizosphere networks, a low AMF relative abundance, as mentioned above, together with a lack of common life history traits between AMF and bacteria, could have caused the underrepresentation of AMF-bacterial cooccurrence. Rhizodeposits, characterized by abundant low-molecular-weight C substrates, drive the formation of a rhizosphere niche shared by taxa that actively consume these substrates (2,66). Bacteria in the rhizosphere form complex networks with each other because of their shared niche space (47). AMF, however, depend on C provided from inside the root  (67), representing a closer association with plants that is fundamentally different from the relationship between plant and rhizosphere bacterial communities. Also, the processes by which AMF sense and actively colonize roots (67) differ from the recruitment of rhizodeposit-consuming microbial communities (68), although AMF's functional association with plants may overlap those of certain plant-growth-promoting rhizobacteria (PGPR), such as providing nutrients and alleviating drought stress (23,69,70). AMF perform these functions by transporting nutrients and water directly to internal root structures via hyphae (100), whereas PGPR function in soil outside roots to free nutrients as well as in the endodermis to regulate plant metabolites (70). In summary, although dependent on root C, AMF occupy a distinct ecological and spatial niche from rhizodeposit-dependent microorganisms.
Our detection of common AMF-bacterial cooccurrence in the detritusphere could also potentially result from functional associations between the two groups. AMF have no known saprotrophic capabilities but appear to be able to alter the rates of litter decomposition through their effects on saprotrophic organisms (7,23). AMF exudates can provide C to bacteria in the detritusphere and promote bacterial litter decomposition. AMF then take up nutrients freed during decomposition. For example, Glomus geosporum and Glomus constrictum hyphae have been associated with bacteria of the genera Cellvibrio, Chondromyces, Flexibacter, Lysobacter, and Pseudomonas, which are known biopolymer degraders (71). Also, the phosphate-solubilizing bacterium Rahnella aquatilis can be induced to mineralize phytate when grown with the fructose-exuding AMF Rhizophagus irregularis (72). In our networks, we found AMF nodes linked to bacteria from the same orders (Myxococcales, Sphingobacteriales, and Xanthomonadales) that include the above-mentioned biopolymer degraders, although the genera of these nodes are unidentified. Material exchanges between AMF and bacteria in soil, which can foster cooccurrence, tend to be more active at extending hyphal tips, which reach far beyond the rhizosphere (67). In a separate study, using the same soil planted with Avena barbata, we found that the AMF Rhizophagus intraradices transferred a significant amount of recent photosynthetic C across a 3-mm air gap to a separate soilcontaining hyphal chamber (A. Kakouridis Based on AMF transport of C to growing tips and where AMF and bacteria more frequently cooccur, we hypothesize that AMF-bacterial associations are most frequent at growing hyphal tips some distance from roots and that AMF hyphae can serve as connections between microbial activity hot spots in rhizosphere and detritusphere soils (73). This hypothesis may represent a critical characteristic of the hyphosphere and is worth further research. In addition, AMF spores are known to harbor intracellular bacterial microbiomes with two main groups of bacteria, Burkholderia-related Betaproteobacteria (74) and Mycoplasma-related Mollicutes (75). In our data, we found bacteria of the order Burkholderiales in networks linking to AMF OTUs, although further genetic information or culture experiments would be required to confirm these relationships in our system. Together, these observations imply that biological processes in the detritusphere could be important drivers of AMF-bacterial cooccurrence.
Cooccurrence between nonmycorrhizal fungi and bacteria. The highly complex cooccurrence network of nonmycorrhizal fungi with bacteria in the rhizosphere could result in large part from niche sharing, a common occurrence among bacterial networks (47). As is typical in natural soil (76)(77)(78), the majority of nonmycorrhizal fungi detected in this study were either saprotrophic fungi or taxa that have saprotrophic capabilities during part of their life cycle or under certain abiotic conditions. The availability of organic C is a major driver of the distribution of both bacteria and saprotrophic fungi in soil. The C-rich rhizosphere could foster niche sharing between rhizodeposit-dependent bacteria and a diverse range of fungi that actively assimilate root exudates (79), debris, and litter, resulting in a pattern of cooccurrence between the two. Organic C availability in bulk soil, on the other hand, tends to be more limited; thus, cooccurrence is less likely to happen than in the rhizosphere. In our second growing season, the presence of litter coincided with a slight increase in the complexity of cooccurrence between nonmycorrhizal fungi and bacteria, possibly because of added amounts of C substrates. But the increase in network complexity with the presence of litter is less obvious than the large and obvious difference in network structures in rhizosphere versus bulk soil, implying that rhizosphere environmental filtering is the main driver of nonmycorrhizal fungal-bacterial cooccurrence.
Given the importance of the rhizosphere environment in driving the cooccurrence of nonmycorrhizal fungi with bacteria, biochemical changes in the rhizosphere associated with plant senescence are likely the cause of the consistent decrease in network complexity at week 12. Approaching senescence, Avena roots exude less, and the composition of exudates changes from sugars and amino acids to stress-and senescencerelated molecules (66,80). Soon thereafter, roots gradually become litter. The resulting altered substrate type can vastly change the environment and microbial assemblage. Root exudates during vegetative growth, like sugars and amino acids, are consumed by versatile microbial taxa, which may promote niche sharing and the cooccurrence of microorganisms; the decline of these microbially preferred exudates and the proportional increase of lignocellulose-rich root litter substrates could become unsuitable for exudate-dependent bacteria and fungi and dissipate shared niches.
Evidence from our nearly complete bacterial genomes, assembled from metagenomes derived from the same experiment, supports the hypothesis that some positive fungal-bacterial links may be explained by shared ecological niches. We traced the movement of plant-derived 13 C into rhizosphere organisms' genomes; the Rhizobium_leguminosarum_59_9 MAG became highly labeled, with 58 atom% excess 13 C (49). This high level of isotope enrichment, in addition to plant symbiosis genetic systems present in the genome, suggests that the corresponding bacteria likely live in close association with the plant and actively take up exudates, thus representing an organism adapted to the rhizosphere. Here, we found that network nodes matching Rhizobium_leguminosarum_59_9 frequently cooccurred with phylogenetically diverse fungal nodes, including several module hubs. Many of the correlating fungal nodes, for example, Chalara heteroderae (81), Acremonium (82), Spizellomycetaceae (83), and Herpotrichiellaceae (83), come from clades with known plant pathogens or pathogens of plant pathogens that interact with plants frequently in the rhizosphere. The large number of fungal nodes that correlated with Rhizobium_leguminosarum_59_9 may represent a shared niche in the rhizosphere. In contrast, Acidobacteria_68_21 is likely a generalist adapted to living in bulk soil and decomposing a wide variety of compounds while protecting its own accumulated biomass against attack from other bacteria, microeukaryotes, and fungi. Based on the bacterium's large metabolic potential, it may live in microhabitats devoid of any reliable C substrate. The fungi linked to Acidobacteria_68_21 may inhabit the same niche and may also be specialized for low-resource microhabitats.
Some other functional potentials harbored by the bacterial genomes that we have sequenced from this soil suggest that some of the fungal-bacterial links involve bacteria with the capacity to break down and consume fungal biomass as well as potential defense mechanisms. Acidobacteria_60_12 was one of the two genomes with a negative link to a fungus. This organism may be inhibitory toward soil fungi; its genome has a large number of antibiotic resistance modules and reactive oxygen species (ROS) protectants (which could shelter the bacterium from fungal defenses) and encodes the production of antifungals that may inhibit or kill organisms, including fungi. The large number of fungal biomass decomposition genes, including chitinases and glucanases (49,51), suggests that this bacterium is involved in fungal carbon recycling, but it is unknown whether it actively parasitizes fungal hyphae.
Due to the harbored adhesin genes, Actinobacteria_70_13 may live physically attached to biological material, although whether the material is living or dead is unclear. It is also likely able to decompose fungal biomass based on its repertoire of CAZy genes that target fungal biopolymers (49,51). Biomass decomposition is known to be a common process in the soil, but we are unable to draw further conclusions from genomic data.
The hyphosphere may be an important niche of fungal-bacterial cooccurrence in soil. In heterogeneous soil habitats, the spatial scale over which organisms interact is an important factor that determines how communities assemble (Fig. 4a). This is reflected in the disproportionately higher number of fungal module hubs and the high frequency of occurrence of the same fungal nodes in multiple networks. Since fungal hyphae are much larger than individual bacterial cells, a single hypha or large AMF spores can occupy and influence a larger space and provide surface contact for a diverse group of bacterial cells, forming a hyphosphere in which diverse microbial processes can happen (84,85). Physically attached hyphae and bacterial cells have been captured previously (74,86,87) and in a parallel study (Fig. 4b) (Kakouridis et al., unpublished). Our network analyses may have captured some of these hyphosphere "hot spots" where bacteria physically gather around fungal hyphae, with the fungi acting as module hubs, linking to a range of bacteria. In fact, different modular compartments may reflect hot spots of biotic interactions (88,89), such as a leaky hyphal tip surrounded by bacteria or a particle of decomposing biomass. Along with the rhizosphere and detritusphere, the hyphosphere could represent a unique soil niche, especially in the context of fungal-bacterial interaction. The cooccurrence patterns between fungi and bacteria observed in this study imply that the hyphosphere could harbor processes involving cross-kingdom biotic interactions, whose importance to soil functioning is a currently murky area that demands future study.
In summary, by constructing cooccurrence networks of soil bacteria and fungi associated with the growth of an annual grass, we found that plant growth greatly impacted the development of network structures over time, and we observed different cooccurrence patterns of AMF with bacteria versus nonmycorrhizal fungi with bacteria. Different ecological niches in the complex soil environment are likely major drivers of (1) Soil close to the root, representing the rhizosphere niche, is C rich because of root rhizodeposits but could be nutrient depleted due to plant uptake of N, P, and other essential elements. A wide range of nonmycorrhizal fungi that directly assimilate plant-derived C share this niche with many rhizosphereadapted bacteria. AMF get C by colonizing the root itself and extend beyond the nutrient depletion zone to search for nutrients, showing a low relative abundance and a lack of niche sharing with rhizosphere bacteria. (2) Litter could serve as an AMF inoculum, and the detritusphere could be where AMF-exuded C primes bacterial decomposition of litter, possibly leading to cooccurrence between AMF and bacteria. (3) The higher relative abundance of AMF in bulk soil extends the space around the root where fresh photosynthetic C can be delivered and paves the way for the creation of a hyphosphere niche for bacteria living off of hyphal exudates. Nonmycorrhizal fungi also cooccur and interact with bacteria but more frequently in the rhizosphere, where substrates are abundant and microbial activity is higher. these cooccurrence patterns. In the rhizosphere, the complex and evolving networks between nonmycorrhizal fungi, mainly saprotrophic, and bacteria are probably the result of niche sharing driven by C availability. The rarity of AMF-bacterial cooccurrence in rhizosphere soil resulted from both a low relative abundance of AMF as well as a lack of common life history and potential niche requirements between AMF and rhizosphere-adapted bacteria. Instead, AMF and bacteria were more frequently linked in bulk soil, especially when it contained root litter, leading us to hypothesize that ecological functions involving the participation of both AMF and bacteria could be important in the detritusphere away from roots. SIP-enabled metagenomic analyses suggest that some links potentially involve bacterial symbiosis with eukaryotes, potentially plants or fungi, and their ability to break down fungal biomass, while others are more likely explained by the sharing of ecological niches.
It is possible that actual interactions caused no detectable covariation in the abundances of the partners at the spatial and temporal scales of our sampling regime, posing a strong argument to study biotic interactions using techniques beyond macro/mesoscale correlations in abundance. However, the mesoscale sampling (rhizosphere, bulk, and hyphosphere) applied in this study may be consistent with the scale of fungal niches.
Although caution needs to be used when interpreting molecular marker-based, reconstructed correlation networks (39,40), the networks from our study highlight the unique characteristics of different soil niches in driving the cooccurrence of soil microorganisms and can serve as a useful first step in the exploration of potential associations within soil microbial communities.

MATERIALS AND METHODS
Greenhouse experiment setup and sampling. The design and setup of the greenhouse experiment have been described previously (48) and are diagramed in Fig. S1a in the supplemental material. In brief, microcosms (11.5 by 2.9 by 25.5 cm) were packed with sieved (,2 mm) surface (0 to 10 cm) soils collected from the Little Buck watershed at the University of California's Hopland Research and Extension Center, where the dominant type of vegetation was Avena species, a common type of Californian annual grass. Avena fatua was grown in these microcosms with controlled temperature, water, day length, and atmospheric CO 2 concentration for two consecutive growing seasons in the Oxford Tract Greenhouse at the University of California, Berkeley. In each growth season, both rhizosphere and bulk soils were sampled at five time points, tracking the growth stages of Avena fatua at preplant (week 0), seedling (week 3), vegetative growth (week 6), flowering (week 9), and senescence (week 12). At the end of the first growing season, senescent shoots were cut and removed, and the dead roots were left with the soil inside the microcosms for a 3-month dry period, mimicking California's summer drought, before the second season started and water was supplied again. Thus, dead roots from the first season remained in the microcosms at the beginning of the second season as litter. The rhizosphere soil was sampled as the soil and aggregates that were still attached to roots after gentle shaking. The bulk soil in the first season was collected from inside a 1-mm-mesh bag that was placed in every microcosm when packing, which contained sieved soils and was designed to prevent root penetration. In the second growing season, no "pure" bulk soil existed. After roots with rhizosphere soil were removed, the remaining soil in the microcosm was referred to as "bulk" soil but with the presence of litter. Thus, the rhizosphere and bulk soils in the two seasons represent four types of soil environment: rhizosphere, bulk soil, rhizosphere plus detritusphere, and detritusphere. About 1 to 2 g of soil was obtained per rhizosphere sample, and 5 g of soil was obtained per bulk soil sample. A total of 16 biological replicates were available for each of the combinations of growth season, growth stage, and soil location, resulting in 288 soil samples for analysis.
Bacterial and fungal community analyses. (i) DNA extraction, library preparation, and sequencing. Soil total DNA was extracted in a cetyltrimethylammonium bromide (CTAB) buffer using a phenolchloroform purification protocol and checked for quality using PicoGreen. Bacterial and archaeal communities were surveyed using amplicon sequencing of the V4 region of 16S rRNA genes, amplified using the primer set 515F (59-GTGCCAGCMGCCGCGGTAA-39) and 806R (59-GGACTACHVGGGTWTCTAAT-39). Fungal communities were investigated by sequencing the amplicons of the fungal ITS2 region using the primer set gITS7F (59-GTGARTCATCGARTCTTTG-39) and ITS4R (59-TCCTCCGCTTATTGATATGC-39). Amplification, library preparation, and sequencing were performed on an Illumina MiSeq 2.0 platform in 250-by-2 paired-end format at the Institute for Environmental Genomics, University of Oklahoma.
(ii) Sequence processing. The data processing of 16S raw sequences was reported in detail previously (48). We inherited the use of operational taxonomic units (OTUs) from this published paper instead of using the higher-resolution amplicon sequencing variants (ASVs) so that the 16S nodes are comparable to the prokaryotic networks that have been discussed previously (48). The ITS raw sequences were processed through the pipeline on the Galaxy platform developed by the Institute for Environmental Genomics, University of Oklahoma. First, primers were trimmed with zero mismatch tolerance, followed by Btrim (90) to remove sequences with quality scores of ,25 (window size of 5) or lengths of ,60 bp.
Fungal-Bacterial Cooccurrence across Soil Niches ® Next, forward and reverse sequences were combined using Flash (91), with a minimum required overlap of 10 bp and a maximum expected overlap length of 250 bp. The maximum allowed ratio of mismatch was set at 0.05 in the overlapping region. Afterward, the combined sequences were further cleaned by removing sequences containing ambiguous base N, deleting sequences with lengths of ,100 bp or .500 bp, and deleting chimeras detected using UCHIME (92) with the UNITE (93) representative/reference sequences as the reference database. The remaining sequences were clustered into OTUs by UCLUST (94) at a similarity threshold of 0.97. The longest sequence in each OTU was chosen as the representative sequence of the OTU.
(iii) Plant sequence removal. To identify and remove plant sequences, all the representative sequences were checked using BLASTN (95), nonfungal sequences (unclassified kingdom, plants, and animals) were discarded, and the remaining ITS representative sequences were all fungal sequences. After the removal of nonfungal sequences, the samples from bulk soil had much higher sequence numbers than those from the rhizosphere due to their close association with plant tissue, a known problem with ITS4 primers. Thus, the OTU table was rarified to 4,484 sequences/sample for all the samples to bring even the depth of fungal sequences. A total of 4,936 OTUs were obtained (Fig. S1b).
(iv) Fungal taxonomic assignment. To assign taxonomic information to each OTU, we used the Training Feature Classifier in QIIME2 (96), trained a classifier with the UNITE database (93) v8 (2 February 2019), and classified the fungal ITS representative sequences afterward. One caveat of this classification method is that the confidence score of classification is given to each sequence only at the species level, while higher taxonomic information could be useful in identifying the ecological guild of the fungal taxa (e.g., AMF or nonmycorrhizal). We refined the output taxonomy from the above-described method using the following criteria. First, we accepted species-level assignment with a confidence score of $0.7. Second, we performed a BLAST comparison of our representative sequences to the UNITE v8 database (used to train the classifier) to acquire all the hits with $70% identity to each representative sequence. Third, for representative sequences with a species-level confidence score of $0.3, if more than 80% of the BLAST hits had the same phylum assignment with the classifier result, their phylum assignments were accepted. Fourth, for representative sequences with a species-level confidence score of $0.5, if more than 80% of the BLAST hits had the same class assignment with the classifier result, their phylum assignments were accepted. Taxonomic assignments that do not meet the above-described criteria are reported as "unidentified." The refined results were further checked by (i) randomly picking representative sequences and performing a BLAST search against the NCBI database, which showed a high confidence of taxonomic assignment, and (ii) manually performing a BLAST search of all the Glomeromycota sequences against the NCBI database, which confirmed the assignments for all AMF sequences. In this way, 79% of the 4,936 OTUs were assigned to a phylum, 47% were assigned to a class, 28% were assigned to a genus, and 22% were assigned to a species.
(v) Fungal ecological guild. To assess the ecological guilds of the fungal nodes in the networks, we applied FUNGuild (76) to the taxonomy assignments for the 563 ITS OTUs. A total of 30.9% (174) of these OTUs were assigned guilds at the confidence levels of "highly probable" and "probable." Among them, 52.3% (91 OTUs) belonged to saprotrophs, 20.1% (35 OTUs) were symbiotrophs (including 22 AMF and 13 endophytes), 13.8% (24 OTUs) were pathogens, and the other 13.8% (24 OTUs) had combined guilds. A total of 65.5% (114) of these OTUs have saprotrophic capability. In addition, 62 OTUs were assigned guilds at the confidence level of "possible," and 97% (60) of them have saprotrophic capability.
Cooccurrence bipartite fungal-bacterial network construction using RMT-based correlation cutoffs. (i) Network correlation cutoff. We used the 16 biological replicates to construct a fungal-bacterial cooccurrence bipartite network for each sampling time point for either bulk or rhizosphere soil. To do this, we first combined the 16S and ITS OTU tables to represent one community and discarded the OTUs with ,10 occurrences in 16 replicates to ensure the accuracy of the data association calculation. A total of 1,645 to 2,467 OTUs were included in the construction of each network. Among them, 84 to 90% were 16S OTUs, and 10 to 16% were ITS OTUs. Next, pairwise correlations of all the OTUs were calculated based on the Pearson correlation of their relative abundances in the 16 replicates. To determine the correlation cutoff for the network link, we applied a random matrix theory (RMT)-based approach (45,50). The RMT method detects the transition of the correlation matrix from a noisy system to an organized nonrandom system as the correlation cutoff increases. The correlation cutoff detected ranged from 0.756 to 0.793 for the 18 networks. Pairs of nodes with a correlation above the cutoff are considered linked in the network. Positive correlations indicate positive links, and vice versa. To ensure that the constructed networks are not sensitive to correlation methods, we also used Spearman correlation to construct parallel networks, whose topology showed a trend similar to those from Pearson correlations (Table S1). In this study, we present only networks constructed based on Pearson correlation. The RMTbased correlation cutoff detection procedure is available through the Molecular Ecological Network Analysis Pipeline (MENAP) (http://ieg4.rccc.ou.edu/MENA/).
(ii) Bipartite subnetworks. The 18 networks contained links of bacteria/archaea-bacteria/archaea, fungi-fungi, and fungi-bacteria/archaea. They were scale free, with the degree distribution following power law distributions with R 2 values from 0.83 to 0.94. These networks had 1,096 to 2,020 nodes, with 1,688 to 3,627 links, and average degrees from 2.27 to 5.93. To specifically target fungal-bacterial/archaeal cross-kingdom cooccurrence, subnetworks were extracted to build bipartite networks that contain subsets of nodes and links specifically connecting fungi and bacteria/archaea. We observed only 14 archaeal nodes with 22 links to fungi across the 18 networks. Since most links are between fungi and bacteria, we refer to these bipartite networks as "fungal-bacterial" networks here. A total of 16% to 28% of the links in the full networks were between a bacterial/archaeal node and a fungal node and consist of the fungal-bacterial networks. A total of 29% to 44% of the nodes in the full networks remained in the bacterial-fungal networks.
(iii) AMF-bacterial versus nonmycorrhizal fungal-bacterial networks. Mycorrhizal and nonmycorrhizal fungi represent different life strategies and ecological guilds, which are particularly important to differentiate in environments such as those near plant roots. Thus, we separated the fungal-bacterial networks into AMF-bacterial and nonmycorrhizal fungal-bacterial networks. AMF-bacterial networks contain all the AMF nodes and bacterial nodes linked to them; nonmycorrhizal fungal-bacterial networks contain all the non-AMF fungal nodes with bacterial nodes linked to them.
(iv) Spearman correlation for network construction. In order to explore whether different correlation methods affect network construction, we also used Spearman correlation to construct networks. We show the topological properties of Spearman correlation-based fungal-bacterial networks in Table S1 to demonstrate that the successional trend of networks derived from Pearson and Spearman correlations over time and space are consistent. All the detailed analyses of nonmycorrhizal fungal-bacterial and AMF-bacterial networks are based on Pearson correlation.
(v) Network topological index calculation. The following topological properties of the fungal-bacterial networks were analyzed using the R packages igraph (97) and bipartite (98).

1.
Number of nodes and links. Each node represents an ITS or 16S OTU, whose relative abundance covaries with at least one other OTU at a correlation coefficient above the RMT-detected cutoff threshold. Each link represents one pairwise correlation between two nodes. The total node number (n) and link number (L) were calculated using the functions gorder and ecount. 2.

3.
Average degree (avgK), the average number of links per node, equal to 2L/n. This was calculated using the degree function.

4.
Average path length or geodesic distance (GD), the average length of the shortest path between nodes. This was calculated using the mean_distance function.

7.
Connectance, or graph density, the realized proportion of all possible links. For bipartite networks in this study, all possible links include only those between a bacterial node and a fungal node.
Connectance and the following indices were all calculated using the network-level function.
(vi) Keystone node identification. Based on network topology, nodes that are important in maintaining the network structure are identified as keystone nodes, including module hubs and connectors, as their removal could cause dramatic changes in the cooccurrence patterns of the involved organisms. Operationally, for modular networks (i.e., networks separable into closely interconnected subcomponents called modules), nodes with high within-module connectivity (z i $ 2.5 [highly connected to other nodes within the same module]) are identified as module hubs, while those with high among-module connectivity (p i $ 0.62 [tend to connect nodes in different modules]) are identified as connectors (89).
(vii) Random network generation. To confirm that the observed network topology represents a nonrandom assembly of microorganisms, 100 random networks were generated by rewiring the links among the nodes and compared with the empirical ones for the 18 bipartite fungal-bacterial networks. A random network is generated according to the steps described below. First, for an empirical network with m bacterial nodes, n fungal nodes, and l total links, max(m,n) links were randomly placed between bacterial and fungal nodes so that no isolated node (node with zero degrees) exists. Second, the remaining links [l 2 max(m,n)] are randomly placed between bacterial and fungal nodes. In this way, m, n, and l are preserved, and the random networks are all bipartite. The average degree and connectance are thus preserved between the random and empirical networks. We calculated the means and standard deviations of the R 2 power law and modularity for the 100 random networks. The R 2 power law of the empirical networks was much higher than that of the random networks, and the empirical networks had much higher modularity than the random ones. The distinct topologies of empirical and random networks imply that Fungal-Bacterial Cooccurrence across Soil Niches ® the empirical networks that we constructed are not likely due to random wiring and can reflect microbial community assembly.
Metagenomic sequencing, bacterial genome resolving, and genome matching to network nodes. For the same greenhouse experiment (48), the soil DNA was sequenced for metagenomes and analyzed in a separate study to construct metagenome-assembled genomes (MAGs) (49). Briefly, a total of 55 partial bacterial MAGs with .70% completeness and ,10% contamination were reconstructed from these sequences, and high-quality 16S rRNA sequences for 13 of these MAGs were obtained (49). In order to match these MAGs to the network nodes, we performed a BLAST search of our representative sequences (16S V4 region) of the bacterial nodes in the network against the 14 high-quality 16S sequences obtained from the MAGs. Five MAGs were found matching with 9 bacterial nodes with .99% similarity of the 16S V4 regions for a length of at least 200 bp.
Statistical analysis. Statistical analyses were performed using R version 3.5.3 (99). To compare the network topological indices between seasons, over plant growth, and between rhizosphere and bulk soils, an analysis of covariance (ANCOVA) model was applied to each network topological parameter using the aov function in the stats package. In the model "topological parameter ; Season*Week*Location," season and location are considered categorical variables, and week is considered a continuous variable. To further determine the pattern that network topological indices change over time with plant growth, a general linear model, "topological parameter ; Week," was fit for each index using the lm function in the stats package.
Data availability. The 16S and ITS amplicon sequencing data are available through the NCBI SRA database under BioProject accession number PRJNA703382. R code for network construction and data analysis is available on GitHub at https://github.com/Mengting-Maggie-Yuan/Fungal-bacterial-network -construction-and-analysis.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.