Effects of soil environmental changes accompanying soil erosion on the soil prokaryotes and fungi of cool temperate forests in Southern Japan

ABSTRACT Soil erosion, which involves the degradation of the physical and chemical properties of soil, is a major threat to the soil environment. Although the effects of soil erosion on the physical or chemical properties of forests have been studied, little has been reported on the soil microbial community, which is likely to affect forest ecosystems. This study aimed to elucidate how the microbiome changed with the soil environment accompanying soil erosion in cool temperate mixed forests in Southern Japan, where soil erosion has been accelerated by the increased population of sika deer. We investigated the soil microbial communities of the different soil erosion intensities at three forest sites. In prokaryotic communities, diversity indices were increased with the sum of the height of exposed roots (SUMH), an index of soil erosion. In fungal communities, the relative abundances of plant pathogenic and wood saprotroph fungi were increased as SUMH increased and those of symbiotrophs and ectomycorrhizal fungi were increased with humus organic matter content, suggesting that the difficulties in establishing plants would be increased as soil erosion progressed because of the changes in the composition and function of fungal communities in eroded areas. Moreover, soil fungal communities had a more complex co-occurrence network than that of prokaryote, suggesting that the effects of soil erosion on fungal network is smaller than those on bacterial network. Changes in the soil environment induced by soil erosion altered the microbiomes in the deeper layers of the soil and had different effect on prokaryotes and fungi.


Introduction
Soil erosion has impacted more than 80% of the world's land and is a major threat to the soil environment (Borrelli et al. 2017).During erosion, litter and the soil surface layer rich in organic matter are runoff, exposing a deeper layer of soil; consequently, the soil environment and ecosystem functions in the region are often damaged (Guerra 1994;Gachene et al. 1997;Li et al. 2021).It has been reported that soil erosion alters the carbon cycle (Lal and Pimentel 2008), including the primary productivity in soil (Gregorich et al. 1998); carbon storage, such as soil organic matter (Guerra 1994;Quinton et al. 2010) and the percentage of organic carbon (Gachene et al. 1997); nutrients (Quinton et al. 2010); aggregateassociated nitrogen (Li et al. 2022); pH (Gachene et al. 1997); reduces the ability of soils to retain water (Li et al. 2021); and impacts soil respiration (Yao et al. 2019).In the agroecosystem, the loss of soil directly affects crop yields and both these factors and the processes of erosion have been intensively studied (Atreya et al. 2006;Colacicco et al. 1989;Snelder and Bryan 1995;Bakker et al. 2007;Qiu et al. 2021).In forest ecosystems, soil movements are controlled by the vegetation covers of the upper canopy layer and understory, tree density, and the slope (Hattori et al. 1992;Nanko et al. 2006).In deciduous forests in Japan, understory vegetation loss and direct rainfall on soil have been attributed to overgrazing by sika deer (Cervus nippon Temminck) (Ohashi et al. 2007;Miura and Tokida 2009).A decrease in the understory vegetation in the forest increases the movement of soil and litter on the forest floor (Chu et al. 2010).Thus, soil erosion can be accelerated following the vegetation loss from overgrazing by sika deer.Although vegetation loss has a negative effect on the abundances of soil microorganisms, including microbes (Zhao et al. 2011;Katayama et al. 2023), little has been reported on the changes in the soil microbial community by the soil erosion or even on the vegetation loss in forest ecosystems.
Microbial communities in the soil play key roles in the ecosystem and its functions because the microbiome is closely linked to plants and ecosystem processes, such as litter decomposition and nutrient cycling, and their functionality increases as soil microbial biodiversity increases (Wagg et al. 2014).In agroecosystem and terrestrial ecosystems, soil erosion has a negative impact on soil bacterial diversity (Qiu et al. 2021) and fungal richness (Du et al. 2021).A previous study reported a positive relationship between soil multifunctionality and bacterial alpha diversity, which decreases as soil erosion progresses (Qiu et al. 2021).The loss of bacterial and fungal diversities, including phylogenetical diversity and species richness, is likely to reduce ecosystem multifunctionality and negatively affecting soil fertility and food production (Delgado-Baquerizo et al. 2016).A previous study on forest ecosystems revealed the effects of environmental changes in the forest on microbial community structure (Nakayama et al. 2019).As shown by these studies, microbial communities and diversity are important information for understanding soil functioning (Wagg et al. 2014;Hu et al. 2021;Qiu et al. 2021) and these communities can be altered by environmental changes, including soil erosion in the forest.
Soil microbial communities in forest ecosystems feature a high abundance of Acidobacteria, Actinobacteria, and Proteobacteria (Uroz et al. 2013).Soil depth is one of the important environmental factors affecting microbial relative abundance and diversity (Brewer et al. 2019;Hao et al. 2020;Nash et al. 2022).The relative abundance of some oligotrophic bacteria, such as Actinobacteria, Chloroflexi, Nitrospirae, and the candidate phyla AD3 and GAL15, consistently increases as soil depth increases (Brewer et al. 2019;Hao et al. 2020).Ectomycorrhizal fungi that help host trees with nutrient uptake (Allen et al. 2003) have lower colonization with soil depth (Nash et al. 2022).These reports indicate that the exposure of deeper soil layer by soil erosion could alter the relative abundance of microbes in the soil surface soil.In addition, the microbial co-occurrence network, which is an important index for understanding the microbial community assembly and stability, may also be changed by soil erosion.The complexity of the fungal network is increased by water erosion in terrestrial ecosystems (Du et al. 2021), whereas a lower network complexity of bacterial network is observed in an eroded agroecosystem (Qiu et al. 2021).There are few studies of the microbial co-occurrence network in the forest and one of these studies (Nakayama et al. 2019) shows that the networks in a natural forest are more complex than those of an artificial forest.Thus, soil environmental changes caused by soil erosion can alter the microbial relative abundance and co-occurrence network in forest ecosystems; however, these changes remain unknown.
The purpose of this study is to elucidate the effects of soil environmental changes accompanying soil erosion on soil prokaryotic and fungal communities and their functionalities in forest ecosystems.To achieve this, we investigated: 1) microbial community, diversity and multifunctionality with the soil environmental changes accompanying soil erosion; 2) environmental factors affecting microbial communities and function; and 3) the co-occurrence networks of microbial communities under different soil erosion intensities.

Study site and soil sampling
This study was conducted at three forest sites in the cool temperate mixed forests in southern Kyushu, Japan: two sites, Sanpo [32°22'N, 131°11'E, 1400 m above sea level (a.s.l.)] and Maruju [32°22ʹN, 131°08ʹE, 1180 m a.s.l.] located in the Shiiba Research Forest of Kyushu University, in Miyazaki Prefecture, and one site, Shiraga [32°09'N, 130° 55'E, 1250 m a.s.l.] located in the national forest, in Kumamoto Prefecture.Soil sampling for this study was conducted in October 2021.To investigate the impact of soil erosion under the relatively homogeneous conditions in the field, we selected 16 individuals of Japanese beech (Fagus crenata (Fagaceae)) at each forest site and assessed the environmental changes around the tree.After removal of the A 0 layer, soil at 10 cm depth was collected at a distance of 1-2 m parallel to the slope from the selected tree.Collection points were within the canopy coverages of the selected tree.Soils were sieved through a 4 mm mesh and kept at −20°C until analysis.Prior to the soil sampling for microbial community analysis, environmental variables, including the canopy tree properties and soil properties, were measured for the same tree by Katayama et al. (2023) in June 2021.Three parameters were measured for the tree properties: diameter at breast height (DBH, cm), tree height (m), and leaf area index (LAI, m 2 m −2 ).Fourteen parameters were measured for soil properties: humus mass (Hmass, g m −2 ), litter mass (Lmass, g m −2 ), soil C (%), soil N (%), soil C/N ratio, soil pH (H 2 O), humus organic matter content (HOMcontent), mineral soil organic matter content (SOMcontent), humus organic matter mass (HOMmass, g m −2 ), bulk density (g cm −3 ), root biomass (g 100 cm −3 ), slope degree (slope, °) and sum and maximum values of the exposed root height from the ground to the root surface (SUMH and MAXH, respectively, cm) (Katayama et al. 2023).Humus and litter were collected from a 20 cm × 20 cm area at each sampling point and humus mass and litter mass were measured after drying at 70°C for 48 h.Soil C, N, and C/N ratio, pH, HOMcontent, SOMcontent, and HOMmass were analyzed using the soil under the A 0 horizon at a depth of 10 cm.Bulk density and root biomass were quantified using a 100 cc soil core sampler after drying at 105°C for 24 h.SUMH and MAXH were calculated using the exposed roots in a 0.7 m × 0.7 m area near the sampling points.In previous studies, the height of the exposed root was often used to estimate erosion intensity (Chartier et al. 2009;Osterkamp et al. 2012) and these indices include the soil erosion process, which comprises soil runoff to a lower altitude and deposition from a higher altitude; large values might indicate that the surface soil has been physically run off.At our sampling points, SUMH and MAXH values were different among the trees within same site and among the three sites.The order of the values of three sites was as follows: Shiraga > Sanpo > Maruju (Figure 1; Katayama et al. 2023).The details of soil sampling design and analysis methods are presented in Katayama et al. (2023).In this study, we used the 17 environmental variables described to analyze the microbial community.the accession number DRA015283.Raw fastq files were filtered using FASTX-Toolkit (ver.0.0.14)(Hannon 2009) to remove Illumina adapters.Quality read filtering was performed using sickle (ver. 1 33) (Joshi and Fass 2011) to trim 3' or 5' ends with low-quality scores.For 16S data, sequences shorter than 130 bases were discarded; for ITS data, sequences shorter than 50 bases were discarded.Then, the paired-ends were merged using FLASH (ver.1.2.11) (Magoč and Salzberg 2011).These sequences were imported to Qiime2 (ver.2021.11)(https://qiime2.org/)for bioinformatics analyses (Bolyen et al. 2019).The qiime2-dada2 plugin was used for denoising and the removal of chimeras (Callahan et al. 2016).

Taxonomic assignment, diversity indices and functional profile annotation of prokaryotic and fungal communities
Taxonomic assignment was performed for the amplicon sequence variants (ASVs) using the qiime2-featureclassifier (Bokulich et al. 2018) by the Greengenes prokaryotic 16S database (version 13.8) and the UNITE fungal ITS database (version 8.2).A database of representative sequences from operational taxonomic units was generated through clustering of Greengenes or UNITE sequences at 97% identity.
Alpha diversity measures were assessed on the ASVs abundances and abundances of ASVs rarefied to an even sequence depth of 17,817 reads for V3V4 and 22,304 reads for ITS1 per sample in Qiime2.Alpha diversity indices, including Chao1, evenness, Faith's phylogenetic diversity (Faith's PD) and Shannon index were calculated using Qiime2.
For the functional profiling of prokaryotic communities, representative sequences that were outputted by Qiime2 using the raw sequence data and PICRUSt2 (Douglas et al. 2020) were utilized to predict the metagenomic functions of the prokaryotic communities via marker gene profiles (Liu et al. 2021;Chen et al. 2022).Before linear-mixed-effects model analysis, centered log-ratio (clr) translation was used for standardization.To represent the multifunctionality of each sample, Z scores of the results of PICRUSt2 (MetaCyc pathway abundances were used) were calculated.The relative abundance of the results of PICRUSt2 was translated to a Z-score using the value of the standard deviation from the mean.The multifunctionality index was calculated as the average Z-score for all pathways measured in each sample.For the functional profiling of fungal communities, rarefied data and FUNGuild tools were used to annotate fungal functions (Nguyen et al. 2016).The fungal functional annotations, including trophic mode, trait, and guild, were assigned for each sample and only confidence scores of "Probable" and "Highly Probable" were used for the subsequent linear-mixed-effects model analysis.

Statistical analysis
The similarity of microbial communities among each sample was analyzed by non-metric multidimensional scaling (NMDS).Sample variation of ASVs was represented by an ordination using an Aitchison distance matrix.The influence of environmental variables on the microbial community structure was tested by the envfit (permutations = 9999) function in the R package vegan (Oksanen et al. 2022).From NMDS envfit function results, significant variables with adjusted p-values of <0.05 were primarily selected.To assess correlations between geographic distance or environmental variables and community composition, Mantel test was run in R package vegan with the threshold of p < 0.05 (Zhang et al. 2014).We calculated the correlation coefficients between community composition and geographical distance estimated from GPS data and selected environmental variables in envfit functions.Then, only the remaining environmental variables with p < 0.05 in the Mantel test were used for further analyses (for details, see the results section).To avoid the multicollinearity for the subsequent linear regression analyses, variables with same or reverse directions in NMDS plots and large correlation coefficient (details in next

Maruju
Sanpo Shiraga paragraph, r > 0.5 or r < −0.5 with adjust p < 0.05) between the two variables were discarded and one variable, with the largest r in the Mantel test, remained.Based on the above analyses, the environmental variables of SUMH, HOMcontent, and pH were selected for both prokaryotic and fungal community analyses.The correlation between environmental variables were calculated using the Spearman's rank correlation coefficient test in the R package psych (Revelle 2017), and "'pairwise'" functions were used to handle missing data.The Benjamini -Hochberg false discovery rate was used to adjust the p-values in the correlation (Benjamini et al. 2006).The results of the correlation matrix were checked by a heatmap graph using the R package pheatmap (Kolde 2019) and corrplot (Wei and Simko 2021).
To test the effects of selected environmental variables on prokaryotic or fungal alpha diversity and the relative abundance of both communities at the taxon levels, linear-mixedeffects models were generated using the R package lmerTest (Kuznetsova et al. 2017).Four alpha diversity indices (Chao1, evenness, Faith's PD and Shannon index) or rarefied prokaryotic or fungal relative abundance (reads/sample) of each taxon were set as response variables and selected environmental variables were set as fixed effects with the site term as random effects.The R package ggplot2 (Wickham 2016) was used for visualization.The effects of prokaryotic alpha diversity on fungal alpha diversity were also analyzed by linear-mixed-effects models.Regarding microbial function, the effects of environmental factors on the multifunctionality and relative abundance of fungal functional categories were also tested by the linear-mixed-effects model.In each model, the average Z-score or relative abundance of fungal functional category was set as the response variable and selected environmental factors were set as fixed effects with the site term as random effects.
To understand the interactions of microbiomes, correlation relationships between the relative abundance of prokaryotes and fungi were calculated using the Spearman's rank correlation coefficient test in the R package psych (method details are presented in the second paragraph of Section 2.4).For further analysis, co-occurrence networks were constructed using the R package psych (Revelle 2017), Hmisc, and igraph (Csardi and Nepusz 2006) based on the Spearman's correlation matrixes.In this analysis, all 48 samples were separated to 3 sites; then, each site of 16 samples was separated again into two batches of eight samples depending on SUMH values (eight samples with highest SUMH values (highly eroded samples) and the remaining eight samples with lowest SUMH (medium-or light-eroded samples)).Only ASVs that occurred with at least 0.1% relative abundance in all eight samples were included in the analysis.Centered log-ratio (clr) translation was used for ASVs abundance treatment, and co-occurrence networks for prokaryotes and fungi were constructed based on the Spearman's correlation using only significant correlations (adjusted p < 0.05) of r > 0.6 or r < −0.6 in the R packages NetCoMi (Peschel et al. 2021) and psych (Revelle 2017) calculated within prokaryotes and fungi, respectively.For prokaryotic -fungal co-occurrence networks, prokaryotic and fungal clr-translated ASVs abundance data were combined and the correlations were then calculated together.The network parameters of the numbers of connected nodes and edges were used to assess soil microbial network complexity.
Node-level topological features, including the degree, eigenvector centrality, and betweenness centrality of our networks were examined.The non-parametric Kruskal -Wallis analysis combined with Dunn's multiple comparison tests was used to test three features in the R packages PMCMR and PMCMRplus (Pohlert 2014(Pohlert , 2022)).The p-values were adjusted using the Bonferroni method.The Louvain method was used for modularity analysis; the eight highest proportions of modularity classes were colored in the plot, and the results of the correlation matrix were visualized using Gephi software (https://gephi.org;Bastian et al. 2009).
NMDS ordination plots showed that prokaryotic and fungal communities were different among the three sites (prokaryotes: R 2 (site) = 0.12, p < 0.001 (PERMANOVA); fungi: R 2 (site) = 0.10, p < 0.001 (PERMANOVA)) and the several variables were significantly (adjusted p < 0.05) associated with prokaryotic and/or fungal communities (Figure 2 and Supplementary Table S1).Based on the distance matrix of the significant environmental variables in NMDS, Mantel tests showed significant correlations (p < 0.05) between selected environmental parameters and microbial communities (Table 1).In particular, SUMH, soil pH, and HOMcontent were significantly related in both prokaryotic and fungal communities.Correlations of HOMcontent (r = 0.39) and all environmental variables that included all significant variables by NMDS (r = 0.36) were higher than geographical distance (r = 0.31) in prokaryotic communities, whereas correlations of all environmental variables (r = 0.20) were lower than geological distance (r = 0.34) in fungal communities.
To avoid the multicollinearity for linear-mixed-effects model analyses, the environmental variables were further selected based on the correlation coefficient and direction of the NMDS plots.In the NMDS plot of prokaryotic communities, the arrow of HOMcontent was in the same direction as the soil C, SOMcontent, root biomass, and slope, and in an opposite direction to bulk density.Because HOMcontent had highest r value in Mantel tests of prokaryotic communities and these four variables had higher correlation coefficient values (absolute value > 0.5) with HOMcontent (Supplementary Figure S1), these four variables were discarded.In the NMDS plot of fungal communities, SUMH, HOMcontent, and soil pH, which were indicated to be significantly correlated with fungal communities by Mantel tests and had different directions in the plot.Thus, SUMH, HOMcontent, and soil pH were selected for the subsequent linear-mixed-effects model analyses in both prokaryotic and fungal communities.The Mantel test was conducted again using the set of SUMH, HOMcontent, and soil pH in prokaryotic and fungal communities and demonstrated that correlation coefficients values (0.49 (p < 0.001) in prokaryotes and 0.28 (p = 0.004) in fungi) were higher than those of all environmental variables (Table 1).

Effects of the soil environmental variables on the alpha diversity of soil microbial communities and relative abundances of each microbial taxon
In prokaryotic composition, Acidobacteria and Proteobacteria showed high relative abundance; in fungal composition, Ascomycota and Basidiomycota presented high relative abundance for all the samples (Supplementary Figure S2).To understand how soil environmental variables affect the alpha diversity of soil microbial communities, a linear-mixed-effects model analysis was conducted (Table 2, Supplementary Figure S3).In prokaryotic communities, the values of Faith's PD significantly increased with SUMH (p = 0.023) and the values of Chao1 diversity significantly decreased as HOMcontent increased (p = 0.036).The evenness and Shannon index of prokaryotic communities did not change with the environmental variables.In contrast, in fungal communities, the values of evenness diversity significantly decreased as HOMcontent increased (p = 0.037) and the value of Shannon diversity followed a downward trend as HOMcontent increased (p = 0.064).Chao1 and Faith's PD were not affected by the environmental variables.Neither diversity index was affected by pH (p > 0.100).Moreover, a linear-mixed-effects model was used to analyze the prokaryotic and fungal communities' diversity indices to elucidate whether the prokaryotic diversity affected fungal diversity.Although there was no significant effect, prokaryotic Faith's PD followed an upward trend with fungal Faith's PD (p = 0.060), and the same trend was also found for prokaryotic and fungal Shannon diversity (p = 0.053) (Supplementary Table S2).
To understand the effects of soil environmental variables on microbial communities at the taxon level, further linearmixed-effects model analysis was conducted.The results of the analysis using prokaryotic taxa as response variables are  shown in Table 3. Relative abundances (reads/sample) of the candidate phyla Dormibacteraeota (AD3), GAL15, and Nitrospirae were significantly increased as SUMH increased, whereas that of Actinobacteria was significantly decreased as SUMH increased.The relative abundance of Chlamydiae, Cyanobacteria, and TM6 significantly increased as HOMcontent increased, whereas that of Chloroflexi, FCPU426, Firmicutes, and WPS-2 significantly decreased as HOMcontent increased.The relative abundance of Chlorobi significantly increased as soil pH increased, but that of Euryarchaeota, Chlamydiae, Verrucomicrobia, and WPS-2 significantly decreased as soil pH increased.In contrast to the significant changes in prokaryotic composition, there was no significant effect of environmental variables on fungal composition, although some marginal significance (p < 0.100) was found (Table 4).

Effect of soil environmental variables on the prokaryotic and fungal functional groups
The relationship between microbial predicted function and environmental variables was investigated by linear-mixedeffects model analyses (Tables 5 and 6).In prokaryote, 53 of level 2 and 447 of level 3 of KEGG categories were identified by PICRUSt2.Five level 2 KEGG categories concerning metabolism significantly increased as SUMH increased, and five categories concerning metabolism significantly increased as HOMcontent increased, whereas eight categories concerning metabolism significantly decreased as pH increased (Table 5).The multifunctionality index significantly increased as Chao1 (estimate = 0.0006, p = 0.030) or Faith's PD (estimate = 0.0119, p < 0.001) increased (Supplementary Table S3).
The effect of soil environmental variables on the fungal functional grouping is presented in Table 6.The fungal composition, with known functional groups based on the FUNGuild, comprised over 45% of the total sequence reads.There are seven types of trophic mode, six types of trait, and sixty types of guild.The relative abundances of plant pathogenic and wood saprotroph fungi were significantly increased as SUMH increased.The relative abundances of symbiotroph and ectomycorrhizal fungi were significantly increased as HOMcontent increased, but the relative abundances of saprotrophicsymbiotrophic and pathotrophic fungi, as well as soft rot fungi, were significantly decreased as HOMcontent increased.The relative abundances of wood saprotrophic, animal pathogen -endophyte-epiphyte -fungal parasiteplant pathogenic -wood saprotroph and plant pathogenic -wood saprotroph fungi were also significantly decreased as HOMcontent increased.There was no significant relationship between fungal predicted functions and soil pH (Table 6).

Co-occurrence network with the changes in SUMH
To understand the interaction between prokaryotes and fungi with changes in SUMH change, a Spearman's correlation heatmap was generated to analyze composition of prokaryotic and fungal communities at the phylum level (Figure 3).The correlation analysis showed that the Mucoromycota phylum was negatively correlated with seven phyla of prokaryote (WPS-2, GAL15, Euryarchaeota, Crenarchaeota, Chlorobi, AD3, and Gemmatimonadetes), but positively correlated with two phyla (Chlamydiae and Planctomycetes).In contrast, the Chytridiomycota phylum was positively correlated with six phyla of prokaryote (Euryarchaeota, Crenarchaeota, Nitrospirae, Chlorobi, AD3, and Chloroflexi), but negatively correlated with one phylum of bacteria (Planctomycetes).
For a further understanding of interactions within prokaryotic or fungal communities in different SUMH conditions, co-occurrence networks were constructed using high (eight highest samples) and low (eight lowest samples) SUMH in each site.The network of prokaryotic and fungal communities demonstrated distinct co-occurrence patterns with different SUMH values (Figure 4).The numbers of connected nodes and edges in the fungal networks were higher than those of prokaryotic networks, except for the low SUMH of the Maruju site.Additionally, the complexity of fungal networks in the high SUMH groups was higher than that in low SUMH groups at all sites.In the prokaryotic network, the numbers of connected nodes and edges were smaller in the high SUMH groups than in low SUMH groups, except for the Sanpo site.We examined node-level topological features, the degree, eigenvector centrality, and betweenness centrality of the obtained networks (Supplementary Figure S4).In the analysis of prokaryotic networks, the Kruskal -Wallis test for all groups revealed significant differences between high and low SUMH groups, as well as among sampling sites following H statistics: H = 109.08(p < 0.001) for the degree, 95.41 (p < 0.001) for eigenvector centrality, and 54.85 (p < 0.001) for betweenness centrality.The Kruskal -Wallis test for all groups of fungal networks showed lower H statistics than that of prokaryotic networks (H = 24.13(p < 0.001) for the degree, 2.61 (p = 0.760) for eigenvector centrality, and 38.35 (p < 0.001) for betweenness centrality).For a modular analysis of the co-occurrence network, the eight largest modules, including the highest number of nodes in each group, were selected for calculating node numbers in different prokaryotic and fungal class (Supplementary Tables S4 and S5).In prokaryotic networks, more than 10% of node numbers for class Acidobacteriia, Alphaproteobacteria, and DA052 were observed as major classes in all groups.In fungal networks, Agaricomycetes was observed to have the highest proportion of node number classes (more than 35%) in all groups.
In prokaryotic -fungal co-occurrence networks, the numbers of all connected nodes and all edges in high SUMH in Maruju were lower than that for the low SUMH values, but the opposite trends occurred for Sanpo and Shiraga (Supplementary Table S6).The number of fungiprokaryotes edges was lower for high SUMH values in Maruju and Shiraga than those for low SUMH values, whereas the value was higher for high SUMH values in Sanpo than that for low SUMH values.The proportions of fungi -prokaryotes edges of all edges were lower than 15% in all groups except the Shiraga low SUMH values group.

Discussion
Little has been reported on the impact of soil erosion on the microbial community in forests, which could provide crucial effects on the microbial functionality and forest ecosystems.This study investigated how the soil environmental changes accompanying soil erosion impacted both prokaryotic and fungal communities by analyzing samples collected from the soils with different erosion intensities at three forest sites.According to the NMDS analysis and the Mantel test in both  prokaryotic and fungal communities, three environmental variables (SUMH, HOMcontent, and soil pH) were selected as the factors affecting the communities (Table 1, Supplementary Table S1).In NMDS plots, the SUMH, which was an index of the erosion, caused the sample variability along the arrow direction, and this trend was confirmed in the samples of three sites (Figure 2), indicating that root exposure caused similar changes in the microbial community in three sites.HOMcontent, which is an index of the organic matter content in the sampling points, also caused the samples to vary among sampling points at three sites.Soil pH followed the trend in between the above two variables.Thus, hereafter, we discussed the effect of SUMH, which is a physical change in the erosion environment, and those of HOMcontent and soil pH, which are chemical changes.

Effect of soil erosion on the prokaryotic community and its functionality
The present study showed that prokaryotic species diversity was increased as SUMH increased and HOMcontent decreased (Table 2), indicating that soil erosion could increase prokaryotic diversity in forest ecosystems.
A previous study in the agroecosystem has reported that soil bacterial diversity significantly decreased along the soil erosion intensities, whereas the lightly eroded areas had higher bacterial diversities than non-eroded areas (Qiu et al. 2021).This indicates that soil bacterial diversity responds to the disturbances differently, and that the peak would occur at the intermediate disturbance levels suggested by the intermediate hypothesis (Connell 1978;Santillan et al. 2019).Our study was conducted in forest ecosystems and whether the soil erosion accompanies the diversity of soil prokaryote in the same way as agroecosystems is unknown, but the increase in prokaryotic species diversity with SUMH will support this hypothesis.Changes in the relative abundance of some oligotrophic bacteria, namely deeper soil-inhabiting bacteria (Table 3), were likely caused by deeper soil layer exposure as soil erosion progressed.SUMH had positive effects on the relative abundance of AD3, Nitrospirae, and GAL15; however, it was not associated with the relative abundance of Acidobacteria and Euryarchaeota (Table 3).AD3, Chloroflexi, Nitrospirae, Euryarchaeota, and GAL15 are considered as adapting to lower nutrient conditions, including resource-limited deeper soil (Krzmarzick et al. 2012;Hug et al. 2013;Ho et al. 2017;Brewer et al. 2019;Hao et al. 2020).Acidobacteria is also considered to be oligotrophic bacteria (Ho et al. 2017); however, the lower relative abundance of Acidobacteria has been reported in deeper soils, which might be caused by increased soil pH owing to the decay of organic matter (Jones et   increased, but unchanged with SUMH (Table 3).Some Chloroflexi also have an oligotrophic property inhabiting deeper soil (Krzmarzick et al. 2012;Hug et al. 2013;Ho et al. 2017;Brewer et al. 2019;Hao et al. 2020) and, in forest ecosystems, Chloroflexi is scarce at high levels of soil nutrient substrates (Yang et al. 2022).Therefore, our results suggest that Chloroflexi could be affected by the runoff of the organic matter of the surface soil along the soil erosion.Given these results, the composition of soil prokaryote has been changed by soil erosion through surface soil and organic matter removal and deeper layer soil exposures.Soil pH is also an important factor in shaping bacterial composition (Wang et al. 2019;Ren et al. 2020).Our results showed that the relative abundance of Chlamydiae and WPS-2 decreased as the pH increased, whereas the relative abundance of Chlorobi increased as the pH increased (Table 3).A previous study supported our results, reporting that WPS-2 has a high relative abundance in forest bare acidic soils (Sheremet et al. 2020).However, other studies report that there is no correlation between phylum Chlamydiae or WPS-2 and pH (Wang et al. 2019;Ren et al. 2020) and a negative correlation between Chlorobi and pH in agricultural soil (Wang et al. 2019).Some studies report that soil pH alone might not control the bacterial community composition in forest soil (Cho et al. 2016;Liu et al. 2020) and that bacterial diversity and community composition are controlled mainly by ecosystem type rather than geography (Fierer and Jackson 2006); our results of the above-mentioned bacterial abundance might be specific results to the forest ecosystem or impacted through the other soil properties.
The present study showed that multifunctionality had a positive linear relationship with prokaryotic diversity (Supplementary Table S3), implying that soil erosion increased multifunctionality.A previous study of agroecosystems reported that soil erosion significantly reduced bacterial diversity and functionality (Qiu et al. 2021).Another study on water erosion showed that the erosional sites had lower enzyme activities and microbial biomass compared with depositional sites (Li et al. 2015).These previous studies were not conducted in a forest ecosystem and the results were not comparable with our study, there was a significant difference in the soil pH change during erosion.In our study, SUMH was positively correlated with soil pH (Supplementary Figure S1), but soil pH change was not observed by Qiu et al. (2021).Considering these results, prokaryotic multifunctionality may be affected by erosion in different ways in different ecosystems and, in the forest ecosystem, erosion-induced variable change, such as pH may be one of the factors.
Our results showed the positive effects of SUMH and HOMcontent and the negative effect of soil pH on several abundances of soil prokaryotic metabolic pathways (Table 5).The abundances of amino acid, carbohydrate, or lipid metabolism pathways were positively related to HOMcontent; however, there was no significant relationship with SUMH.A previous study reported that the abundances of carbohydrates and amino acid metabolism pathways decreased in deeper soil, whereas the abundances of lipid and energy metabolism pathways increased in deeper soil (Koner et al. 2022), indicating that the effects of organic matter loss in this study were similar to the changes in the abundances of metabolism pathways in deeper soil (Koner et al. 2022).Our study also found that soil pH was also an important factor in the abundances of metabolism pathway changes.A previous study showed that soil pH was more important than soil organic carbon in changing prokaryotic functions in agricultural soils (Wang et al. 2019).Our results were partially similar to this previous study and the increase in soil pH during soil erosion decreased prokaryotic metabolic activities.There are limited studies on the changes in prokaryotic functions with the soil environmental variables and further analyses are required to understand the impact of soil erosion on prokaryotic ecosystems.

How soil erosion affects the fungal community and its functional groups
The effects of soil environmental variables on fungal diversity and community were clearly different to those on prokaryotic diversity and community.The present investigations on fungal diversity (Table 2) and the relative abundance of each phylum of soil fungal community (Table 4) revealed that humic organic matter content decreased the fungal community evenness, but the soil environmental variables did not change other diversity indices and the relative abundances of each phylum.A previous study has reported that soil fungal communities near roots were unaffected by environmental factors and were mainly influenced by geographic distance (Zhang et al. 2017), indicating that soil fungal communities in the forest have a high tolerance for soil environmental changes.Our results suggested that the fungal diversity indices of forest soil might be less affected by the soil erosion than the prokaryotic diversity indices.
Fungal functional groups have been affected by the SUMH and HOMcontent and the results indicated soil erosion would impact fungal communities in forest ecosystems.The relative abundance of the plant pathogen and wood saprotroph increased with an increase in SUMH, and pathotrophs, as well as soft rot fungi, were significantly decreased by HOMcontent (Table 6).Pathotrophic fungi are generally considered to perform negative effects on the performance of plants (Anthony et al. 2017) and the present results suggested that an increase in SUMH would increase the proportions of plant pathogenic fungi, and that less organic matter might increase fungi-induced negative effects on plants.Additionally, the relative abundances of symbiotrophs, such as ectomycorrhizal fungi, were increased as HOMcontent increased (Table 6).Ectomycorrhizal fungi actively decompose organic matter to acquire nitrogen in surface soil (Hobbie and Horton 2007;Lindahl et al. 2021) and transfer the nitrogen to the plant.Deeper soil also contained less organic matter (Edmondson et al. 2012), and the colonization of ectomycorrhizal fungi declined as soil depth increased (Nash et al. 2022).The elucidated negative relationship between HOMcontent and ectomycorrhizal fungi suggested that soil erosion would cause a reduction in organic matter and a decrease in the abundance of symbiotrophs, which could lead to negative effects on plant nutrient acquisitions, growth, reproduction, and seedling establishment.Given the above results, forest areas with little organic matter followed by soil erosion may have difficulties in plant seedling recruitment, survival and growth in the forest and further vegetation loss would be induced by the fungal compositional changes in the forest.

Network complexity of the microbiome at different SUMH conditions and differences between prokaryotes and fungi
The co-occurrence networks obtained in this study demonstrated that fungal networks' complexity was higher than prokaryotic networks' complexity at all sites (Figure 4).A cooccurrence network provides interpretable information about community stability and assembly (Freilich et al. 2018;de Vries et al. 2018;Wagg et al. 2019;Gao et al. 2022), although the environmental preferences of each taxon may possibly interfere with biotic interactions in the co-occurrence networks and the obtained network may not always demonstrate a true biotic association.A previous study reported that, during environmental changes such as drought, bacterial networks would be lost complexity whereas fungal networks would not be lost (de Vries et al. 2018).Studies about fungal communities also report that fungal communities are more resilient to environmental disturbances (Osburn et al. 2019) and that soil resource limitations may enhance fungal community stability and increase nutrient transfer efficiency that can be attributed to the increased network complexity (Morrien et al. 2017;Du et al. 2021).The present analyses of the relative abundances of each phylum also showed that the number of prokaryotes was more influenced by the soil environmental changes (SUMH and HOMcontent) than fungi (Tables 3 and 4).These results suggest that the effects of the soil erosion on the fungal networks were smaller than prokaryotic networks.
In the present module analysis, the bacterial classes Acidobacteriia, Alphaproteobacteria, and DA052 and the fungal class Agaricomycetes (Supplementary Tables S4 and S5) were the main members of the eight largest modules.The importance of these classes in network stability is also reported in previous studies.Class Alphaproteobacteria, belonging to the phylum Proteobacteria, is considered critical to the N transformation processes (Nakayama et al. 2021).Several classes of Agaricomycetes impact the composition of microbial and plant communities (Miyauchi et al. 2020) and these are the main microbiomes for maintaining network stability in the same ecosystem despite environmental changes (Chen et al. 2020(Chen et al. , 2020;;Ritter et al. 2021).Thus, the main classes revealed in the present module analysis might be important to sustain network stability during soil erosion in the forest.However, as we do not have further evidence to explain how microbiomes live in soil erosion areas, further studies such as metagenomic and metatranscriptomic analyses could give insight into the network stabilities.
Prokaryotic -fungal co-occurrence networks can be influenced by environmental changes, such as drought stress and the conversion from natural forest to plantation (de Vries et al. 2018;Nakayama et al. 2019); however, the present prokaryotic -fungal co-occurrence network analysis showed different responses at the three sites for high and low SUMH values (Supplementary Table S6).Recent studies have shown that microbial co-occurrence network connectivity influences the stability of microbial community (de Vries et al. 2018;Wagg et al. 2019;Gao et al. 2022), and decreasing the number of fungal and prokaryotic edges results in a loss of network complexity (de Vries et al. 2018;Nakayama et al. 2019).In our prokaryotic -fungal co-occurrence network, similar results were observed only in Maruju and values of nodes and edges in other sites were not lower for high SUMH values.The three sites have experienced different intensities of soil erosion because the ranges of SUMH and HOMcontent values were different (Supplementary Figure S3; Katayama et al. 2023) and the erosional history at each site may affect the prokaryotic -fungal co-occurrence network.The Maruju site had weakly experienced soil erosional changes because the SUMH value was lower than that in Sanpo and Shiraga and the HOMcontent was higher.Thus, the co-occurrence network differences between high and low SUMH value may represent the early stages of soil erosional changes in the network.Our study was conducted only once at the three different sites, thus, continuous monitoring is necessary to reveal the soil microbial diversity, function, and networks impacted by the soil environment changes in the different soil erosion stages.

Conclusions
In the present study, we investigated the impacts of soil environmental changes during the soil erosion on soil microbial communities in cool temperate forests in Southern Japan.SUMH, HOMcontent, and soil pH were the major drivers that changed microbial diversity, community, and network in this forest ecosystem.In particular, prokaryotic communities were changed significantly in the composition and function and relative abundance of some oligotrophic bacteria that were often observed in deeper soil.Those changes were associated with the SUMH value, suggesting that prokaryotic change was caused by deeper layer soil exposure as soil erosion progressed.However, soil fungal communities were less influenced by the soil environmental variables and had a more complex cooccurrence network than that of prokaryote, suggesting that the effects of soil erosion on fungal network was smaller than prokaryotic network in forest ecosystems.However, the functions of the fungal communities were altered in the soils, with the abundance of plant-pathogens and poor symbiotrophs in the highly eroded areas, indicating that soil erosion will prevent the seedling establishment and lead to further vegetation loss.These insights will be helpful to understand and conserve forest ecosystems during environmental changes.

Figure 1 .
Figure 1.Photos of three study sites.Upper photos present each view of forests.Lower photos demonstrated highly eroded areas around the roots of Japanese beech in each forest.

Figure 2 .
Figure 2. NMDS ordination plots of the prokaryotic (a) and fungal (b) communities of the three sites.Arrows indicate significant environmental variables as analyzed by the envfit function (adjusted p < 0.05).

Figure 3 .
Figure 3. Spearman's correlation heatmap between the relative abundance of the main phylum of fungi (row) and prokaryotes (column).Correlation coefficient values, r, are presented in different colors on the right side color range.Adjusted p-value: * < 0.05 and ** < 0.01.

Figure 4 .
Figure 4. Separate prokaryotic and fungal co-occurrence networks at three sites.Nodes represent individual ASVs; edges represent significant positive or negative Spearman's correlations (r > 0.6 or < −0.6 and adjusted p < 0.05).Different colors, including purple, green, blue, black, orange, red, dark green, and tan, represent the eight largest modules.

Table 1 .
Mantel test summary statistics for selection of environmental variables.

Table 2 .
Results of linear-mixed-effects model analysis of prokaryotic and fungal alpha diversity with three selected environmental variables.

Table 3 .
Results of linear-mixed-effects model analysis of the relative abundance of prokaryotic phyla for three selected environmental variables.

Table 4 .
Results of linear-mixed-effects model analysis of the relative abundance of fungal phyla for three selected environmental variables.

Table 5 .
Results of linear-mixed-effects model analysis of relative abundance of KEGG categories predicted by PICRUSt2 for three selected environmental variables.

Table 6 .
Results of linear-mixed-effects model analysis of the relative abundance of fungal functional group predicted by FUNGuild for three selected environmental variables.