Biological soil crusts on agricultural soils of mesic regions promote microbial cross-kingdom co-occurrences and nutrient retention

Introduction Biological soil crusts (biocrusts) are known as biological hotspots on undisturbed, nutrient-poor bare soil surfaces and until now, are mostly observed in (semi-) arid regions but are currently poorly understood in agricultural systems. This is a crucial knowledge gap because managed sites of mesic regions can quickly cover large areas. Thus, we addressed the questions (i) if biocrusts from agricultural sites of mesic regions also increase nutrients and microbial biomass as their (semi-) arid counterparts, and (ii) how microbial community assemblage in those biocrusts is influenced by disturbances like different fertilization and tillage regimes. Methods We compared phototrophic biomass, nutrient concentrations as well as the abundance, diversity and co-occurrence of Archaea, Bacteria, and Fungi in biocrusts and bare soils at a site with low agricultural soil quality. Results and Discussion Biocrusts built up significant quantities of phototrophic and microbial biomass and stored more nutrients compared to bare soils independent of the fertilizer applied and the tillage management. Surprisingly, particularly low abundant Actinobacteria were highly connected in the networks of biocrusts. In contrast, Cyanobacteria were rarely connected, which indicates reduced importance within the microbial community of the biocrusts. However, in bare soil networks, Cyanobacteria were the most connected bacterial group and, hence, might play a role in early biocrust formation due to their ability to, e.g., fix nitrogen and thus induce hotspot-like properties. The microbial community composition differed and network complexity was reduced by conventional tillage. Mineral and organic fertilizers led to networks that are more complex with a higher percentage of positive correlations favoring microbe-microbe interactions. Our study demonstrates that biocrusts represent a microbial hotspot on soil surfaces under agricultural use, which may have important implications for sustainable management of such soils in the future.

Introduction: Biological soil crusts (biocrusts) are known as biological hotspots on undisturbed, nutrient-poor bare soil surfaces and until now, are mostly observed in (semi-) arid regions but are currently poorly understood in agricultural systems. This is a crucial knowledge gap because managed sites of mesic regions can quickly cover large areas. Thus, we addressed the questions (i) if biocrusts from agricultural sites of mesic regions also increase nutrients and microbial biomass as their (semi-) arid counterparts, and (ii) how microbial community assemblage in those biocrusts is influenced by disturbances like different fertilization and tillage regimes.
Methods: We compared phototrophic biomass, nutrient concentrations as well as the abundance, diversity and co-occurrence of Archaea, Bacteria, and Fungi in biocrusts and bare soils at a site with low agricultural soil quality.
Results and Discussion: Biocrusts built up significant quantities of phototrophic and microbial biomass and stored more nutrients compared to bare soils independent of the fertilizer applied and the tillage management. Surprisingly, particularly low abundant Actinobacteria were highly connected in the networks of biocrusts. In contrast, Cyanobacteria were rarely connected, which indicates reduced importance within the microbial community of the biocrusts. However, in bare soil networks, Cyanobacteria were the most connected bacterial group and, hence, might play a role in early biocrust formation due to their ability to, e.g., fix nitrogen and thus induce hotspot-like properties. The microbial community composition differed and network complexity was reduced by conventional tillage. Mineral and organic fertilizers led to networks that are more complex with a higher percentage of positive correlations favoring microbe-microbe interactions. Our study demonstrates that biocrusts represent a microbial hotspot on soil surfaces under agricultural use, which may have important implications for sustainable management of such soils in the future. KEYWORDS biocrust, amplicon, hotspot, managed site, Cyanobacteria, microbial co-occurrence network, mesophilic region

Introduction
Biological soil crusts (biocrusts) are an assemblage of soil particles, photoautotrophic primary producers, heterotrophic microorganisms and microfauna on or within the first few millimeters of the soil surface (Weber et al., 2022). They are particularly abundant in arid and nutrient-poor environments (Belnap, 2003;Maier et al., 2018). In those ecosystems biocrusts fulfill important functions related to carbon and nitrogen fixation, the storage of water and nutrients as well as the induction of soil formation and stabilization (Belnap et al., 2001;Costa et al., 2018). Phototrophic biota such as Cyanobacteria, microalgae, or mosses are essential key-players to provide such functions (Miralles et al., 2012). For example, Cyanobacteria can stabilize the soil matrix due to their filamentous growth (Chamizo et al., 2018;Jung et al., 2018) and the production of sticky exopolysaccharides (Cania et al., 2019a). The fixation of carbon and nitrogen attracts heterotrophic microbes including Bacteria, Fungi, Archaea and Protists Maier et al., 2018;Roshan et al., 2021). As a consequence, Cyanobacteria are often described as keystone taxa of biocrusts in studies, which investigated microbial community composition and performed correlation network analyses (Chilton et al., 2018;Pombubpa et al., 2020).
There is increasing evidence that biocrust formation is not confined to nutrient-poor and arid regions. Recent studies identified biocrusts in managed ecosystems of mesic regions (Gall et al., 2022) like forests (Baumann et al., 2017;Glaser et al., 2018;Kurth et al., 2021;Glaser et al., 2022a). In contrast to arid biocrusts, they occur here as ephemeral stages. In forests, it was demonstrated that the biomass of biocrusts quickly increased in spring before herbal plant growth or after disturbance events like tree cutting or wind driven tree fall, which gives phototrophic biota of biocrusts a selective advantage (Kurth et al., 2021). Similar to forests, agroecosystems also provide potential niches for biocrust development, such as the time between harvest and sowing or between the rows of broad leave crops like corn, potatoes, or sugar beets. Besides the potential time and space for biocrust development in agroecosystems, multiple studies have reported that single components of agricultural management including tillage (Curaqueo et al., 2010;Wang et al., 2016;Wagg et al., 2018) and fertilization are detrimental for microbial community assembly, which might hamper the development of biocrust communities Schulz et al., 2013;Maier et al., 2018). Additionally, the previously described central role of Cyanobacteria in microbial networks might be obsolete in biocrusts of fertilized soils, because of the high external nutrient input. However, 'on-farm' studies, which investigate the combined influence of tillage intensity and fertilizer type or amount on microbial community composition, are still missing. Taking the positive effect of biocrusts on many ecosystem functions and their potential importance for sustainable agricultural management into account, there is a strong need to overcome this limitation.
Thus, this study aimed to investigate the combined effect of mineral fertilization and organic treatments as well as tillage intensity on: (1) the composition and correlation of microbial communities in biocrusts of an agricultural field and (2) the ability of those biocrusts to further increase the amount of available nutrients. Samples were taken under the auspices of a long-term fertilization experiment in Germany, which combined different tillage intensities (minimal, reduced, conventional tillage) with different fertilization treatments (different levels of mineral nitrogen fertilizer, with or without crop residue retention). Long-term data revealed already that reduced tillage and crop residue retention increased carbon stocks and sugar beet yields in this experiment (Armbruster et al., 2009(Armbruster et al., , 2012. Sampling was conducted in autumn at the end of sugar beet cultivation but before harvest, which allowed a maximum period for biocrust development. Samples were taken between sugar beet rows. We analyzed nutrient concentrations, the community composition of Archaea, Bacteria and Fungi as well as their co-occurrence and mutual exclusion pattern in biocrusts in comparison to the surrounding bare soil.

Experimental field and sampling procedure of biocrusts
The current study was conducted 2016 at the agricultural experimental farm "Rinkenbergerhof " in Germany (49°21′34.0"N 8°25′14.8″E), which belongs to the Agricultural Investigation and Research Institute Speyer ( Figure 1). Samples were taken in frame of the long-term "International Organic Nitrogen Fertilization Experiment (IOSDV), " which is carried out since 1983. The site is characterized by a mean annual temperature of 10.0°C and a mean annual precipitation of 593 mm. The soil of the IOSDV experiment has been described as Cambisol consisting of equal amounts of silt and sand and 9% of clay. It has a field capacity of 10%. The German arable assessment ("Ackerzahl") evaluated the soil with 25 to 30, which represents sites with a low soil quality (Bischoff and Emmerling, 1997;Armbruster et al., 2012;Schmid et al., 2018). Phosphorus, potassium and magnesium were applied in equal amounts to all plots (31 kg of phosphorus ha −1 yr. −1 , 121 kg of potassium ha −1 yr. −1 and 28 kg of magnesium ha −1 yr. −1 ). In 2016 during sugar beet cultivation pesticides including fungicides were used based on the particular needs and irrigation of 110 mm was applied to avoid drought damage of the sugar beet plants.
The IOSDV experiment investigates the combination of mineral fertilization and organic treatments (since 1983) as well as tillage intensity (since 2004) in a 3 years crop rotation of sugar beet, winter wheat, and winter barley (Bischoff and Emmerling, 1997;Armbruster et al., 2012;Schmid et al., 2018). The experiment is designed in a fullfactorial block design with three replicates per treatment and a plot size per treatment of 6 m * 7.5 m. In frame of the experiment, we sampled the following treatments: (1) 120 and 240 kg N ha −1 yr. −1 mineral nitrogen fertilization (calcium ammonium nitrate). The applied rates are typical rates used in low-and high input agriculture in the region of the study, respectively. (2) We compared the additional organic treatment (+org) to control plots (−org). This was characterized by the retention of crop residues after harvest and the cultivation of a cover crop (Raphanus sativus var. oleiformisis) every third year after winter barley cultivation. (3) Tillage intensity was investigated by sampling plots with reduced (rT) and conventional tillage (cT). During rT the soil is broken up by a cultivator to a depth of 10 cm without turning the soil and under cT, the soil is plowed to 30 cm. Because of the different treatments, carbon stocks varied between the treatments. Sampling took place on October 5 in 2016. The sampling time was chosen to allow for a maximal development of biocrusts and prior to the harvest of sugar beets, which would destroy the biocrusts, because of the heavy machinery used. On some plots, the biocrusts were so prevalent, that reference sampling of bare soil was challenging, which might have been caused by high precipitation of 25 mm in the 2 weeks before sampling. However, during the sampling and right before there was no rain. In total, 24 biocrusts and 24 bare soil samples were taken, which consisted of three replicates of the eight treatment combinations: cT 120 −org, cT 120 +org, cT 240 −org, cT 240 +org, rT 120 −org, rT 120 +org, rT 240 −org, rT 240 +org (see Figure 1). Biocrusts were visually identified as green covered soils according to the biocrust definition of Weber et al. (2022). Biocrusts were sampled from the top millimeters of the mineral soil between the sugar beet rows and were lifted with a spatula and sampled as a coherent layer of approximately 3 mm thickness. Most of the biocrusts were located in tractor traffic lanes and in little grooves (Supplementary Figure S1). Green biofilms on macroscopic organic matter like decomposed leaves were not considered as biocrusts. Samples from biocrust-free areas (no visible phototrophic biomass), referred to as bare soil, were taken from the top 3 millimeters. A composite of five biocrust and bare soil samples per plot was homogenized to reach sufficient amounts of sampling material. All samples were sieved to a particle size of 2 mm. Approximately, 12 g of fresh soil was stored at 4°C for chemical analysis and ca. 2 g was shock frozen on dry ice in the field and stored at −80°C for subsequent molecular analyses.

Chemical and physical soil parameters
Soil pH was analyzed as described in DIN ISO 10390 (International Organization for Standardization, 2005) with the electrode SenTix 61 and pH meter (inoLab pH 720 Level 1, Wissenschaftlich-Technische Werkstätten, Weilheim, DE) in 0.01 M CaCl2 extract of 2 g of fresh soil (Stempfhuber et al., 2015). Nitrate, ammonium, dissolved organic nitrogen (DON) and carbon (DOC), were extracted with 0.01 M CaCl2 solution (in a 1:4 ratio) from 3 g of fresh soil by 45 min overhead shaking at 67 rpm (Reax 2, Heidolph Instruments, Schwabach, Germany), followed by filtration through Whatman™ 595 1/2 filters (Sigma-Aldrich, St. Louis, MO, United States) and determined by a segmented flow analyzer (Skalar SANPlus 5100 with autosampler 1050, Skalar analytic, DE, EU). DOC and DON were analyzed on DIMATOC2000 (DIMATEC Analysentechnik, Essen, DE). To calculate dissolved inorganic carbon (DIC), one drop of 32% HCl was added to the filtrate prior to analysis and the difference to DOC without HCl was used to calculate DIC (Brankatschk et al., 2011). Chlorophyll a was used as a proxy for phototrophic biomass and determined according to Ritchie (2008). The procedure was as follows, 0.7 g frozen soil was extracted in 3 mL 96% aqueous ethanol (v/v) and incubated for 30 min at 78°C. Afterwards, the extract was centrifuged at 5°C at 6,000 rpm. The absorption of the supernatant was measured at wavelengths of 632 nm, 649 nm, 665 nm, 696 nm using the spectrophotometer UV-2401PC (Shimadzu, Kyoto, JPN). Extraction was repeated until no chlorophyll could be detected in the supernatant anymore and the content was summed up for all extraction steps.

DNA extraction
DNA was extracted from 0.3 g of soil using the NucleoSpin ® Soil kit (Macherey-Nagel, Düren, DE) according to the manufacturer's manual using lysis buffer SL2 and 150 μL of enhancer. DNA quality was assessed using a UV-VIS spectrophotometer (NanoDrop ® ND-1000 spectrophotometer, Thermo Fisher Scientific, Waltham, Massachusetts, United States). DNA concentration was determined using the Quant-IT™ Picogreen ® dsDNA Assay Kit (Thermo Fischer Scientific). A negative extraction control sample without soil was processed as negative control to check for possible contaminations of chemicals from the kit used for nucleic acid extraction.

Quantification of prokaryotic and fungal biomass
Real-time qPCR was used to determine the abundance of Bacteria, Archaea and Fungi, which was used as a proxy for their biomass. For Bacteria (Bach et al., 2002) and Archaea (Bano et al., 2003;FIGURE 1 Picture of a biocrusts from a plot with conventional tillage, without crop residues and a mineral fertilization amount of 240 kg N ha −1 yr. −1 , taken on October 5th, 2016. Integrated is a map of Germany showing the location of the experimental field. The table below the figure shows the details of treatment combinations for tillage, mineral fertilizer and organic treatment. Tillage management is compared between conventional (cT) and reduced tillage (rT). The effect of mineral fertilizer amounts is evaluated in levels of 120 and 240 kg N ha −1 yr. −1 . Organic treatment as crop residues (+org) are compared to a control plot (−org) where crop residues were removed after harvesting.
SYBR Green ® based assays (Applied Biosystems, Foster City, CA, United States) were performed on a 7,300 Real-Time PCR System (Applied Biosystems). Details about forward (F) and reverse (R) primer, reaction conditions, and calibration standards are summarized in Supplementary Table S1. Primers were purchased from Metabion (Planegg, Germany) and bovine serum albumin (BSA) from Sigma-Aldrich (Missouri, United States). To exclude inhibitory effects, dilution tests were performed. Standard series (10 6 to 10 2 gene copies μl −1 ) and samples diluted to 1/32 were included in each run. To check for possible contamination of chemicals used, a negative control without DNA of samples was included, as well. To evaluate the quality of the qPCR, melting curve analyses were performed and randomly chosen samples were checked by electrophoresis on a 1.5% agarose gel. The amplification efficiency was calculated by ε = 10 (−1/slope) −1 with the slope of the standard series, and was 78-84% for all genes. The r 2 of the standard curve was >0.987.

Diversity of prokaryotes and Fungi
As no specific primers exist, which separately target Bacteria and Archaea for Illumina MiSeq ® sequencing (Illumina, San Diego, CA, United States), we used the universal primers Arch0519 F (Klindworth et al., 2013) and Pro 805 R (Herlemann et al., 2011), which were optimized for the simultaneous sequencing of both prokaryotes, Bacteria and Archaea, and target the V4 region of the 16S rRNA gene. For Fungi, the ITS 3 primer mix and ITS 4 primer mix (Tedersoo et al., 2015) were used. Primer sequences are given in Supplementary Table S2 were pooled for sequencing on the MiSeq ® instrument with 2 × 300 bp (Illumina) using MiSeq Reagent kit v3 (600 cycles) and spiked with 20% PhiX as a control. Samples with less than 10,000 reads were re-sequenced. Sequencing adapters were removed using AdapterRemoval (Schubert et al., 2016). DADA2 package Version 1.8.0 (Callahan et al., 2016) in R Version 3.6.1 (R Core Team, 2022) was used for length and quality filtering including PhiX and chimera removal. For prokaroytes, forward and reverse reads were trimmed at 10 and 250 bp and 10 and 200 bp, respectively. For Fungi, the parameters were set to 10 and 275 bp and 10 and 225 bp, respectively. Alignment of mate pairs, inferring into amplicon sequencing variants (ASVs) and merging of sequencing runs was done prior to taxonomy assignment. For prokaroytes, this was done against the Silva database Version 132  and for Fungi, against the Unite database Version 7.1 [ -11-20, Kõljalg et al. (2013]. After taxonomy assignment, reads present in samples of blank extraction or PCR negative controls were removed (11 of 9,637 ASVs for prokaroytes, 6 of 3,715 ASVs for Fungi) as well as ASVs assigned as Chloroplasts or Mitochondria (on average 11.99%). To differentiate between Chloroplast sequences (representative chloroplast reads from plants, Algae and Bacteria) and cyanobacterial ASVs detected in our study, a phylogenetic tree was calculated and only sequences clearly assigned as cyanobacteria were processed (Supplementary Figure S2). Random subsampling was performed to the lowest number of reads (prokaroytes: 19,722 reads, Fungi: 22,731 reads) using the function rarefy_even_depth() of the R package phyloseq Version 1.30.0 (McMurdie and Holmes, 2013) with the random seed set to 3,006. As they were no longer present after subsampling, 448 ASVs out of 9,626 ASVs for prokaroytes and 283 ASVs out of 3,709 ASVs for Fungi have been removed.
The sequence data was submitted to NCBI via the Sequence Read Archive (SRA) and is available under the accession number PRJNA646655.

Statistical data analysis
Data analysis was performed with R Version 3.6.1 (R Core Team, 2022). The vegan package Version 2.5-7 was used to calculate rarefaction curves, alpha diversity [Shannon diversity (S) (Shannon, 1948), richness as number of ASVs (R) and Pielou's eveness (P) (Pielou, 1966)] of the subsampled data. To test for the normal distribution of the data residual vs. fitted plots and sample quantile vs. theoretical quantile plots were tested for normal distribution and homogenous variance to verify the models. If needed, data was transformed to meet normal distribution with the 1 + log transformation for abiotic parameters, gene abundances, and alpha diversity; square root transformation for taxa on family level with abundances of more than 2% for at least one replicate of Bacteria and 3% of Fungi. Linear models were applied to detect significant differences according to biocrust presence or treatment for soil parameters, absolute abundances of microbial groups, alpha diversity (S, R, P) and relative abundance of microbial families obtained from sequencing results. To disentangle treatment effects (tillage, mineral fertilization amount, and organic treatment), biocrust and bare soil samples were separated.
Frontiers in Microbiology 05 frontiersin.org Difference in β-diversity of the treatments was calculated by Bray-Curtis-distance, plotted as nonmetric multidimensional scaling (NMDS) plots (Oksanen et al., 2017) and tested for significance by PERMANOVA. The barplots and heatmaps were created using ggplot2 package Version 3.3.5 (Wickham, 2009). CoNet (Faust and Raes, 2016) as an add-on of Cytoscape Version 3.7.2 (Shannon et al., 2003) was used to calculate correlation networks on relative abundance data. For network calculations, ASVs unclassified at the phylum, class and order level were removed. Relative abundance was calculated separately for prokaryotes and Fungi and data sets were merged afterwards. Networks were calculated separately for the two sample (biocrust and bare soil) and management types (tillage: cT and rT, mineral fertilizer: 120 and 240, organic treatment: -org and + org) to detect differences in co-occurrence patterns (N = 12). To obtain correlations only valid for all treatments, ASVs present in less than 10 samples were excluded from analysis. Significant correlations (co-occurrences vs. mutual exclusion) were determined based on Pearson and Spearman correlation as well as Bray-Curtis and Kulback-Leibler dissimilarity. At least two of the analyses need to support the link connecting two nodes. 1,000 permutations and bootstrap scores were generated with Brown value of p merging (Brown, 1975). Gephi 0.9.2 (Bastian et al., 2009) was used for visualization of undirected networks with the Fruchterman-Reingold layout (Fruchterman and Reingold, 1991). Correlation partners are visualized with the chordDiagram function of circlize package Version 0.4.13 in R (Gu et al., 2014). Network hubs were defined as highest connected nodes (Tipton et al., 2018), where hubs not specified at the family level and lower were ignored.

Basic properties of biocrusts
The content of phototrophic biomass measured as chlorophyll a was 6-times higher in biocrusts compared to bare soil (F = 183.7, p < 0.001) reaching 17.87 ± 7.2 μg g −1 soil dry weight (dw) in biocrusts compared to 2.98 ± 1.2 μg g −1 dw in bare soil ( Figure 2; Table 1;  Supplementary Table S3). The management had no significant influence on the chlorophyll a content. The microbial biomass, determined as the amount of extracted DNA (F = 18.5, p < 0.001) and the abundance of Bacteria, Archaea and Fungi was significantly higher (18.7 ≥ F ≤ 106.8, p < 0.001) in biocrusts compared to bare soils (Supplementary Table S4). Bacterial 16S rRNA genes (2.3*10 10 ± 8.4*10 9 copies g −1 dw) were significantly (F = 277, p < 0.001) more abundant, than Fungal ITS sequences (2.4*10 9 ± 1.3*10 9 copies g −1 dw) and archaeal 16S rRNA genes (1.1*10 7 ± 8.9*10 6 copies g −1 dw) ( Figure 3; Table 1;  Supplementary Table S5). In biocrusts, the bacterial and fungal abundance was not influenced by management. Contrastingly in bare soils, tillage, mineral fertilizer amount and organic treatments altered archaeal, bacterial and fungal biomass. For example, in bare soils fungal biomass was significantly higher with organic treatments (F = 6.7, p = 0.02). Further, an interacting treatment effect for mineral fertilizer and organic treatments were observed and in samples with 120 kg N ha −1 yr. −1 the fungal abundance was significantly reduced with organic treatment (F = 6.1, p = 0.025). Although, a similar pattern was observed for bacterial abundance this was not significant (F = 3.6, p = 0.059). Archaeal biomass was reduced in biocrusts and bare soils in the treatment with mineral fertilizer of 240 kg N ha −1 yr. −1 (F > 10.8, p ≤ 0.05).

Microbial diversity in biocrusts
In total, 3,601,616 reads were obtained from 16S rRNA gene amplicon sequencing and 4,995,629 reads from ITS amplicon sequencing. The average loss of reads because of the different filter steps during data processing was 18.9% for prokaryotes and 29.4% for Fungi (Supplementary Tables S6, S7). However, after subsampling, all rarefaction curves still reached a plateau indicating a sufficient sampling depth for further analysis (Supplementary Figure S3). Details of all sequencing runs including number of demultiplexed, filtered and merged reads are summarized in the Supplementary Tables S6, S7. For α-diversity indices of prokaryotes, no difference between bare soils and biocrusts was detected (F < 0.6, p > 0.4). In contrast, fungal diversity and richness were significantly (F > 8.8, p < 0.005) lower in biocrusts (S = 3.97 ± 0.35, R = 274 ± 88) compared to bare soils (S = 4.35 ± 0.43, R = 390 ± 110), while evenness was also not affected. Regarding the management, the evenness of prokaryotes was reduced in cT compared to rT in both compartments (F > 6, p < 0.03) (Supplementary Tables S8, S9), while tillage had contrasting effects on fungal α-diversity in bare soils and biocrusts. For example, rT compared to cT caused a decrease in fungal richness in bare soils (F = 11.9, p = 0.003) and an increase in biocrusts (F > 6.6, p < 0.02). Interacting treatment effects (p < 0.05) were detected in biocrusts for archaeal/bacterial richness (Tillage * Mineral Fertililzer Amount * Organic Treatment) and fungal diversity and richness (Tillage * Organic Treatment) (Supplementary Tables S8, S9).
Beta-diversity analysis revealed significant differences between biocrust and bare soil microbial communities (p < 0.02) as well as between the different treatments (p < 0.05) ( Table 2; Supplementary Figures S4, S5). Regarding the prokaryotic community, the presence of biocrusts reduced treatment effects on the community composition. The fungal community composition was affected by all treatments no matter if biocrusts or bare soils were considered (p < 0.05).
Overall, 19 bacterial families were significantly (p < 0.03) different between biocrusts and bare soils ( Figure 4B). Eleven of those were higher in relative abundance in biocrust samples, including many reads which were assigned as α-and γ-Proteobacteria, while reads corresponding to Acidobacteria and Actinobacteria were lower abundant in biocrusts compared to bare soils. Interestingly, the relative abundance of ASVs linked to cyanobacterial families was similar in biocrusts (16.8 ± 9.0%) and bare soils (16.0 ± 13.3%).
The amount of mineral fertilizer addition significantly (4.6 < F < 23.0, p < 0.05) changed the relative abundance of 54% of the highly abundant bacterial families in biocrusts and bare soils. The identity of the affected families was mostly the same in biocrusts and bare soils and included the highly abundant Burkholderiaceae and Chitinophagaceae. Regarding Cyanobacteria, the relative abundance of Coleofasciculaceae and Nostocaceae was seven times higher in Copies of target region [per gram soil dry weight (dw)] 16S rRNA gene for Archaea and Bacteria and ITS region for Fungi in bar plots as the average of replicates (n = 3) and error bars displaying the standard deviation. Treatments are given for tillage [conventional (cT) vs. reduced tillage (rT)], organic treatment [without (−org) vs. with crop residues (+org)] and mineral fertilization amount (120 vs. 240 kg N/ha·a).
Frontiers in Microbiology 08 frontiersin.org treatments with 120 kg N ha −1 y −1 compared to 240 kg N ha −1 y −1 in both biocrusts and bare soils (F > 5.8, p < 0.02). In contrast, the effect of crop residues was less pronounced, especially in biocrusts only 10% (compared to 26% in bare soils) of the dominant bacterial families were positively affected and had higher abundances due to the organic treatment. Of those, Pseudonocardiaceae, Xanthobacteraceae, Chthoniobacteraceae, and Rubritaleaceae were affected by mineral fertilizer and organic treatmentor the interaction of both (Pseudonocardiaceae, Chthoniobacteraceae) ( Figure 4B). Tillage significantly (15.1 < F < 65.3, p < 0.04) changed the relative abundance of 15 bacterial families in biocrusts, but only seven in bare soils. For example, in biocrusts Chitinophagaceae were increased by cT (4.7 ± 1.4% to 9.4 ± 4.3%) (F ≥ 23.0, p < 0.001) while the "Unknown Family" of Oxyphotobacteria (including also species of Leptolyngbya in Silva database v132) decreased under rT (8.8 ± 3.3% to 5.3 ± 3.7%) (F ≥ 8.0, p ≤ 0.01). Interacting treatment effects were mainly observed, where one of the factors alone significantly influenced the abundance of bacterial families (Supplementary Table S10). Only the family env. OPS_17 (Bacteroida of Bacteroidetes) in biocrusts was significantly (F = 6.4, p = 0.027) affected by all treatments in combination but not by one single one. Longimicrobiaceae in biocrusts were significantly (F = 6.0, p = 0.027) affected by tillage and in combination with mineral fertilization amount without these two affecting the family individually. In general, more interacting effects were observed in biocrusts (n = 11) than in bare soils (n = 5).
Two hundred and twenty different fungal families were detected. In our analyses we focused on the 29 most abundant ones (at least 2% abundance in one of the samples) ( Figure 4C). Predominant families were Nectriaceae (bare soils 22.6 ± 6.9%, biocrusts 23.1 ± 7.2%), Pleosporaceae (bare soils 14.4 ± 5.6%, biocrusts 19.4 ± 7.1%), Heat maps of (A) archaeal, (B) bacterial, and (C) fungal community composition in bare soils and biocrusts at the family level are displayed in the left two columns (for Bacteria (B) and Fungi (C) with an abundance above 2% in at least one replicate). The average of field replicates (n = 3) is presented. Values range from white (0%), orange, and dark red to the highest abundance in purple (A) 1%, (B) 16%, (C) 35%). Significant variations (p < 0.05) are shown in the right columns for the biocrust effect (bare soil vs. biocrust) (orange), mineral fertilizer amount (light blue), organic treatment (dark blue), and tillage (green) separately for bare soils and biocrusts. The families additionally influenced by interacting treatment effects are supplied in Supplementary Tables S10, S11.

Cross-kingdom correlation networks in biocrusts
Network analysis represents co-occurrence and mutual exclusion patterns of the microbial community based on significant correlations (edges) between taxa (nodes) detected in >80% of the samples. Networks were calculated separately for both compartments (biocrusts and bare soils) and management effects (cT vs. rT, 120 vs. 240 N kg ha −1 yr. −1 , −org vs. +org), resulting in 12 correlation networks based each on 12 samples given in co-occurrence and mutual exclusion (Figures 5, 6). The total number of nodes (in average for bare soils: 209, biocrusts: 152) and edges (bare soils: 367, biocrusts: 238) was higher in bare soils for all treatments (Table 3). Especially rT resulted in the highest numbers of edges and nodes in bare soils reaching 645 and 255, respectively. The same positive effect of rT on the number of edges and nodes was observed for biocrusts (nodes: 223, edges: 417). In biocrusts, the detrimental effect of cT were most pronounced as the number of edges and nodes dropped 1.5 times more compared to bare soils. Regarding the nature of the correlations (co-occurrence vs. mutual exclusion) an interesting pattern emerged based on the comparison of fertilization treatments and tillage: biocrusts had a higher proportion of co-occurrences compared to bare soils in all treatments related to mineral fertilization or organic treatment peaking in the treatments with highest application rates like 240 (77.5%) and + org (76.8%). In contrast, rT and cT networks revealed a higher proportion of co-occurrences in bare soils with the highest value found under rT (77.3%).
Bacteria were the dominating kingdom in our network analysis forming up to 82% of nodes. In 11 of 12 networks, the ratio of nodes assigned to Bacteria and Fungi ranged between 2 and 3.6, with cT being the only exception. Here, the Bacteria:Fungi ratio dropped to 1.4 in bare soils, while it strongly increased in biocrusts to a ratio of 5. This could be assigned to a reduction of Tremellomycetes and Dothideomycetes ASVs in the network. With the increasing share of Fungi in the bare soil network, bacterial-fungal correlations increased under cT, mostly between different Basidiomycota and Actinobacteria, Oxyphotobacteria and Phycisphaerae.
In total, 31 different families of Archaea, Bacteria and Fungi in biocrusts and 25 in bare soils were identified as network hubs (highest degree), where 12 families appeared in both compartments as hubs (Figures 5, 6, highlighted in bold in Supplementary Table S12). Out of the 31 hub families, 19 and 13 were specific in biocrusts and bare soils, respectively. Most of these unique hubs only appeared in one or two networks. The only exceptions of the specific families were Blastocatellaceae in bare soils, which were identified as hubs in 240, − org, cT, and rT networks and Rhodobacteraceae in biocrusts, which formed hubs in networks 120, −org, +org, and cT.

Biocrusts are richer in nutrients and microbial abundance than bare agricultural soils
Biocrusts have previously been described in arid or oligotrophic environments where the accumulation of nutrients is one of their characteristic properties Schulz et al., 2013;Maier et al., 2018). To investigate whether this also holds true in fertilized, mesic agricultural ecosystems was one major aim of this study. We observed that under sugar beet cultivation in mesic regions, biocrusts could form large patches (Figure 1; Supplementary Figure S1) with high contents of phototrophic biomass, which was comparable to those of arid biocrusts (Román et al., 2019). Biocrust formation was accompanied by increased DOC, nitrate and DON concentrations as well as higher microbial biomass compared to bare soil implying that these typical biocrust properties are also present in these mesic biocrusts. However, we cannot exclude that this effect changes over the course of a season or at different soil types, which retain nutrients better than these poor soils with high amounts of sand (Kurth et al., 2021). The increase in nutrient concentrations might be attributed to two different mechanisms. First, due to the activity of the phototrophic biomass, CO 2 is fixed and increased DOC concentrations. This is supported by the positive correlation of DOC (Supplementary Figure S6) and chlorophyll a concentration in the biocrusts. Additionally, the significantly higher relative abundance of potential N 2 fixing bacteria like Rhizobiaceae might further facilitate nitrogen input as shown for other biocrusts previously (Lan et al., 2011;Maier et al., 2018;Román et al., 2019). Second, the polymeric matrix of biocrusts traps or retains nutrients more efficiently (Costa et al., 2018;Rossi et al., 2018;Cania et al., 2019a). This is corroborated by the higher relative abundance of Bacteria potentially involved in the production of extracellular polysaccharides in the biocrusts like Burkholderiaceae, Flavobacteriaceae, Chitinophagaceae, Leptolyngbyaceae, and Rhizobiaceae (Cania et al., 2019a;Vuko et al., 2020). Although microbial biomass and nutrient concentrations were generally higher in biocrusts, alpha-diversity and network density were lower (Figures 3, 5, 6; Supplementary Tables S6, S7). These findings contradict previous biocrust studies, where microbial diversity increased with ongoing biocrust development (Chilton et al., 2018;Maier et al., 2018). However, these studies were conducted in arid or nutrient scarce environments, where biocrusts served as habitable micro-ecosystems and paved the way for the attraction of further microbes. In line with our finding, Glaser et al. (2022b) also found reduced diversity in nutrient-poor, mesic dunes. This data, as well as ours from agricultural biocrusts of mesic regions, support findings from other nutrient-rich hotspots in managed ecosystems like the rhizosphere, drilosphere (Uksa et al., 2015) and biocrusts from mesic forests (Glaser et al., 2022a), where only a subset of the diverse bulk soil community was selected by the specific hotspot. Comparable to those hotspots, correlation networks were dominated by Sphingomonadaceae (α-Proteobacteria), Burkholderiaceae (γ-Proteobacteria), Chitinophagaceae (Bacteroidia), Nocardioidaceae (Actinobacteria), and Pleosporaceae (Dothideomycetes) (Supplementary Table S12). These families were described as key microbiota during organic matter decomposition (Banerjee et al., 2016), particularly the degradation of chitin and cellulose (Coenye, 2014;Glaeser and Kämpfer, 2014;Maier et al., 2018;Wieczorek et al., 2019). Thus, they may be advantageous for organic C turnover in biocrusts, which were introduced by organic treatment or activities of phototrophic microbes. Furthermore, Sphingomonadaceae, Burkholderiaceae and Chitinophagaceae or Pleosporaceae and Nectriaceae were further detected as dominating families in exo-and lipopolysaccharides production, especially in developing biocrusts, during soil formation and under conventional tillage (Osińska-Jaroszuk et al., 2015;Xiao and Zheng, 2016;Cania et al., 2019a,b;Vuko et al., 2020). Actinobacteria were important nodes in biocrust networks. However, they were not among the high abundant taxa, which highlights that rarely abundant taxa may contribute essential functions to microbial communities (Shi et al., 2016;Benjamino et al., 2018). Indeed, they were correlated to other proteo-, cyanobacterial or fungal classes but most remarkably was the high proportion of their self-co-occurrences. Though, they share the same niches and promote the growth of other groups of their own phylum (Berry and Widder, 2014;de Vries and Wallenstein, 2017). The actinobacterial network participation was much less under conventional tillage compared to reduced tillage. We observed a drop of 30.7% in bare soils and of 11.6% in biocrusts (see Figures 5, 6). Disturbances were shown to change network composition (Berga et al., 2012) and tillage was shown to affect Actinobacteria due to their hyphal growth form (Cania et al., 2020). Very interesting is that this drop is much lower in biocrusts, which is why we conclude that Actinobacteria are somehow protected against tillage in biocrusts. Important functions covered by Actinobacteria might be their ability to fight fungal infections in crops (Barka et al., 2016) or the degradation of pesticides (Mawang et al., 2021). Further functional analysis would be necessary to answer questions about their beneficial potential in agriculture. Nevertheless, this finding opens the debate for further research about the potential and risks of biocrusts in agricultural ecosystems. Networks in Fruchterman-Reingold layout based on significant ASV correlation analysis for bare soil and biocrust for each management variation (n = 12). Nodes are colored by phylum and the color code is given below the network figures, where also the abundance of phyla within the networks is given. Red connections display negative and green connections display positive correlations. The number of correlations defines the size of one node. The lower the p value of the correlation, the thicker the line between two nodes.
Frontiers in Microbiology 12 frontiersin.org Despite the lower complexity of biocrust correlation networks, the share of co-occurrences was higher in biocrusts compared to bare soil. The polymeric matrix of biocrusts facilitates not only nutrient trapping but also the exchange among microbes by, for example, facilitating movement and horizontal gene transfer interaction among microbes (Costa et al., 2018;Rossi et al., 2018;Cania et al., 2019a). Moreover, the complexity of extracellular polysaccharides and proteins promotes the association of microbes, which can decompose polymers . This is further underlined in our study by the finding that the highest frequency of co-occurrences is in the +org treatment, which adds further complex and recalcitrant carbon compounds to the biocrusts (Wang et al., 2017). Indeed, it could be demonstrated earlier that addition of organic treatment increases network complexity in soil (Schmid et al., 2018). The organic amendment specifically increased co-occurrences of bacterial-fungal partners (Figures 5, 6) like Alphaproteobacteria, Dothideomycetes, and Pezizomycotina, which had been identified as keystone taxa in organic matter decomposition previously (Hartmann et al., 2015).

Differential response of Cyanobacteria in bare soil and biocrusts
Cyanobacteria belonged to the dominating Bacteria (Figures 4,  5)-surprisingly, not only in the biocrusts but also in bare soils. This contradicts other agricultural studies where they have not been detected or only displayed a minor proportion of the bulk soil community (Uksa et al., 2015;Cano-Díaz et al., 2019). In particular, Microcoleus was one of the abundant observed genera in our study. Microcoleus is known as a key taxon for biocrust formation in drylands due to bundle sheath and the production of sticky extracellular polymeric substances (Garcia-Pichel and Wojciechowski, 2009;Pereira et al., 2013) and was shown to attract copiotrophic microbes by releasing photosynthesized carbon into the cyanosphere (Couradeau et al., 2019). This was further accompanied by chlorophyll a concentration in the bare soils in the same range as detected for Cyanobacteria-dominated biocrusts in arid areas (Román et al., 2019). Also, we identified Cyanobacteria as network hubs with a proportion of up to 12.5% of all network edges in bare soil networks (Supplementary Table S12), which suggests an early crustal stage, as observed in other early-stage biocrusts (Chilton et al., 2018;Maier et al., 2018). These connections might be beneficial for copiotrophic Bacteria (Burkholderiaceae, Rhodobacteraceae) and Fungi (Cucurbitariaceae), when no additional organic carbon was provided by crop residues (Couradeau et al., 2019;Pombubpa et al., 2020). In contrast, bare soil treatments with organic fertilization promoted negative correlations between Cyanobacteria and Gammaproteobacteria, Dothideo-and Sordariomycetes (Mueller et al., 2015;Wen et al., 2020), which might be attributed to the growth of heterotrophic Bacteria responding to the carbon input (Fierer et al., 2007).
Interestingly, the biocrusts of the cT treatment and the bare soils show similar network pattern with highest cyanobacterial relative abundance and connectivity (Figures 4-6; Supplementary Table S10). As revealed by different studies, soil surface disturbance, like tillage, is a critical factor for biocrust development and sets biocrusts back to an initial development stage, where Cyanobacteria play a crucial role (Lange et al., 1997;Kuske et al., 2012;Bu et al., 2013;Ferrenberg et al., 2015;Steven et al., 2015;Belnap et al., 2016;Maier et al., 2018). Biocrusts that were destroyed by soil turning due to tillage before seeding, quickly re-establish with a fast succession. This regeneration of biocrusts after soil surface disturbance observed in this study Correlation partners on class level based on the network analysis (n = 12) are shown as shares on total edges (in %) for bare soils and biocrusts for each management variation [mineral fertilization amount (120 vs. 240 kg N/ha·a), organic treatment (−org vs. +org), and tillage (conventional (cT) vs. reduced (rT) tillage)]. Red connections display negative and green connections display positive correlation. Microbiology   13 frontiersin.org   happened much faster than observed in (semi-) arid areas , where recovery rarely takes less than 5 years, and can need several decades depending on the type of disturbance. Nevertheless, to disentangle specific patterns of biocrust establishment on mesic, managed ecosystems, future studies need to analyze disturbance effects in more detail, in more systems and throughout the year.

Conclusion
Our study demonstrates that despite, and sometimes because of, intensive management like tillage and fertilization, biocrusts develop on agricultural fields and build up substantial phototrophic biomass over one growing season. Increased retention of nutrients and watermediated by biocrusts might improve water and nutrient availability at the beginning of the crop-growing season, possibly with positive feedback on crop growth. This was associated by a reduction of microbial diversity and the promotion of cross-kingdom co-occurrences, indicating that biocrusts become an additional hotspot for activity in soils comparable to rhizosphere or drilosphere. Surprisingly, Cyanobacteria played a negligible role in biocrust networks in autumn before harvest but dominated networks in bare soil communities. Thus, we conclude that they may serve as a seed for biocrust formation, at least in agricultural soils with low quality in terms of organic matter content or waterholding capacity. However, the effect of tillage on the biocrust formation pattern and biocrust properties requires further investigation. Future studies should be accompanied by measurements recording seasonal changes or resilience measurements toward short-term effects like summer drought or heavy rain events.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article, section "Material and Methods". The raw data for soil properties generated and analyzed for this study will be made available by the corresponding author upon request.

Author contributions
SS, MS, and UK contributed to the conceptualization and design of this study. MAr manages the fields of LUFA Speyer and assisted in sampling preparation and performance. Sampling was performed by SS. Material preparation, data collection, and analysis were performed by JK. Sequencing was performed by SK. The pipeline for quality filtering of sequencing and taxonomic assignment was provided by CS. The pipeline for network analysis was provided by GV. The first version of the manuscript was written by JK and all authors commented on previous versions of the manuscript. All authors contributed to the article and approved the submitted version.

Funding
This work has been supported by DFG Nachwuchsakademie (SCHU 2907/4-1) and funded additionally by the BMBF BonaRes project Inplamint (031B1062C). We acknowledge support from the Open Access Publishing Fund of the Technical University Munich.