*Article* **Effects of Trace Metals and Municipal Wastewater on the Ephemeroptera, Plecoptera, and Trichoptera of a Stream Community**

**Marek Let \*, Jan Cern ˇ ý, Petra Nováková, Filip Ložek and Martin Bláha \***

Faculty of Fisheries and Protection of Waters, South Bohemian Research Center of Aquaculture and Biodiversity of Hydrocenoses, University of South Bohemia in Cesk ˇ é Budˇejovice, Zátiší 728/II, 389 25 Vod ˇnany, Czech Republic; cernyj18@frov.jcu.cz (J.C.); novakovapetra@frov.jcu.cz (P.N.); lozekf@frov.jcu.cz (F.L.) ˇ **\*** Correspondence: mlet@frov.jcu.cz (M.L.); blaha@frov.jcu.cz (M.B.)

**Simple Summary:** Mayflies (Ephemeroptera), stoneflies (Plecoptera), and caddisflies (Trichoptera) (EPT) are aquatic insects that are well known to the general public and are commonly used as indicators of environmental quality in water management. Knowledge of how EPT communities react to human-induced gradients in real environments can be important, for example, during the assessment of the implications of newly planned or currently active human disturbances for natural or cultural landscapes. We sampled a stream ecosystem affected by mining and smelting industries and communal wastewaters with pronounced concentrations of cadmium, lead, and zinc, as well as high levels of pesticides, pharmaceuticals, illegal drugs, and sewage-derived organic matter. Changes in other environmental factors such as increases in temperature were also studied at the affected sites. The abundance and species richness of stoneflies fell rapidly at the study sites. The richness of mayfly families also declined, from four to one, even though overall mayfly abundance was not affected. Conversely, the abundances of caddisflies were higher at the affected sites, and their richness did not decrease. This study will provide feedback for ecotoxicologists who perform better controlled and manipulated tests in laboratories, although any such test results are limited by simple artificial environments.

**Abstract:** Abundances of EPT larvae sampled in a Central European locality affected by mining and smelting, as well as by the continual inflow of treated communal wastewaters (WWs), were recorded. High concentrations of trace metals in water (maximum 1200 <sup>μ</sup>g·L–1 for zinc) and sediments (maximum 140,000 mg·kg–1 in dry weight for lead) were found at the most contaminated sites. The highest loads of pesticides, pharmaceuticals, and illegal drugs were found under the WW effluent. Other associated factors such as the physicochemical parameters of the water and alterations to microhabitats were also evaluated and taken into account. Although EPT richness was lower at affected sites, abundances did not fall. Stoneflies were dominant at unaffected sites, while caddisflies dominated at affected sites. Only baetid mayflies were detected at the sites contaminated by trace metals and WWs; ephemerellid, heptageniid, and leptophlebiid mayflies were absent from these sites. The site contaminated by trace metals was also inhabited by numerous limnephilid caddisflies, in which limb malformations were detected in up to 11.8% of all specimens of a single taxon. Downstream from the entrance of the WWs, the locality was dominated by hydropsychid caddisflies. The increasing prevalence of predator or passive filter-feeding strategies in these EPT communities was significantly related to increasing water conductivity and acute ecosystemic exposure to 'poorly treated' WWs.

**Keywords:** anthropogenic disturbances; aquatic insect; environmental gradients; heavy metals; industrial pollution; wastewater treatment plant

**Citation:** Let, M.; Cerný, J.; ˇ Nováková, P.; Ložek, F.; Bláha, M. Effects of Trace Metals and Municipal Wastewater on the Ephemeroptera, Plecoptera, and Trichoptera of a Stream Community. *Biology* **2022**, *11*, 648. https://doi.org/10.3390/ biology11050648

Academic Editors: Daniel Puppe, Panayiotis Dimitrakopoulos and Baorong Lu

Received: 8 April 2022 Accepted: 21 April 2022 Published: 24 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

The aquatic insect orders of mayflies (Ephemeroptera), stoneflies (Plecoptera), and caddisflies (Trichoptera) (henceforth, EPT) are frequently used indicators of environmental quality, particularly in running waters [1,2]. Their larvae dominate in unpolluted headwaters where proportions of fine sediment deposits are low and are essential for correct nutrient flow cycling in these environments [3–5]. The contribution of EPT taxa adapted to conditions in downstream ecosystems (lower rhithral and potamal zones) is likewise important, although populations of these specialists, namely certain species of mayflies and stoneflies, are prone to be decimated by the cumulative effects of intensive human activities in lowland areas [6–8].

Given their rapid expansion through the natural environment, the ecological impact of human-induced gradients needs to be monitored, and possible threats are objectively predicted. The contamination of stream ecosystems by toxicants is a global issue; nevertheless, engineering work severely modifies the morphology of whole catchment areas and causes substantial changes in hydrological regimes, habitat structure, and water quality [8–10]. Hence, communities in running waters in cultural landscapes are habitually subjected to multiple anthropogenic factors [11]. Therefore, the manipulative testing of stressors within controlled conditions is now seen as increasingly important. Despite this, natural experiments (observational studies) are still needed to reflect the validity of derived results and vice versa [12].

The experiment described in this work explored an EPT community in a temperate zone in the European Central Highlands ecoregion [13], specifically, in a catchment area severely affected by industrial pollution (trace metals) and treated municipal wastewaters (WWs) whose levels are relatively high in terms of current European standards [14,15]. Measured confounding environmental variables (EVs)—concentrations of pesticides, physicochemical parameters of water, choriotopic composition, and current velocity—were included in the analyses of the changes in community composition. The shift in artificially/hierarchically assigned values for feeding strategies in the EPT community along with measured gradients of the EVs was also tested.

The contamination of streams by trace metals (often generalised as 'heavy metals') is a global phenomenon, and many studies tackling the responses of macroinvertebrate communities have identified mining as the cause of this type of pollution [16–19]. Increased concentrations of cadmium, copper, lead, and zinc—the so-called 'bivalent metals' that bind to both sulphur/nitrogen and oxygen groups—are most often reported as the cause of shifts in stream macroinvertebrate compositions in real environments or as the potential origin of impaired bottom-up control in ecosystems [20,21]. It is easy to misinterpret an oversimple comparison of concentrations measured in a real environment with the lethal values estimated for several EPT taxa under controlled conditions since many contradictions between observational studies and laboratory or mesocosm experiments exist in assessments of the toxicity of trace metals for aquatic insects [22].

Effluent from wastewater treatment plants (WWTPs) often represents an 'entrance gate' for many pollutants, including sewage-derived particulate organic matter (SDPOM), nutrients, pesticides, active pharmaceutical compounds (PhACs), synthetic organic compounds, phthalates, and trace metals into streams [14,23–27]. Over the past two decades, there has been a general decrease in organic matter and nutrient content in outputs from modernised WWTPs in the European Union [28]. Nonetheless, the continual release of compounds that are resistant to conventional purification technologies from WWTPs still occurs [29,30]. In addition, uncontrolled outflows of raw sewage from WWTPs and sewer infrastructure, especially during storm events, persist and represent a further challenge to the protection of natural water resources [31].

This article attempts to decipher the relationships between these human-induced gradients and EPT taxa and their feeding strategies which could be beneficial and permeable in such an impaired environment. It provides a list of taxa that can be considered tolerant to particular types of pollution, and the possible susceptibility of missing taxa is compared and discussed with previous reports. Since EPT taxa usually represent the essential part of macroinvertebrates in streams [3–5] and the shift in feeding strategies may reflect the change in available food resources, the article also provides partial evidence about the human-induced changes in the functioning of lotic ecosystems, namely, in nutrient cycling and energy flow [32].

#### **2. Materials and Methods**

#### *2.1. Locality Description*

The two surveyed localities—(i) Obecnický brook and (ii) Litavka river—are situated in central Bohemia (Czech Republic, Central Europe; spring to confluence: (i) N 49.7186939, E 13.8732253–N 49.7083778, E 13.9831206 and (ii) N 49.6563228, E 13.8557881–N 49.9602636, E 14.0847414; Figure 1). EPT larvae were sampled at four sites:


iron, lead, zinc, and uranium) that has been active for several centuries [15]. An active industrial complex that recycles, above all, lead waste occupies part of the river's alluvial plain. Heaps of waste material (especially sodium slag) are still stored several meters from the riverbank. Nevertheless, all mining activity has ceased. Here, the effects of agricultural disturbance are also likely to occur as arable fields cover approximately half of the deforested catchment area (excluding human settlements). The rest of the catchment is covered by pastures and meadows. The effect of the municipal WWs is presumably low because of a relatively small human population upstream. The riparian vegetation has been substantially reduced upstream along the Litavka river (third-order watercourse according to Strahler's system), which has been channelled in places between artificial banks. In addition, several shallow reservoirs (up to 7 ha), including some small sludge deposits, are situated upstream. We sampled the upper channelled parts with old reinforced banks and little riparian vegetation, as well as a lower unchannelled, naturally heterogenic stretch (morphologically similar to those described for Site 2) with banks lined with old alders and willows (Figure 1). The fish stocks here consist of common minnow (*Phoxinus phoxinus*), brown trout, and European perch (in descending order of abundance; Table S9). Insects prevailed in the macrozoobenthic samples, and small oligochaetes were relatively abundant, but crustaceans and molluscs were almost totally absent from the macrozoobenthic samples;

(iv) Site 4 was a stretch of the Litavka river downstream from its confluence with the Pˇríbramský brook (fourth-order watercourse according to Strahler's system; from N 49.7113508, E 14.0093011–N 49.7198211, E 14.0133511; Figures 1 and S2). This site was expected to be the most affected by anthropogenic activities because of the combination of industrial pollution from mining and smelting and contamination by treated municipal WWs from the town of Pˇríbram. Treated WWs are continually pumped into the Pˇríbramský brook approximately 900 m upstream from its confluence with the Litavka river. According to Grabicova et al. [14], Site 4 has the greatest concentrations of psychoactive PhACs of all the inspected localities, probably because of the low dilution possibilities (i.e., a small watercourse receiving relatively large amounts of municipal WWs). As at Site 3, we sampled an upper channelled stretch with artificial banks and little riparian vegetation, as well as a lower unchanneled section. Nonetheless, this morphologically heterogenic stretch only had a narrow and irregular fringe of riparian trees (alders and willows), possibly because of the dynamic migration of the river's course across its alluvial plain [15]. The fish stock consists of common minnow, common roach (*Rutilus rutilus*), chub (*Squalius cephalus*), gudgeon (*Gobio gobio*), brown trout, European perch, common bream (*Abramis brama*), common rudd (*Scardinius erythrophthalmus*), stone loach, and rainbow trout (*Oncorhynchus mykiss*) (in descending order of abundance; Table S9). The occurrence of limnophilic fish was due to the presence of aquaculture ponds upstream. There was no clear dominance by insects in the macrozoobenthic samples since small oligochaetes, leeches (e.g., *Erpobdella octoculata*, *Helobdella stagnalis*), and the water louse (*Asellus aquaticus*) were very abundant. Molluscs were detected in low numbers.

**Figure 1.** Map of the study area. The yellow circles indicate the sampling points. A set of three sampling points represents a whole plot. The contoured orange-coloured areas in the close-up view of Site 3 (lower left-hand corner) are slag heaps and waste material generated by the smelting industry. The river's flow direction is indicated by blue arrows.

#### *2.2. Sampling and Analysis of Samples*

The whole sampling campaign was conducted in 2020, specifically on 11 May, 24 June, 11 August, and 24 September. The macrozoobenthos were sampled using a Surber sampler (30 × 30-cm frame, 500-μm mesh) at four sampling sites and three points (triplicates) at each site (see part 2.1. 'Locality description' and Figure 1). Each triplicate consisted of three subplots within a whole plot or site. The location of the whole plot was chosen at random during each sampling campaign. A diversification of samples within each whole plot was made by sampling riffles (shallow and with high water velocity), inflows into pools (deep and with high water velocity) and the backs of pools (deep and with low water velocity). Samples were processed using a round steel sieve (40 cm diameter, 500 μm mesh) and preserved in 70% technical ethanol. Organisms were sorted in the laboratory from the material taken from the stream bottom. Mayfly, stonefly, and caddisfly larvae were determined to the deepest possible taxonomic level using a binocular microscope and stereomicroscope (Olympus SZX16) and determination keys; their frequencies (abundances) in the samples were recorded. The methodology used to analyse the water samples for pesticides, PhACs and drugs, and the water and sediment samples for trace metals is described in detail in the Supplementary Materials (Section Materials and Methods).

Environmental data were taken for both the whole and subplot levels. The environmental variables (EVs) differing only at the whole-plot level (hereafter referred to as 'whole-plot EVs') consisted of (i) physicochemical parameters—conductivity, pH, oxygen concentrations, and water temperature—measured using a HACH® multi-meter, (ii) pesticide and PhAC concentrations in water (measured from samples taken on 11 August 2020), (iii) cadmium, lead, and zinc concentrations measured in sediments and water (measured from samples taken on 24 September 2020), and (iv) daily volumes of sewage (untreated or 'poorly treated' municipal WWs) recorded up to 15 days prior to each sampling day at the Pˇríbram WWTP (data provided by 1. SˇcV JSC, Pˇríbram, Czech Republic). Detailed monthly data of water and sediment analyses from Sites 1, 2, and 3 (2009–2020) and data of actual water discharge from the same sites were provided by the Czech Hydrometeorological Institute (Prague, Czech Republic) and Povodí Vltavy, State Enterprise (Prague, Czech Republic). In addition, the water temperature was continuously monitored during the sampling campaign by four data loggers (TFA) placed approximately in the middle of each site.

The EVs also differed at the subplot level (hereafter referred to as 'subplot EVs') in terms of their (i) choriotopic compositions, which were empirically estimated for each Surber sampling spot; (ii) current velocities [m·s–1] (minimum, maximum and average), measured continuously for 0.5 min above the Surber sampling spot with a flowmeter (MiniController MC20 with the Flowprobe for MiniWater20, Schiltknecht Messtechnik AG, Gossau, Switzerland) placed close to the bottom, in the water column and close to the water surface; (iii) widths and depths at each subplot; and (iv) classification as riffles, pools, inflows into pools, or backs of pools.

#### *2.3. Statistical Analyses*

Data were analysed using Canoco 5 version 5.12 (written by Cajo J.F. ter Braak and Petr Šmilauer; Wageningen, the Netherlands, and P. Šmilauer, Czech Republic) [35] and R-studio version 4.1.1 (written by RStudio team; Boston, MA, United States) [36] software. The analysis of shifts in feeding strategies was based on the preferences of individual species according to Graf et al. [37], Buffagni et al. [38], and Graf et al. [39]. The analyses with species community composition as response variables were conducted using unimodal Canonical correspondence analysis (CCA), while analyses with feeding strategies as response variables were performed using linear Redundancy analysis (RDA). Species abundances were log + 1 transformed before the CCA. Detrending was not used because no dependency between positions of samples on the first or second ordination axes was observed. Feeding strategies were averaged by species composition-weighted standardisation (hereafter referred to as 'averaged feeding strategies'). Only centring by species was applied using RDA.

The significances of axes constrained by (i) site identity, (ii) the identity of sampling times, (iii) whole-plot EVs, and (iv) subplot EVs (see part 2.2. 'Sampling and sample analysis') were tested using Monte-Carlo permutation tests (999 permutations used). Permutations were carried out according to the experiment's hierarchical design: (i) when the whole-plot EVs or the identity of sites were used as explanatory variables, both whole-plots and subplots in the four blocks defined by the four sampling times used as covariates were permuted using cyclic shifts; (ii) when the subplot EVs were used as explanatory variables, only the subplots in the four blocks defined by sampling times were permuted using cyclic shifts; and (iii) when the identity of sampling times was used as an explanatory variable, both whole plots and subplots in the four blocks defined by the four sites were permuted using cyclic shifts. A false discovery rate was used for adjusting *p* values. Forward selection (FS) was used for selecting significant whole- and subplot EVs. 'Hybrid analysis' was chosen when tests on visualised constrained axes were not significant (*p* > 0.05). The significance of differences in choriotopic compositions between sites was tested in a similar way; individual values were transformed by the angular (inverse sine) function and the transformed values were used as response variables. This transformation was also applied when the choriotopic compositions were used as explanatory variables.

Generalised linear models (GLMs) were employed to model trends in shifts of multiple averaged feeding strategies along environmental gradients or along the resulting ordination axes. The significance of differences between sites and sampling times in EPT abundances and richness was tested using GLMs with negative binomial and Poisson distributions, respectively (the quasi-likelihood estimation method was applied in the case of the Poisson distributions). The significance of improvements using Generalised linear mixed-effects models (GLMMs) rather than GLMs was tested by likelihood ratio tests. A post hoc test with Tukey's method for *p* value adjustment was applied.

#### **3. Results**

#### *3.1. Environmental Conditions*

The differences between sites are summarised in Table 1. Several environmental gradients, including increasing temperature, pH, and conductivity and decreasing oxygen concentrations, changed along with the longitudinal downstream profile (from Site 1 to Site 4). For more information about the basic physicochemical parameters at other sites in the longitudinal profile, see Supplementary Materials (Section Materials and Methods; Table S1). The greatest total concentration of cadmium, zinc, and lead in the waters was detected at Site 3, followed by Site 4, while the highest loads of pesticides and PhACs were found at Site 4 below the discharge of effluent from the WWTP (Table 1). Treated WWs constituted 1/13–1/3 of all discharges at Site 4 (Table S2). The periodical releasing of untreated or 'poorly treated' WWs containing organic material and coarse communal waste (mainly hygienic products) also occurred. The monthly released volumes of this sewage ranged from 299 to 130,016 m3 in 2020 (Table 2); however, considerable volumes also leaked from sewers upstream from the WWTP but were not included in the analyses. For more information about contamination at other sites along the longitudinal profile, see Supplementary Materials (Section Materials and Methods; Tables S3–S6).

The average current velocities at individual sites were as follows: Site 4 >Site 3 > Site 1 > Site 2 (in descending order from the fastest to the slowest; Table 3). Differences in the choriotopic compositions at sites tested using a RDA were marginal (test on first axis: pseudo-*F* = 1.2, *p* = 0.085; test on all axes: pseudo-*F* = 2.8, *p* = 0.003; Figure S4); only Site 4 and Site 2 were significantly different from the others when tested with the RDA (explained variation = 7.5%, pseudo-*F* = 3.8, *p* (adj.) = 0.025 and explained variation = 6.0%, pseudo-*F* = 2.9, *p* (adj.) = 0.049, respectively). Samples from Site 1 compared with those from Site 2 had a much greater proportion of smaller mineral particles (size < 6.3 cm— 'microlithal', 'akal', and 'psammal') and 'mosses' (Marchantiophyta only at the Site 1),

and a lower proportion of light deposits ('coarse particulate organic matter'—CPOM and 'fine particulate organic matter'—FPOM), debris, and xylal (Table 4). Samples from Site 3 resembled most those from Site 1, with similar proportions of smaller mineral particles, deposits, debris, and xylal (Table 4). Finally, Site 4 differed most from the others because of the highest proportions of deposits and debris and the occurrence of filamentous algae and anthropogenic trash (especially wet wipes; Table 4).

**Table 1.** Resulted mean physicochemical parameters measured at the given sites and concentrations of metals, pesticides, and active pharmaceutical compounds (PhACs) in water samples taken at these sites.


**Table 2.** Monthly volumes of untreated or 'poorly treated' wastewater released from the wastewater treatment plant into Site 4 in 2019 and 2020.


**Table 3.** Resulted mean current velocities measured for the macrozoobenthic samples from each site.



**Table 4.** Resulted mean choriotopic percentual compositions empirically estimated for each macrozoobenthic sample at all given sites.

#### *3.2. EPT Abundance and Richness*

In total, 87 EPT taxa were detected in the samples: 11 mayflies, 29 stoneflies, and 47 caddisflies (Table S7). EPT abundances were not substantially different between sites (Figure 2A; Table 5) or sampling times (Table 5). Nevertheless, total EPT richness estimated for individual sites did show a gradual decreasing and significant trend from Site 1 to Site 4 (Figure 2B; Table 5). Sampling time had a nonsignificant effect on EPT richness (Table 5).

The greatest differences between sites in terms of abundance or richness were observed in stoneflies followed by mayflies (Table 5), both of which were impaired at disturbed Sites 3 and 4 (Figure 2C–E). Conversely, caddisflies were more abundant at disturbed Sites 3 and 4 (Figure 2C), although caddisfly richness and family richness did not significantly differ between sites (Figure 2D,E; Table 5). Significant differences between the

four sampling times in abundance or richness were detected in mayflies and stoneflies but not in caddisflies (Table 5). The abundance of mayfly nymphs was lowest in August and highest in May, whilst the opposite was found in stoneflies: the highest stonefly abundance was detected in August and the lowest in May.

**Figure 2.** Bar graphs visualising (**A**,**C**) 'abundance' (number of individuals), (**B**,**D**) 'richness' (number of taxa), (**E**) 'family richness' (number of families), and (**F**–**H**) abundances and richness of individual detected families of Ephemeroptera, Plecoptera, and Trichoptera, respectively, summarised for each 'site' from the samples taken at all sampling times. Values in the bar graphs (**F**–**H**) are transformed by a natural logarithm, and the 'richness' is displayed above each column. E—Ephemeroptera, P—Plecoptera, T—Trichoptera.

The stonefly and caddisfly orders significantly differed in abundance between sites (Table 5). Stonefly abundances rapidly declined at disturbed Sites 3 and 4 compared with Sites 1 and 2, whilst caddisfly abundances were lower at Sites 1 and 2 (Figure 2C). Mayfly abundance did not significantly differ between sites (Figure 2C, Table 5) but marginally differed between sampling times (Table 5). No significant interaction was detected between the four-level factor site and sampling time.



The total abundances and species richness shown in Figure 2F,G indicate that practically only members of the Baetidae family (namely *Baetis fuscatus*, *B. rhodani*, *B.* aff. *scambus*, and *B. vernus*) and the Leuctridae family (stoneflies) were detected at disturbed Sites 3 and 4. At these sites, *Habrophlebia lauta* (Leptophlebidae) and *Seratella ignita* (Ephemerellidae) all but completely disappeared, while *Paraleptophlebia submarginata* (Leptophlebiidae), *Ecdyonurus torrentis*, and members of the genus *Rhithrogena* (both Heptageniidae) were totally absent, although these latter taxa were detected upstream at the interconnected Sites 1 and 3 (Figure 2F, Table S7). Of the stoneflies, only *Leuctra albida*, *L. fusca*, and *L. geniculata* were detected at Site 3 (*L. geniculata* was only detected at this site, Table S7) and only *L. fusca* in low abundances and body sizes was detected at Site 4 (Figure 2G, Table S7). *Siphonoperla torrentium* (Chloroperlidae), several species belonging to the genera *Protonemura*, *Amphinemura*, and *Nemoura* (Nemouridae), all very abundant at Sites 1 and 2, along with *Diura bicaudata* and *Perlodes* aff. *microcephalus* (Perlodidae) were completely absent in the samples taken from disturbed Sites 3 and 4 (Figure 2G, Table S7). Nevertheless, compared with Site 1, the richness of Leuctridae decreased also at Site 2 below the water reservoir (Figure 2G), and, for instance, *Leuctra nigra*, frequent in samples from Site 1, was missing from samples from Site 2 (Table S7).

Of the caddisflies, Site 1 had a greater richness of trichopterans than Site 2; the abundant family Philopotamidae, for example, was no more detected at Sites 2, 3, and 4; finally, members of the Goeridae were detected only sporadically at particular sites (Figure 2H). On the other hand, the richness of Hydropsychidae was lowest at Site 1. *Odontocerum albicorne* (Odontoceridae) and *Sericostoma flavicorne*/*personatum* (Sericostomatidae), very abundant at both Sites 1 and 2, were completely missing from the samples from Sites 3 and 4 (Figure 2H, Table S7). One of the dominant families at Site 3 was Limnephilidae, represented above all by species from the genera *Potamophylax*, *Halesus*, and *Chaetopteryx* (in descending order of total abundance; Figure 2H, Table S7). One of the dominant families at Site 4 below the WWTP was Hydropsychidae, represented by *Hydropsyche siltalai*, *H. angustipennis*, and *H. incognita*/*pellucidula* (in descending order of total abundance; Figure 2H, Table S7); the population of limnephilids at Site 4 was substantially lower (Figure 2H). At both disturbed Sites 3 and 4, the family Polycentropodidae was much more abundant than at Sites 1 and 2 (Figure 2H). A substantially greater proportion in the total abundances of *Cyrnus trimaculatus* than *Polycentropus flavomaculatus* was observed at Site 4 than at Site 3 (2.54 and 0.03, respectively; Table S7). The Leptoceridae represented by *Athripsodes bilineatus* was most abundant at Site 3, while the Psychomyiidae, represented only by *Psychomyia pusilla*, were detected only at Site 4 (Figure 2H; Table S7). The total abundances of the family Rhyacophilidae were almost identical at all sites (Figure 2H). Even though the clear morphological determination of larvae from the group *Rhyacophila* sensu stricto is not always possible [40], *Rhyacophila* cf. *aurata* was dominant in samples from Sites 3 and 4 (especially Site 3), whilst *Rhyacophila nubila* gr. was commonest in samples from Sites 1 and 2.

#### *3.3. Shift in EPT Community Composition along Environmental Gradients*

The differences in species compositions between sites were significant (test on first axis: pseudo-*F* = 1.9, *p* = 0.001; test on all axes: pseudo-*F* = 4.0, *p* = 0.001; Figure S5), as were the differences in composition of species between sampling times (but only the composition of specimens > 0.5 mm in size; test on first axis: pseudo-*F* = 0.7, *p* = 0.024; test on all axes: pseudo-*F* = 1.8, *p* = 0.001; Figure S6). The whole- and subplot EVs significantly correlated with the shift in species composition (*p* < 0.05); those selected by FS are shown in Figures 3 and 4, respectively.

**Figure 3.** Species + whole-plot environmental variables (EVs) biplots: the first 33 best fitting species are labelled (Ephemeroptera, Plecoptera, and Trichoptera labelled in red, green, and blue, respectively). A partial canonical correspondence analysis (CCA) was used: (**A**) first and second ordination axes; (**B**) second and third ordination axes. Ordination axes were constrained by the EVs selected by forward selection, which accounted for 30.48% of the explained variation. Test on first axis: pseudo-*F* = 0.5, *p* = 0.001; test on all axes: pseudo-*F* = 2.3, *p* = 0.001. Conditional effects (from the greatest to the lowest explained variation): 'Log-transformed sum of cadmium, lead and zinc (Met)': explained variation = 11.7%, pseudo-*F* = 5.7, *p* (adj.) = 0.003; 'Conductivity (Cond)': explained variation = 6.1%, pseudo-*F* = 3.1, *p* (adj.) = 0.020; 'Oxygen concentration (Oxy)': explained variation = 4.0%, pseudo-*F* = 2.1, *p* (adj.) = 0.012; 'Log-transformed sum pesticides (Pest)': explained variation = 3.0%, pseudo-*F* = 1.6, *p* (adj.) = 0.045; 'pH' explained variation = 2.9%, pseudo-*F* = 1.6, *p* (adj.) = 0.047; 'Total volume of released sewage for three days prior to sampling time (Sew)': explained variation = 2.9%, pseudo-*F* = 1.6, *p* (adj.) = 0.036. Abbreviations: *AmpBor*—*Amphinemura borealis*, *AmpSul*—*A. sulcicollis*, *Ath-Bil*—*Athripsodes bilineatus*, *BaeAffSc*—*Baetis* aff. *scambus*, *BaeFus*—*B. fuscatus*, *ChaeCfVi*—*Chaetopteryx* cf. *villosa*, *CyrTri*—*Cyrnus trimaculatus*, *DruAnn*—*Drusus annulatus*, *HabLau*—*Habrophlebia lauta*, *HaleCfTe*—*Halesus* cf. *tesselatus*, *HydInf*—*Hydatophylax infumatus*, *HydAng*—*Hydropsyche angustipennis*, *HydInc*—*H. incognita*, *HydInc/Pel*—*H. incognita*/*pellucidula*, *HydPel*—*H. pellucidula*, *HydSax*— *H. saxonica*, *HydSil*—*H. siltalai*, *LeuAlb*—*Leuctra albida*, *LeucCfAr*—*L.* cf. *armata*, *LeuBra*—*L. braueri*,

*LeucCfPs*—*L.* cf. *pseudocingulata*, *LeucSp*—*Leuctra* sp., *LeuFus*—*L. fusca*, *LeuGen*—*L. geniculata*, *LeuIne*— *L. inermis*, *LeuNig—L. nigra*, *NemoCfUn*—*Nemoura* cf. *uncinata*, *NemoSp*—*Nemoura* sp., *NeuBim*— *Neureclipsis bimaculata*, *OdoAlb—Odontocerum albicorne*, *PhiLud*—*Philopotamus ludificatus*, *PolFla*— *Polycentropus flavomaculatus*, *PotaCfCi*—*Potamophylax* cf. *cingulatus*, *PotAffRot*—*P.* aff. *rotundipennis*, *PotRot*—*P. rotundipennis*, *ProHra*—*Protonemura hrabei*, *ProMey*—*P. meyeri/nitida*, *ProNit*—*P. nitida*, *ProtCfLa*—*P.* cf. *lateralis*, *PsyPus*—*Psychomyia pusilla*, *RhiSem*—*Rhithrogena semicolorata*, *RhitCfPu*—*R.* cf. *puytoraci*, *RhyaCfAu*—*Rhyacophila* cf. *aurata*, *RhyAffDo*—*R.* aff. *dorsalis*, *RhyAffNu*—*R.* aff. *nubila*, S*erFla/Per*—*Sericostoma flavicorne*/*personatum*, *SerIgn*—*Serratella ignita*.

**Figure 4.** Species + subplot environmental variables (EVs) biplot: the first 24 best-fitting species are shown (Ephemeroptera, Plecoptera, and Trichoptera labelled in red, green, and blue, respectively). A partial Canonical correspondence analysis (CCA) was used. Ordination axes were constrained by the EVs selected by forward selection, which accounted for 10.54% of the explained variation. Test on first axis: pseudo-*F* = 0.1, *p* = 0.024; test on all axes: pseudo-*F* = 1.4, *p* = 0.001. Conditional effects (from the greatest to the lowest explained variation): 'Arcsine-transformed proportion of mosses in choriotope (Mosses)': explained variation = 5.3%, pseudo-*F* = 2.4, *p* (adj.) = 0.015 and 'Mean current velocity (Velo)': explained variation = 5.2%, pseudo-*F* = 2.5, *p* (adj.) = 0.023. Abbreviations: *AmpSul*—*Amphinemura sulcicollis*, *BaeFus*—*B. fuscatus*, *HabLau*—*Habrophlebia lauta*, *HaleCfTe*—*Halesus* cf. *tesselatus*, *HydInc/Pel*—*Hydropsyche incognita*/*pellucidula*, *HydSax—H. saxonica*, *HydSil*—*H. siltalai*, *LeuBra*—*L. braueri*, *LeucCfAr*—*L.* cf. *armata*, *LeuNig*—*L. nigra*, *LitNig*—*Lithax niger*, *NeuBim*—*Neureclipsis bimaculata*, *OdoAlb—Odontocerum albicorne*, *PhiLud*—*Philopotamus ludificatus*, *PleCon*—*Plectrocnemia conspersa*, *PolFla*—*Polycentropus flavomaculatus*, *ProHra*—*Protonemura hrabei*, *RhitCfPu*—*Rhithrogena* cf. *puytoraci*, *RhyaCfAu*—*Rhyacophila* cf. *aurata*, *RhyAffNu*—*R.* aff. *nubila*, *RhyaSp*—*Rhyacophila* sp. (early instar larvae), *RhyNubGr*—*R. nubila* gr., S*erFla/Per*—*Sericostoma flavicorne*/*personatum*.

#### *3.4. Shift in Averaged Feeding Strategies along Environmental Gradients*

The differences in compositions of the averaged feeding strategies between sites were significant in general; however, only the composition at Site 4 was significantly different from all other sites (*p* < 0.05; Figure 5A). There was also a significant shift in the whole-plot EVs (*p* < 0.05). Significant relationships with 'conductivity' (32.5% of explained variability), 'total volume of released sewage for three days prior to sampling time' (5.9% of explained variability), 'oxygen concentration' (2.8% of explained variability), and 'total concentration of PhACs' (2.3% of explained variability) were detected (*p* < 0.05). 'Passive filter feeder' and 'predator' feeding strategies were dominant in the WW-polluted environment; the 'grazer and scraper' feeding strategy seemed to be limited by the most polluted WW; and the 'gatherer/collector' and 'shredders' feeding strategies were generally suppressed in the most polluted environment (Figures 5B and S7B,C). The shift in the subplot EVs was not significant (test on first axis: pseudo-*F* = 0.6, *p* = 0.379; test on all axes: pseudo-*F* = 1.9, *p* = 0.168).

**Figure 5.** (**A**) Averaged feeding strategies pie chart and (**B**) Averaged feeding strategies + whole-plot environmental variables (EVs) biplot. Partial Redundancy analyses (RDA) were used. First ordination axes were constrained: (**A**) by the factor 'Site' that accounted for 30.98% of the explained variation and (**B**) by the EVs selected by forward selection, which accounted for 38.91% of the explained variation. (**A**) Test on first axis: pseudo-*F* = 5.4, *p* = 0.003; test on all axes: pseudo-*F* = 6.1, *p* = 0.026; (**B**) test on first axis: pseudo-*F* = 2.5, *p* = 0.006; test on all axes: pseudo-*F* = 5.7, *p* = 0.009. Conditional effects (from the greatest to the lowest explained variation): (**A**) 'Site 4 (S4)' explained variation = 22.0%, pseudo-*F* = 12.1, *p* (adj.) = 0.027; 'Site 3 (S3)' explained variation = 8.6%, pseudo-*F* = 5.2, *p* (adj.) = 0.090; and 'Site 2 (S2)' explained variation = 0.3%, pseudo-*F* = 0.2, *p* (adj.) = 0.943. (**B**) 'Conductivity (Cond)': explained variation = 32.5%, pseudo-*F* = 20.7, *p* (adj.) = 0.005 and 'Total volume of released sewage for three days prior to sampling time (Sew)': explained variation = 6.4%, pseudo-*F* = 4.4, *p* (adj.) = 0.045.

#### *3.5. Malformations and Mortality of Caddisflies*

A substantial proportion of final instar caddisfly larvae in *Halesus* cf. *tesselatus*, *Potampohylax* aff. *rotundipennis*, and *R.* cf. *aurata* sampled at Site 3 had malformed limbs (11.8%, 7.0%, and 7.7%, respectively). Peculiarly fused terminal parts of walking legs (between tibiae and tarsal claws) or overall shortened walking legs were observed in the two given limnephilid taxa; deformed and shortened anal claws were observed in *R.* cf. *aurata* (Figure 6). Malformed specimens of *H.* cf. *tesselatus* were also observed at Site 4 (28.6%), although at the latter site, the total number of sampled specimens was substantially lower than at Site 3 (7 vs. 76 specimens at Site 3). There were no malformed specimens in the

abundant *R.* cf. *aurata* in Site 4; however, there was one malformed specimen out of a total of 20 specimens of *R. nubila* gr. at Site 4.

**Figure 6.** Malformed limbs of *Halesus* cf. *tesselatus* (two pictures on the left) and *Rhyacophila* cf. *aurata* (right side) sampled at Sites 3 and 4.

We also observed a drift of dead limnephilid pupae at Site 3 in September during the macrozoobenthic sampling. We found dead and decaying larvae and pupae of *Drusus annulatus*, *H.* cf. *tesselatus*, and Limnephilidae gen. sp. inc. in the given samples (14 dead limnephilid specimens and 7 living limnephilid specimens in total). Several dead and decaying pupae of *Hydropsyche* sp. were also found in the samples from Site 2 taken in June.

#### **4. Discussion**

#### *4.1. Variability in the EPT Community along an Environmental Gradient*

The EPT community can significantly reflect the quality of the environment (Figures 3–5). Even though our study only consisted of a small-scale natural experiment (i.e., we only sampled two mutually interlinked watercourses), the results match, in general, those of a study of whole macroinvertebrate communities carried out in 23 Swiss streams [41]. As in our work, in the Swiss study, about 30% of community variability was explained by water quality; conductivity (which correlates to total dissolved solids) together with oxygen concentrations explained 10.1% in our study, while concentrations of pesticides explained 3.0% of community variability in both studies (Figure 3).

The changes in EPT abundance and richness related to pollution by trace metals observed in our study generally resemble reports from other sites, the greatest similarity being in the susceptibility of mayflies—and especially heptageniids—to this type of pollution (Figure 2). Sites 3 and 4 were both warmer—maximums up to 23 ◦C—than unpolluted Sites 1 and 2 (Table 1; Figure S3), a relevant finding given that the temperature gradient, which is often related to the toxicity of chemical compounds and the solubility of oxygen in water, is considered a good predictor of change in the composition of EPT communities [42]. The increasing temperature gradient observed when moving downstream is undoubtedly natural; however, the difference in the temperature regime between the contaminated and uncontaminated sites is also likely to be the product of the canalisation, the removal of the riparian vegetation, and the presence of shallow reservoirs upstream on the Litavka

river (see part 2.1 'Locality description' [43]). That many species of mayflies and stoneflies would avoid high temperatures, mine effluent and certain trace metals was predicted beforehand [17–19,44–46]. Our results show that baetids were the most resistant (or resilient) to the trace metal contamination, and within this group, the dominant taxa were the *Baetis fuscatus* group (i.e., *B. fuscatus* and *B.* aff. *scambus,* whose separation was not possible since the taxon identified and presented here as *B.* aff. *scambus* shared features of both species). The total abundance of *B. fuscatus* gr. compensated for the decrease in abundances of other mayfly taxa so that this small-sized baetid could act as a proxy for other mayflies in human-disturbed sites. The other most commonly detected baetids, *B. rhodani* and *B. vernus*, did not show any substantial changes in abundance between sites. Nevertheless, trace metal pollution can have an effect later on in the life cycle of baetids and other aquatic insects by rapidly causing increased mortality in emerging imagines (or subimagines) [47,48].

Stoneflies are generally the aquatic insects most sensitive to conditions modified by humans [49,50]. There was significantly decreased richness in stoneflies—even at Site 2 compared with Site 1—and their abundance only significantly fell at Sites 3 and 4 (Figure 2; Table 5). We interpret their decrease in richness at Site 2, situated below the reservoir, as the result of the discontinuity in the river network, i.e., a change in the deposition regime (a lower proportion of smaller mineral particles and a greater proportion of light deposits; Table 4), and the separation of Site 2 from the network of tributaries in the crenal zone situated upstream from the reservoir (Figure 1). However, there are several indications that stoneflies are less sensitive to metals than mayflies. Many stoneflies (as well as caddisflies) are more tolerant to increased acidity than many mayflies [51]. Lower pH leads to natural access to metals (e.g., aluminium) originally bound to substrates [52,53]. Observational work studying the effects of contamination from mining has reported relatively high resistance in some stoneflies, including *Amphinemura sulcicollis*, *Leuctra fusca*, and *L. inermis*, taxa that are similar to those found in our study [16,19]. Little evidence exists regarding the levels of trace metals concentrations that are lethal for stoneflies [54], although certain North American stoneflies have been found to withstand (in significantly reduced numbers) concentrations of 11, 120, and 1100 <sup>μ</sup>g·L–1 of cadmium, copper, and zinc, respectively, for 10 days in microcosm conditions. By comparison, numbers of heptageniids and *Baetis tricaudatus* fell significantly in only half of these concentrations [46].

Therefore, the significant reduction in stonefly abundance and richness at both contaminated sites cannot be conclusively attributed to pollution by trace metals. The absence there of most of the stoneflies detected upstream could also be due to water warming since most stoneflies prefer lower temperatures than those recorded at Sites 3 and 4 [39,42]. However, the interaction of multiple effects, including habitat degradation and several types of pollution, are assumed to take place.

According to our data, the less sensitive stoneflies were able to withstand the heavy industrial trace metal pollution but were then eliminated by the municipal WW, despite the advanced technology used in the WW treatment [14]. A synergic effect of multiple stressors is likely to occur. We can demonstrate that specimens of all sizes of *Leuctra fusca* gr. (species identified as *L. albida* and *L. fusca* in this study [55]) and *L. geniculata* occurred in the water contaminated by trace metals at Site 3 (Table S7) but not in the water contaminated by municipal WWs. Only *L. fusca* gr. (the only species identified as *L. fusca*) of small sizes were detected at Site 4, contaminated by both trace metals and municipal WWs. It is important to note that the organic pollution produced by the WWs sometimes reaches high levels. According to unpublished data provided by the Czech Hydrometeorological Institute and Povodí Vltavy, State Enterprise, the average estimated five-day biochemical oxygen demand (BOD5) was equal to 2.39 mg·L–1 at Site 4; however, the estimated maximum was 17.00 mg·L–1 just four months before the beginning of our sampling campaign. In the Pˇríbramský brook (Auxiliary Site 4 downstream from the effluent from the WWTP, see Supplementary Materials), the estimated BOD5 at the same time was

24.00 mg·L–1 (the average and maximum of the estimated BOD5 were 4.76 and 60.00 mg·L–1 in 2013–2020, respectively).

Caddisflies were the least harmed insects in the EPT community. Their richness was not significantly different between sites, and their abundance was significantly higher at Site 4 than at Sites 1 and 2 (Figure 2; Table 5). Many caddisflies are acid-tolerant or acid environment specialists (e.g., *Chaetopterygopsis maclachlani* and *Rhyacophila polonica* detected during our study [51,56]), and some are reported to be tolerant to trace metal pollution (e.g., *Hydropsyche* sp., leptocerids, *Rhyacophila* sp. [17,18,46]). Ubiquitous caddisfly taxa can occupy warmer sites polluted by organic materials with low oxygen concentrations (e.g., potamal hydropsychids and *Psychomyia pusilla* [56,57]). Nevertheless, caddisflies are generally very sensitive to pesticides and avoid sedimentation and droughts, as does the EPT group in general [11,58,59].

The vast majority of EPT taxa situated on the right-hand side of the CCA biplot (Figure 3A) were caddisflies. The relative abundances of these caddisflies were positively correlated to the human-induced gradients. The caddisfly communities at Sites 1 and 2 were quite similar; however, net-spinning philopotamids were lacking from Site 2 and were replaced by net-spinning hydropsychids and polycentropodids (Figure 2H). Changing conditions along the longitudinal profile (e.g., food source and current velocity) probably played an important role in shaping the net-spinning caddisfly community. Philopotamids prefer diatoms and fine detritus and have been reported to be quite sensitive to organic pollution [60,61]). According to our data, philopotamids were found at the subplots with the greatest maximum water velocity (0.73–0.97 m·s–1). Slower water speeds caused by water abstraction [60] combined with the natural longitudinal stream gradient [37,60] favoured hydropsychids and polycentropodids. Unlike philopotamids, hydropsychids and polycentropodids are either facultatively or obligatory carnivorous [37,62].

The altered conditions at Sites 3 and 4 probably prevent the occurrence of the eruciform caddisflies *O. albicorne* and *S. flavicorne*/*personatum*, which were dominant in the samples from Sites 1 and 2. The warmer temperatures (and poorer oxygen conditions) are expected to have a negative effect. Haidekker and Hering [42] have reported that *O. albicorne* tolerates well mean summer temperatures up to 16 ◦C. However, this figure was exceeded at Site 4 by 1.24 ◦C (Figure S2), and so, based on individual measurements (Table 1), we assume that a similar temperature regime exists at Site 3. Furthermore, both these taxa are partially endobenthic and are more exposed to the toxic compounds that accumulate in sedimented particles [63]; the long larval development time in *S. flavicorne*/*personatum* [64] will also prolong the exposure of its larvae to toxicants. On the other hand, high abundances of large eruciform caddisflies—limnephilids of the genus *Halesus* and *Potamophylax*—were found at Site 3. These taxa are probably tolerant to trace metal pollution. Nonetheless, we found malformed walking legs in up to 11.8% of individuals that, in light of previous studies, could be a sign of trace metal contamination [65] or other persistent and hazardous compounds generated by the industrial complex upstream from Site 3 [66].

The lower abundances and richness of limnephilids observed at Site 4 downstream from the WWTP effluent could be a consequence of the interaction between trace metals and WW pollution (or the WW pollution alone). Hydropsychids, polycentropodids, rhyacophilids, and psychomyiids (*P. pusilla*) were most abundant downstream from the WW effluent. The taxa are shown in the upper right-hand corner of Figure 3A were positively related to the 'total volume of released sewage three days prior to the sampling' (hereafter referred to as '3-day sewage volumes'), which probably contained a high proportion of SDPOM—a mixture of organic detritus and microorganisms such as bacteria and algae that is a high-quality food resource downstream from WW effluents [23]. On the other hand, the released sewage could trigger a drift of benthic prey organisms with greater oxygen demands (e.g., mayfly nymphs), and SDPOM could be used by periphyton and macrophytes (aquatic and riparian), which can serve as microhabitats or food. For example, the carnivorous *C. trimaculatus* with a more pronounced microhabitat preference for macrophytes [37] was more common than *P. flavomaculatus* at Site 4 than Site 3. The most

abundant hydropsychids can use multiple food resources and, aside from predating on small benthic organisms and consuming detritus, they also graze on algae [37,67]. Besides crustaceans, insects (including EPT), and terrestrial prey, detritus also constitute a substantial portion of gut content in polycentropodids [68]. It cannot be rejected that hydropsychids and polycentropodids could have a detrimental impact on the populations of other EPT taxa. Nevertheless, their net spinning increases heterogeneity of current velocities and nets represent valuable microhabitats, e.g., positive relationships to ephemerellid mayfly and chironomids were reported [69,70]. We associate the high abundance of early instars of *H. siltalai* in September at Sites 3 and 4 with the presence of aquatic mosses (Figure 4).

The proportion of mosses in the choriotope (explained community variation = 5.3%) and mean current velocity (explained variation community variation = 5.2%) had the highest parsimonious and significant power to explain EPT community composition at the subplot level (Figure 4). These two variables were positively correlated because mosses were detected mainly in shallow, fast-flowing habitats. Mosses (as well as roots) provide better vertical zonation possibilities for aquatic organisms, as well as food and shelter [71]. The information provided by Buffagni et al. [38], Graf et al. [37], and Graf et al. [39] does not correspond particularly well to our results: of the species best corresponding to the presence of aquatic mosses only *Hydropsyche saxonica, Plectrocnemia conspersa*, and *Protonemura hrabei* were reported to have a certain preference for macrophytes (mosses included). Nevertheless, Krno [72] suggests that the goerid *Lithax niger* is a bryophyte dweller and reports that *Rhyacophila obliterata* can only live on emergent mosses, whilst *R. nubila*, *R. polonica* can also live on submerged mosses. Similar findings have been reported by Glime [71]. Conversely, Graf et al. [37] consider all goerids and rhyacophilids to inhabit micro-/mesolithal and macro-/megalithal and suggest that they are 'specialised for lithal'.

A high mean velocity best corresponded to the highest relative abundance of *H. pellucidula*/*incognita* (the larvae of these two species are often not separable [40]): fitted optimum over 0.60 m·s–1; *R.* cf. *aurata* and *R. nubila* gr.: fitted optimum around 0.50 m·s–1; and *H. siltalai*: fitted optimum 0.47 m·s–1 (Figure 4). The mayfly *B. fuscatus* had a fitted optimum around 0.35 m·s–1 (Figure 4). By contrast, according to our data, *O. albicorne* had a fitted optimal mean current velocity around 0.14 m·s–1, *S. flavicorne*/*personatum* 0.10 m·s–1, and the mayfly *H. lauta* less than 0.10 m·s–1.

#### *4.2. Shift in Averaged Feeding Strategies in the EPT Community along Environmental Gradients*

The analysis of shifts in ecological traits and, above all, their relationship to environmental characteristics has recently become a key topic in community ecology [73]. Substantial changes in the composition of feeding strategies (functional feeding groups) due to disturbances are assumed to occur in the macrozoobenthos [74,75]. Significant differences in the composition of average EPT feeding strategies were detected between localities (Figure 5A). Conductivity and 3-day sewage volumes were found to be the most parsimonious significant predictors in our dataset (Figure 5B). Both these EVs were positively related to the input of various materials and were positively related to both passive filter feeder and predator feeding strategies (represented almost only by caddisflies; Figure S7B). It suggests in higher susceptibility of EPT shredders to the pollution or unavailability of the CPOM (main food resource utilised by shredders) edible for these organisms (e.g., crystalline deposits containing trace metals can precipitate on the external surface of alder leaves because of fungal grown [21]). Even though shredder abundance will be compensated by numerous *Asellus aquaticus* and gatherers/collectors by chironomids and oligochaetes (see part 2.1. 'Locality description'), the abundance of winged adults of the aquatic insect directly utilising CPOM may be substantially lowered and could cause an alteration in nutrient and energy flows through the environment.

Despite that, 3-day sewage volumes positively correlated to the responses of passive filter feeders, predators, and shredders and negatively to the responses of gatherers/collectors (Figure S8A), the greatest 15-day sewage volumes (Figure S8B) positively correlated to the responses of gatherers/collectors and grazers and scrapers. Nevertheless, this trend

was not strongly associated with the '15-day sewage volumes' because the dominance of gatherers/collectors and grazers and scrapers was caused by a high abundance of nymphs of the mayfly *B. fuscatus* gr., which probably colonised this site after drifting from an upstream stretch after a high discharge event (up to 8.5 times the average discharge). This hydrological situation occurred at the same time as the release of a high volume of sewage. The sewage release ended four days before the sampling time. The greatest similarity in the overall EPT community of samples to the samples from Site 3 is shown in the overlap of the ellipses in Figure S4. This similarity was subsequentially lost when the heavy rains stopped: the greatest dissimilarity between samples from individual sites was observed in August, followed by September. Nonzero values of high 3-day sewage volumes were linked to EPT compositions sampled in both these months, as well as in May. The significant linkage between 3-day sewage volumes and the composition of feeding strategies, as a result of FS in RDA (Figure 5B), supports the presumption that the population of campodeiform net-spinning caddisflies is directly or indirectly dependent on the supply of 'poorly treated' WWs.

#### **5. Conclusions**

This study highlights the impacts on the EPT community of trace metal pollution originating from mining and smelting and effluent from a municipal WWTP, which releases continuously treated WW and, periodical, 'poorly' treated WW. According to our data and previously published studies, trace metal pollution can have quite different consequences for the EPT community than, for example, pesticides and several caddisflies, especially limnephilids, were found to be very abundant at the affected site. Additionally, malformations in caddisfly larvae and mortality of caddisfly pupae were observed. Even though the affected sites were warmer, the absence of ephemerellid, heptageniid, and leptophlebiid mayflies is attributable to trace metal pollution, given the results from previous studies. The factors responsible for the decline in the stonefly community are less clear, and, except for higher water temperature, other kinds of pollutants and habitat degradation can be assumed to have a negative impact. Nevertheless, *L. fusca* gr. and *L. geniculata* can be considered tolerant to trace metal pollution but sensitive to WW pollution or its combination with trace metals. The municipal WW combined with trace metal contamination led to the greatest changes in the EPT community, including changes in the composition of feeding strategies. Caddisfly passive filter feeders and predators—mainly hydropsychids, polycentropodids, and rhyacophilids—dominated the most polluted environment. EPT shredders and collectors/gatherers dominated in unaffected sites. Our results demonstrate how the EPT community reacts to human-induced gradients in natural environments, which is important knowledge for assessing the implication of planned or current human disturbances in natural or cultural landscapes.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/biology11050648/s1, Materials and Methods; Results, Figure S1: photographic documentation of Site 1 and Site 2, Figure S2: photographic documentation of Site 3, Site 4 and Auxiliary Site 4, Figure S3: summer temperature regime monitored by data loggers (TFA), Table S1: physicochemical parameters of water at longitudinal profile, Table S2: parameters of discharge at monitored sites, Table S3: contamination of water and sediments by trace metals at longitudinal profile, Table S4: contamination of water by trace metals at Auxiliary Site 3 and Site 4 during 2018–2020, Table S5: pesticides and pharmaceutical active compounds in water samples taken at the longitudinal profile, Table S6: list of pesticide and pharmaceutical active compounds analysed in water samples with particular limits of quantification (LOQs), Figure S4: classified sample diagram resulted from choriotopic compositions using the RDA (sampling sites), Table S7: list of EPT taxa detected at given sites and times, Figure S5: classified CCA sample diagram resulted from species compositions and CCA species pie chart (sampling sites), Figure S6: CCA Species pie chart (sampling times), Figure S7: plotted significant averaged feeding strategies response curves, Figure S8: plotted significant averaged feeding strategies response curves against '3-d Sew' and '15-d Sew', Table S8: the statistics of relationship strengths for plots in Figures S7 and S8, Table S9: abundance and biomass

of fishes detected at the particulate sites during 20 May and 8 October in 2020, Figure S9: plotted significant relationship between EPT family richness and Metal index. References [14,76–78] are cited in the supplementary materials.

**Author Contributions:** Conceptualization, M.L. and M.B.; methodology, M.L.; formal analysis, M.L.; investigation, M.L., J.C., P.N., F.L., and M.B.; writing—original draft preparation, M.L.; writing— ˇ review and editing, M.B. and P.N.; supervision, M.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** Project was supported by the Czech Science Foundation (20-16111S) and by the Ministry of Education, Youth and Sports of the Czech Republic (project CENAKVA no. LM2018099).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** We would like to thank Michael Lockwood for English language editing, Wei Guo, Pavel Franta, and Jiˇrí Jakš for their contribution during field works, Olga Valentová for analysis of water hardness, the Czech Hydrometeorological Institute (Prague, Czech Republic), Povodí Vltavy, State Enterprise (Prague, Czech Republic), and 1. SˇcV JSC (Pˇríbram, Czech Republic) for providing environmental data, and finally to the Local Organisations of the Czech Anglers Union Hoˇrovice and Pˇríbram for enabling the fish stock monitoring.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Zhuo Li 1,2, Xiaowei Liu 1,2, Minghui Zhang 1,2 and Fu Xing 1,2,\***


**Simple Summary:** Understanding relationships between biodiversity and ecosystem functions is important in the context of global plant diversity loss. We evaluated the relationships between soil bacterial and fungal diversity, rare microbial taxa, and soil multifunctionality in a semi-arid grassland with varied plant diversity levels. The fungal richness rather than bacterial richness was positively related to soil multifunctionality. The relative abundance of saprotrophs was positively correlated with soil multifunctionality, and the relative abundance of pathogens was negatively correlated with soil multifunctionality. Furthermore, the rare fungal taxa played a disproportionate role in regulating soil multifunctionality. The shift of plant biomass allocation patterns increased plant below-ground biomass in the high diversity plant assemblages, which can alleviate soil microbial carbon limitations and enhance the fungal richness, thus promoting soil multifunctionality. This study provides a new perspective for evaluating the relative roles of fungal and bacterial diversity in maintaining multiple soil functions under global plant diversity loss scenarios.

**Abstract:** Loss in plant diversity is expected to impact biodiversity and ecosystem functioning (BEF) in terrestrial ecosystems. Soil microbes play essential roles in regulating ecosystem functions. However, the important roles and differences in bacterial and fungal diversity and rare microbial taxa in driving soil multifunctionality based on plant diversity remain poorly understood in grassland ecosystems. Here, we carried out an experiment in six study sites with varied plant diversity levels to evaluate the relationships between soil bacterial and fungal diversity, rare taxa, and soil multifunctionality in a semi-arid grassland. We used Illumina HiSeq sequencing to determine soil bacterial and fungal diversity and evaluated soil functions associated with the nutrient cycle. We found that high diversity plant assemblages had a higher ratio of below-ground biomass to aboveground biomass, soil multifunctionality, and lower microbial carbon limitation than those with low diversity. Moreover, the fungal richness was negatively and significantly associated with microbial carbon limitations. The fungal richness was positively related to soil multifunctionality, but the bacterial richness was not. We also found that the relative abundance of saprotrophs was positively correlated with soil multifunctionality, and the relative abundance of pathogens was negatively correlated with soil multifunctionality. In addition, the rare fungal taxa played a disproportionate role in regulating soil multifunctionality. Structural equation modeling showed that the shift of plant biomass allocation patterns increased plant below-ground biomass in the highly diverse plant plots, which can alleviate soil microbial carbon limitations and enhance the fungal richness, thus promoting soil multifunctionality. Overall, these findings expand our comprehensive understanding of the critical role of soil fungal diversity and rare taxa in regulating soil multifunctionality under global plant diversity loss scenarios.

**Keywords:** plant diversity loss; soil multifunctionality; fungal diversity; saprotrophic fungi; rare microbial taxa

**Citation:** Li, Z.; Liu, X.; Zhang, M.; Xing, F. Plant Diversity and Fungal Richness Regulate the Changes in Soil Multifunctionality in a Semi-Arid Grassland. *Biology* **2022**, *11*, 870. https://doi.org/10.3390/ biology11060870

Academic Editors: Daniel Puppe, Panayiotis Dimitrakopoulos and Baorong Lu

Received: 7 May 2022 Accepted: 3 June 2022 Published: 6 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Soil supports a wide range of ecological functions and services that are important to human welfare, and soil multifunctionality can be used as a comprehensive index to evaluate soil quality [1–4]. Multifunctionality is an essential biological and management concept that describes the ability of an ecosystem to maintain multiple ecological functions simultaneously [1,2,5]. As significant drivers of ecosystem functions, soil microbial communities can be regulated by plant diversity via changing the nutrient availability and microenvironmental conditions [6–8]. However, global climate changes and intensive anthropogenic activities have led to the loss of plant diversity in grassland ecosystems [2,9–11]. Therefore, there remains an urgent need to understand the mechanisms by which soil microorganisms drive soil multifunctionality in the context of global plant diversity loss.

Below-ground microbes comprise a large portion of global terrestrial ecosystem life's genetic diversity and support important ecosystem functions and services, such as biomass production, nutrient cycling, and decomposition of organic matter [10,12,13]. Additionality, soil microorganisms are an essential link between the above-ground and below-ground components of terrestrial ecosystems [14–17]. Bacteria and fungi constitute vital parts of the soil microbiome [18]. Nevertheless, bacteria and fungi play different roles in regulating ecosystem processes [14]. In detail, bacteria are the essential drivers of soil N cycle processes, such as N2 fixation, nitrification, and denitrification [19–21]. However, other studies have suggested that soil fungal communities might dominate in driving soil functions [14]. Soil fungi can decompose recalcitrant plant litter efficiently and form links with roots to capture and transport subsurface carbon (C) [22]. For instance, saprophytic fungi are primary mediators for the decomposition of plant litter, and their mycelium networks across the soil litter interface and networks are a highly dynamic channel through which nutrients can be easily distributed [23]. Despite the above studies, most research has concentrated on how soil microbes affect individual or specific nutrient cycling functions. However, few studies have evaluated the relative significance of bacterial and fungal diversity in regulating soil multifunctionality in the context of plant diversity in a semi-arid grassland.

Furthermore, the effects of rare microbial taxa on soil functions have often been ignored in most previous soil multiple functions studies because they have mainly concentrated on predominant taxa; a majority of rare microbial taxa were frequently eliminated from original datasets. However, rare microbial taxa have important ecological roles because these may be activated to perform critical functions under environmental changes [24,25]. For instance, rare taxa were found to play essential roles in nutrient cycling after disturbances in aquatic ecosystems [26]. Rare microbial taxa have also been used as indicators of changes in ecosystem functions under long-term greenhouse cultivation conditions in subtropical agricultural soils [27]. Therefore, rare microbial taxa can represent a microbial "seed bank" that may be activated when the environment is disturbed [24,28]. The intensity and frequency of plant diversity loss due to global climate change are expected to increase in the future [9,29,30]; therefore, it is of great significance to link the ecological functions of rare microbial taxa with soil multifunctionality. Unfortunately, relevant knowledge is still scarce, especially in semi-arid grassland ecosystems. Such knowledge is critical for developing a management framework to maintain rare taxa involved in functionality and decrease the impact of future climate change on a semi-arid grassland.

Loss in plant diversity of terrestrial ecosystems might aggravate C and other nutrient limitations [30–32]. Heterotrophic organisms depend on plants to obtain C substrates; therefore, there is ample evidence that soil microbes are usually limited by C [33]. Microorganisms play an essential role in ecosystem functions [6,18,34]; thus, more research is undoubtedly required to evaluate how the shift in microbial nutrient limitations affects microbes under changes in plant diversity. Soil enzymes produced by microbes transform substrates in soil C and nutrient cycles [35]. Thus, enzymatic stoichiometry provides a method for this study [35]. Several studies have shown that less diverse plant communities allocate fewer photosynthates to below-ground ecosystems by decreasing the exudation of root exudates [32]. Root exudates are an important C source for microbial growth [36]; therefore, a less diverse plant community might intensify microbial C limitation [37]. C limitation could affect microbial biosynthesis processes [35]; thus, lower below-ground biomass allocation would reduce soil functions such as microbial growth and respiration by increasing microbial C limitation [36,38]. The above results show that plant diversity affects soil microbial communities by trophic interactions and altering soil biochemical processes [12,13]. Thus, the loss of plant diversity may alter the relationships between plant and soil microbial communities. Considering that soil microbial communities mediate ecosystem functions, it is essential to elucidate the mechanisms by which plant-diversity-induced changes in microbial C limitation regulate soil multifunctionality through impacting soil microorganisms in a meadow grassland.

In addition, ecosystem multifunctionality is different from ecosystem services multifunctionality in that the former represents the overall performance of an ecosystem without considering stakeholders; whereas the latter is defined as the ability of ecosystems to synergistically provide various ecosystem functions and services that translate into a variety of social benefits and welfare [1]. Therefore, investigations of ecosystem service multifunctionality have significant importance in future semi-arid grassland ecosystems management and sustainable development. In this study, a meadow grassland in northeast China was selected to discuss the soil multifunctionality and microbial driving mechanisms under a plant diversity loss scenario. The objectives of this study were to 1) evaluate the important roles of soil bacterial and fungal diversity and rare taxa in driving soil multifunctionality under different plant diversity conditions; 2) reveal the underlying mechanisms by which plant-diversity-induced alterations of microbial C limitation regulate soil multifunctionality through soil microorganisms; 3) provides a scientific reference for the evaluation of semi-arid grassland ecosystem services multifunctionality.

#### **2. Materials and Methods**

#### *2.1. Study Site*

The experiment was conducted at the Tongyu Observatory in a semi-arid climate and environment (44◦ 25 N, 122◦ 52 E). This site is located within the Songnen grassland ecosystem of northeastern China (Figure 1a). The vegetation is a meadow steppe, situated at the eastern end of the Eurasian grassland belt with *Leymus chinensis* (Trin.) Tzvelev (https://www.ipni.org/) (accessed on 1 June 2022) is the dominant species. The soil texture is clay loam, according to the International Society of Soil Science Standard. The study area has a temperate continental monsoon climate. The mean annual precipitation is ~404 mm, and more than 80% of the rainfall is concentrated during the growing season (from May to September). The mean annual temperature is approximately 5.7 ◦C.

#### *2.2. Experimental Design and Sampling*

Within approximately 1 km<sup>2</sup> of the study area, we randomly selected three single plant species patches i.e., *Carex duriuscula* C.A.Mey., *Lespedeza hedysaroides* (Pall.) Kitag., and *Calamagrostis rigidula* A.I.Baranov and Skvortsov assemblages and three multiple species coexisting plant assemblages with different dominant *Lespedeza daurica* (Laxm.) Schindl., *L. chinensis*, and *Hierochloe glabra* Trin. naturally existing in the grassland. Single species assemblages and multiple species assemblages represent low and high diversity levels, respectively. The plant species composition is shown in Table S1. A vegetation survey was conducted in mid-August 2018 when the standing biomass reached its maximum. For the six study sites selected above, the interval between every two sites was greater than 100 m. For plant and soil sampling, five plots (1 m × 1 m) with an interval of more than 10 m were randomly selected at each study site (Figure 1b). In this study, soil characteristics of low diversity plant assemblages and high diversity plant assemblages were similar (*p* > 0.05) (Figure S1); thus, soils were comparable across different plant diversity plots. We selected six sites with similar soil characteristics but different levels of plant diversity. It is reasonable to obtain data based on one-time sampling to explore the impact of different

levels of plant diversity on soil, microorganisms, and soil functions. Additionally, this method has been widely adopted by many studies in the field of ecology [32,39].

**Figure 1.** Locations of the study site (**a**) and experimental design for plant and soil sampling (**b**) in the Songnen grassland, China. Three single plant species patches, i.e., *Carex duriuscula* C.A.Mey. (CD), *Lespedeza hedysaroides* (Pall.) Kitag. (LH), and *Calamagrostis rigidula* A.I.Baranov and Skvortsov (CR) assemblages and three multiple species coexisting plant assemblages with different dominant *Lespedeza daurica* (Laxm.) Schindl. (LD.ass), *Leymus chinensis* (Trin.) Tzvelev (LC.ass), and *Hierochloe glabra* Trin. (HG.ass) naturally existing in the grassland. For the six study sites selected above, the interval between every two sites was greater than 100 m. For plant and soil sampling, five plots (1 m × 1 m) with an interval of more than 10 m were randomly selected at each study site. Five soil cores were randomly selected in each plot by a soil auger.

All the plant species were identified and recorded in each plot; then, soil and plant samples were collected simultaneously from each plot. Five soil cores (0–10 cm) were randomly selected in each plot by a soil auger (5 cm diameter), then mixed to form a composite sample. Next, the soil samples were sieved (2 mm) to remove any roots and stones. We harvested above-ground living plants to evaluate each plot's above-ground biomass (AGB). The below-ground biomass (BGB) of each plot was measured by root biomass from three soil cores with a diameter of 10 cm and a depth of 30 cm. All the soil samples were placed in well-sealed zippered bags and transported to the laboratory within 24 h in cooling boxes.

#### *2.3. Plant and Soil Samples Analysis*

The plant AGB and BGB were determined by plant materials heated at 105 ◦C for 30 min to rapidly cease metabolic activities [40] and then oven-dried at 65 ◦C for 48 h until the weight remained constant [41]. The soil samples were divided into three subsamples for physiochemical analyses and the assessment of microbial communities. One aliquot of soil was stored at 4 ◦C and used to determine the soil water content (SWC), soil available N (AN) content (NH4 +-N and NO3 −-N), soil net nitrification rate (Nn), and the net N mineralization rate (Nm). Another aliquot of soil was stored at −80 ◦C and used for the analysis of high-throughput sequencing and to assess activities of α-1,4-glucosidase (αG), β-1,4-glucosidase (βG), β-1,4-xylosidase (βX), β-D-cellobiohydrolase (CBH), leucine aminopeptidase (LAP), β-1,4-N-acetylglucosaminidase (NAG), and alkaline phosphatase (ALP). Finally, the third aliquot of soil was air-dried at room temperature for soil pH, electrical conductivity (EC), soil total organic carbon (TOC), total soil nitrogen (TN), total soil phosphorus (TP), and available soil phosphorus (AP) analyses.

The soil pH was detected in deionized water (soil: water = 1:5, w:v) with a portable pH meter (Leichi PHBJ-260, Shanghai, China), and the soil EC was measured with an electronic conductivity meter (Leici DDS307, Shanghai, China) [42]. The SWC was measured as the weight loss recorded after the fresh soils had been oven-dried to a constant weight at 105 ◦C [10]. The soil TN content was measured with an elemental analyzer (vario EL cube, Elementar, Langenselbold, Germany). Briefly, 30 mg of the air-dried soil was put into a tinfoil cup and tightly wrapped in soil and then put into an automatic sampling tray; the results were analyzed on the machine. Soil TOC was determined with an elemental analyzer (Isoprime 100, Isoprime Ltd., Manchester, UK). The NH4 +-N and NO3 −-N concentrations were measured by extraction with 2 M KCl (soil: water = 1: 5, w:v) and analyzed with a continuous flow analyzer (Futura, Alliance-AMS) [41]. The Nn and Nm were calculated according to alterations in the concentrations of NH4 +-N and NO3 −-N before and after incubation [43]. For soil TP, 0.5 g air-dried soil and ball milling were digested with HClO4 (7.7 mL, 75%) at 203 ◦C for 75 min [44]. Soil AP was extracted with 0.5 M NaHCO3, and the molybdenum blue colorimetric method was used to analyze [45]. All enzyme activities were measured using a 4-methylumbelliferyl (MUB) substrate, except for 7-amino-4-methyl-coumarin (7-AMC) for the LAP [46]. Seven enzyme activities were determined in black 96-well microplates. Four replicate wells were set up for each enzyme test sample (200 μL slurry + 50 μL substrate) and corresponding substrate control (200 μL buffer + 50 μL substrate). Four replicate wells were set up for standard fluorescence (200 μL buffer + 50 μL standard), slurry control (200 μL slurry + 50 μL buffer), and quench standards (200 μL slurry + 50 μL standard). The assay plate was incubated in a dark environment at 25 ◦C for 3 h. Fluorescence was measured using a microplate fluorometer (TECAN infinite F200, Tecan Group, Switzerland) with excitation and emission filters of 360 nm and 460 nm, respectively.

#### *2.4. Assessment of Microbial Communities*

According to the manufacturer's protocol, total genomic DNA was extracted from 0.3 g of soil with the PowerSoil DNA Isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA). The concentration and purity of the extracted DNA were measured with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA ). DNA quality was evaluated by 1% agarose gel electrophoresis. We acquired information on the diversity and composition of soil bacteria and fungi by performing 16S rRNA and ITS genes amplicon sequencing and the 338F/806R (5 -ACTCCTACGGGAGGCAGCA-3 , 5 - GGACTACHVGGGTWTCTAAT-3 ) and ITS1F/ITS2R (5 -CTTGGTCATTTAGAGGAAGTAA-3 , 5 -GCTGCGTTCTTCATCGATGC-3 ) primer pairs, respectively. The total volume of the PCR amplification for bacteria and fungi was 50 μL, including 10 μL buffer, 0.2 μL Q5 High-Fidelity DNA Polymerase, 10 μL High GC Enhancer, 1 μL dNTP, a 10 μM concentration of each primer, and 60 ng genome DNA. The thermal cycling conditions were 95 ◦C denaturation for 5 min, 15 cycles at 95 ◦C for 1 min, 50 ◦C for 1 min, 72 ◦C for 1 min, and 72 ◦C final extension for 7 min. Ultimately, we quantified all PCR products by Quant-iT™ dsDNA HS Reagent and pooled them together. Furthermore, we conducted a high-throughput sequencing analysis of 16S rRNA and ITS genes by the Illumina

Hiseq 2500 platform (2 × 250 paired ends) at Biomarker Technologies Corporation (Beijing, China). Bacterial and fungal sequences were quality-filtered with the QIIME software package and merged using the FLASH software package [47]. The operational taxonomic units (OTUs) were defined by 97% similarity. Representative sequences of bacteria and fungi were annotated with the SILVA and UNITE databases [48]. The final total data set retained 1852 and 1325 OTU numbers and 2,057,790 and 2,071,174 clean reads for bacteria and fungi, respectively. The raw bacterial and fungal reads were deposited into the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under accession numbers PRJNA810930 and PRJNA810946, respectively.

#### *2.5. Definition of Abundant and Rare Taxa*

Previous studies have widely used relative abundance as a metric to describe microbial taxa or OTUs in their environment. Thus, relative abundance is useful for classifying abundant or rare taxa in microbial communities [24,49]. We identified relative abundance thresholds as 1% for abundant taxa and 0.01% for rare taxa [26,27,50]. These classifications ignored intermediate taxa (relative abundance between 0.01 and 1%) and oscillatory taxa (rare and abundant under different environment conditions) [27]. All OTUs were classified into six categories based on the criteria used in recent studies [26,27,50]: always abundant taxa (AAT, relative abundance ≥ 1%); conditionally abundant taxa (CAT, relative abundance ≥ 0.01% in all samples and ≥ 1% in some samples); always rare taxa (ART, relative abundance < 0.01% in all samples); conditionally rare taxa (CRT, relative abundance < 0.01% in some samples but never >1% in any samples); moderate taxa (MT, relative abundance between 0.01% and 1% in all samples); and conditionally rare and abundant taxa (CRAT, relative abundance < 0.01% in some samples and ≥ 1% in some samples). Finally, AAT and CAT were classified as abundant taxa, and ART and CRT were classified as rare taxa [26,27,50].

#### *2.6. Assessment of Multifunctionality*

Soil multifunctionality was quantified based on ten soil functions widely adopted in previous studies [51–55]. These parameters are related to organic matter decomposition, climate regulation, and nutrient cycling, including αG, βG, βX, CBH, LAP, NAG, Nn, Nm, ALP, and AP. Of them, αG, βG, βX, and CBH represent the C cycle, LAP, NAG, Nn, and Nm represent the N cycle, and ALP and AP represent the P cycle (Table S2). The averaging and threshold approaches have commonly been used to estimate the relationship between biodiversity and multifunctionality [2,34,53,56–58]. In this study, we evaluated the soil multifunctionality index using three methods, which contained averaging approach [56], single-threshold, and multiple-threshold approach [57,58]. To obtain an average soil multifunctionality index, each function was standardized by Z-score transformation and then averaged [2,56]. We used a single-threshold approach to calculate the number of soil functions exceeding a given threshold (25%, 50%, 75%, and 90%) [2,54]. However, the multiple-threshold method can demonstrate the effect of bacterial and fungal richness on soil multifunctionality in the full threshold range.

#### *2.7. Enzymatic Stoichiometry*

We constructed the vector analysis of soil enzymatic stoichiometry to evaluate soil microbial C limitation levels [59]. A longer vector length represents a stronger C limitation.

$$\text{Vector length (unitless)} = \sqrt{(\ln \beta \text{G/h} (\text{NAG} + \text{LAP}))^2 + (\ln \beta \text{G/h} \text{LAP})^2}$$

#### *2.8. Statistical Analysis*

We used an independent sample *t*-test to evaluate the differences in plant and soil properties and soil multifunctionality between different levels of plant diversity. When the data did not meet the requirements for normality (Shapiro Wilk test) and homogeneity of variance (Bartlett test), log transformations were performed on the data. We calculated

the microbial richness with the "*vegan*" package [60]. Ordinary least squares (OLS) linear regression was conducted to explore the relationships between microbial C limitation and microbial richness; OLS linear regression also was conducted to explore the relationships between bacterial richness, fungal richness, relative pathogen abundance, saprotrophic relative abundance, relative symbiont abundance, and soil multifunctionality [14,15,34,61]. Three methods of soil multifunctionality were evaluated with the "multifunc" package [58]. According to recent studies [14,47,61,62], functions were obtained based on fungal taxa by the FUNGuild database (http://www.stbates.org/guilds/app.php) (accessed on 6 January 2019). Random forest modeling analysis was applied with the "randomForest" package [63] to identify the major statistically significant microbial taxa (OTU level) in impacting soil multifunctionality. Then, seven microbial taxa, including four bacterial taxa and three fungal taxa, were selected in the random forest modeling. The "A3" package [64] was used to assess the importance of each predictor on soil multifunctionality. In addition, we used piecewise structural equation modeling (SEM) with the "piecewiseSEM" [65], "nlme" and "lme4" packages [66] to investigate the potential direct and indirect effects of plant diversity on soil multifunctionality. The prior model was constructed based on current ecological knowledge of the grassland ecosystems [7,36,53,67] to evaluate the link between soil biodiversity and multifunctionality (averaging), assuming that plant diversity altered plant biomass allocation patterns and microbial C limitation, thus promoting soil microbial richness and soil multifunctionality (Figure S2a). Fisher's C test (0 ≤ Fisher's C ≤ 2 *df* and 0.05 < *p* ≤ 1.00) was used to verify the rationality of the modeling results. Among significant models, the one with the lowest Akaike information criterion (AIC) value was selected for the final SEM analysis [68] (Figure S2b–d). All statistical analyses and visualizations of this study were performed in R v.4.1.2 software.

#### **3. Results**

#### *3.1. Plant and Soil Properties and Soil Multifunctionality*

The ratio of BGB to AGB was significantly higher with high plant diversity than with low plant diversity (Figure 2a; *p* < 0.001). Microbial C limitation was higher with low diversity than with high diversity (Figure 2b; *p* < 0.01). The results showed that plant diversity significantly altered soil multifunctionality (Figure 2c), i.e., high diversity significantly increased soil multifunctionality compared with low plant diversity (Figure 2c; *p* < 0.001). SWC, soil pH, EC, AN, and TOC showed no significant differences between the two plant diversities (Figure S3a–d,g). Soil TN and TP were higher with low diversity than with high diversity (Figure S3e,f; *p* < 0.05 and *p* < 0.01). In addition, the soil C/N was much higher with high diversity than with low diversity (Figure S3h; *p* < 0.001).

#### *3.2. Soil Microbial Community Composition*

The rarefaction curves of OTU richness of bacterial and fungal communities almost approached saturation (Figure S4a,b), showing that the amount of data of sequenced reads were reasonable. The bacterial OTUs were assigned to 23 phyla, 76 classes, 110 orders, 188 families, 271 genera, and 226 species; similarly, the fungal OTUs were assigned to 10 phyla, 23 classes, 51 orders, 103 families, 177 genera, and 152 species. The bacterial communities were dominated by Acidobacteria, Proteobacteria, and Actinobacteria (Figure S5). A significantly higher relative abundance of Acidobacteria was observed with low plant diversity than with high plant diversity (Figure S5a; *p* < 0.001). The relative abundance of Actinobacteria and Chloroflexi was higher with high diversity than with low diversity (Figure S5b,d; *p* < 0.001). The fungal communities were dominated by Ascomycota and Basidiomycota (Figure S5). The bacterial communities were dominated by the genera *RB41*, *uncultured\_bacterium\_c\_Subgroup\_6*, and *uncultured\_bacterium\_f\_Gemmatimonadaceae* (Figure 3a). However, the fungal communities were dominated by the genera *Ceratobasidium*, *Mortierella*, and *Cladosporium* (Figure 3b). The richness (OTU richness) of fungal and bacterial were significantly greater with high diversity than with low diversity (Figure S6; *p* < 0.001, *p* < 0.05).

**Figure 2.** The ratio of BGB to AGB (**a**), microbial C limitation (**b**), and soil multifunctionality (**c**) in response to plant diversity. \*\* *p* < 0.01 and \*\*\* *p* < 0.001 (*t*-test). AGB: above-ground biomass, BGB: below-ground biomass.

**Figure 3.** Relative abundance of bacterial (**a**) and fungal (**b**) dominant genera in different plant diversity levels, "Others" represents all genera with relative abundance < 1%.

#### *3.3. Microbial Taxa in Relation to Soil Multifunctionality*

The OLS regression models showed that microbial C limitation was significantly and negatively related to fungal richness (Figure 4b; *R*<sup>2</sup> = 0.309, *p* < 0.001). Furthermore, the results showed that fungal richness was significantly and positively associated with average soil multifunctionality (Figure 5b; *R*<sup>2</sup> = 0.233, *p* = 0.003). However, bacterial richness showed no correlation with soil multifunctionality (Figure 5a). Based on the positive correlation between fungal richness and soil multifunctionality, we aimed to identify the relationships between fungal guilds and soil multifunctionality. A negative linear correlation was found between the relative abundance of pathogens and soil multifunctionality (Figure 5c; *R*<sup>2</sup> = 0.140, *p* = 0.023). However, we found a significantly positive relationship between the relative abundance of saprotrophs and soil multifunctionality (Figure 5d; *R*<sup>2</sup> = 0.165, *p* = 0.014).

**Figure 4.** Relationships between soil microbial carbon (C) limitation and bacterial richness (**a**) and fungal richness (**b**). The black lines represent the fitted ordinary least squares (OLS) linear regressions. Red circles represent high plant diversity, while blue ones indicate low plant diversity.

Soil multifunctionality was positively and significantly related to fungal richness rather than bacterial richness. This result was confirmed using the single-threshold method (Figure S7) and the multiple-threshold method (Figure 6). Threshold methods were also used to assess whether multiple functions were performed at high levels simultaneously. Specifically, we performed single-threshold analysis; positive and significant correlations were found between the bacterial richness and soil multifunctionality at given thresholds of 75% and 90% (Figure S7g,h; *R*<sup>2</sup> = 0.108, *p* = 0.042, *R*<sup>2</sup> = 0.105, *p* = 0.044). However, there were significant and positive correlations between the fungal richness and soil multifunctionality at given thresholds of 25%, 50%, 75%, and 90% (Figure S7a–d; *R*<sup>2</sup> = 0.220, *p* = 0.005, *R*<sup>2</sup> = 0.111, *p* = 0.039, *R*<sup>2</sup> = 0.313, *p* < 0.001, *R*<sup>2</sup> = 0.292, *p* = 0.001). However, the multiplethreshold method does not require a threshold to be set and studies a continuous threshold gradient (Figure 6a,b). The minimum threshold (Tmin) for fungi was 16%, the lowest threshold at which fungal richness began to have a positive effect on soil multifunctionality. When the threshold was 36%, the maximum richness effect (Rmde) was 0.016, i.e., the relationship strength of fungal richness with the strongest positive effect, indicating that the addition of one species of fungus could increase the function by 0.016 (Figure 6d).

**Figure 5.** Relationships between soil multifunctionality and bacterial richness (**a**), fungal richness (**b**), pathogens relative abundance (**c**), saprotrophs relative abundance (**d**), and symbionts relative abundance (**e**). The black lines represent the fitted ordinary least squares (OLS) linear regressions. Red circles represent high plant diversity, while blue ones indicate low plant diversity.

**Figure 6.** Effects of bacterial (**a**) and fungal (**b**) richness on the number of functions above thresholds. Lines represent the slope between soil microbial richness and the number of functions greater than or equal to a threshold value ranging from 5% to 99% of the maximum for each function. The dotted curves indicate the changes in the number of functions per unit increment of the richness of bacteria (**c**) and fungi (**d**). Tmin is the minimum threshold that soil multifunctionality becomes influenced by changes in microbial richness, and Rmed is the realized maximum effect of richness on soil multifunctionality.

#### *3.4. Microbial Taxa Predicting Soil Multifunctionality*

Bacterial rare taxa accounted for 81.64% of total OTUs and 40.90% of the relative abundance, respectively (Figure 7a). In comparison, abundant bacterial taxa accounted for 1.13% of total OTUs and 17.55% of the relative abundance, respectively (Figure 7a). In addition, bacterial conditionally rare and abundant taxa accounted for 0.11% of total OTUs and 0.67% of the relative abundance, respectively (Figure 7a). Bacterial moderate taxa accounted for 17.12% of total OTUs and 40.86% of the relative abundance, respectively (Figure 7a). Fungal rare taxa accounted for 81.66% of total OTUs and 23.93% of the relative abundance, respectively (Figure 7a). Fungal abundant taxa accounted for 0.91% of total OTUs and 11.52% of the relative abundance, respectively (Figure 7a). Finally, fungal conditionally rare and abundant taxa accounted for 17.28% of total OTUs and 64.02% of the relative abundance, respectively (Figure 7a). Fungal moderate taxa accounted for 0.15% of total OTUs and 0.54% of the relative abundance, respectively (Figure 7a). Random forest modeling results indicated that fungal ART (*p* = 0.012) was the most important microbial taxa in predicting soil multifunctionality (Figure 7b).

**Figure 7.** The percentage of OTUs and relative abundance of six categories of microbial taxa (**a**). The random forest regression model indicates the main microbial drivers of soil multifunctionality (**b**). MSE is the mean square error. This accuracy importance measure is computed for each tree and averaged over the forest (5000 trees) and determined by the increase in the MSE. \* *p* < 0.05 on the bar shows that the associated microbial taxa have a significant effect on soil multifunctionality. ART: always rare taxa; CRT: conditionally rare taxa; AAT: always abundant taxa; CAT: conditionally abundant taxa; CRAT: conditionally rare and abundant taxa; MT: moderate taxa.

#### *3.5. Direct and Indirect Effects of Plant Diversity on Soil Multifunctionality*

The SEM analysis explained 53% of the total variation in soil multifunctionality (Figure 8a). Plant diversity indirectly affected soil multifunctionality by influencing the ratio of BGB to AGB, microbial C limitation, and fungal richness but not bacterial richness (Figure 8a and Figure S3). Plant diversity, the ratio of BGB to AGB, and fungal richness displayed positive effects on soil multifunctionality (Figure 8b). However, microbial C limitation exhibited a negative effect on soil multifunctionality (Figure 8b).

**Figure 8.** The piecewise structural equation model was used to test the direct and indirect relationships between plant diversity and soil multifunctionality (**a**). The corresponding values of the solid line arrows are the width of normalized path system arrows that reflect the size of the normalized path coefficient. The blue and red arrows show significant positive and negative correlations, respectively (\* *p* < 0.05, \*\* *p* < 0.01, and \*\*\* *p* < 0.001). Dashed lines show non-significant relationships. The values above each variable represent the explanatory degree (*R*2) of each variable in the model. Standardized total effects (direct plus indirect effects) derived from the piecewise structural equation model (**b**).

#### **4. Discussion**

Our study provides evidence that fungal richness rather than bacterial richness was significantly related to soil multifunctionality in a semi-arid grassland in northeast China. Saprotrophic fungi and rare fungal taxa were essential for maintaining soil functions. Furthermore, under high diversity plant assemblages, the changes in plant biomass allocation patterns increased plant below-ground biomass, which can alleviate microbial C limitation and thus enhance the fungal richness, ultimately promoting soil multifunctionality. Such results suggest that the above-ground and below-ground biodiversity, as well as rare fungal taxa, are vital to maintaining ecosystem functions in a semi-arid grassland ecosystem.

Soil microorganisms are some of the most sensitive components to precipitation change in semi-arid grasslands [69]. Therefore, biodiversity in semi-arid zones can be seriously compromised if rainfall changes. For example, reducing precipitation significantly increased fungal diversity [70]. However, one study showed a significant positive correlation between fungal diversity and precipitation [71]. Previous studies have suggested that bacterial communities are more sensitive to changes in precipitation than fungal communities in a semi-arid grassland [72]. However, we found that fungal richness and saprotrophic fungi were principal biotic factors in regulating soil multifunctionality in a semi-arid grassland. Fungi can create an environment around themselves by secreting polysaccharides from their hyphae to prevent dehydration [73]. Moreover, substrate diffusion restrictions might force the soil fungal mycelium network to expand, which is conducive to absorbing water and nutrients [74]. Therefore, fungi are more resistant to drought than bacteria, which may play an important role in arid and semi-arid grasslands [75]. For instance, soil fungi are crucial to organic matter decomposition, and root-associated fungi are important regulators of ecosystem C dynamics [22]. Moreover, fungal richness was significantly positively correlated with denitrification activity, indicating that fungi could promote soil N cycling [14]. Our results also showed that saprotrophic fungi significantly affected soil multifunctionality (Figure 5d). Saprotrophic fungi mainly grow throughout the soil litter interface and are involved in plant litter decomposition [76]. Previous studies have also suggested that saprotrophic fungi might affect soil C storage through interactions with ectomycorrhizal fungi [77]. Moreover, free-living saprotrophs generally play an essential ecological role in dead plant material because they can derive C by propagules or hyphae from dead organic material [78]. Saprotrophic fungi can also obtain fresh nutrients from recalcitrant organic polymers using extracellular enzymes and the nonenzymatic Fenton reaction [79,80]. For example, *Podospora anserina* is a potent plant biomass degrader and efficiently utilizes lignocellulose as a C source through dedicated lignin degradation enzymes [81,82]. Therefore, the ability of saprophytic fungi to translocate C resources implies that saprophytic fungi are essential agents of nutrient redistribution in soils [83]. Furthermore, the hyphal network formed by saprotrophic fungi is involved in the formation of soil aggregates, which is of great significance for soil water retention and erosion resistance [84]. Collectively, the present study suggested the non-negligible roles played by the fungal richness and saprophytic fungi in regulating soil functions.

An interesting result was that fungal ART was identified as the main driver of soil multifunctionality (Figure 7b). Previous studies have suggested that rare fungal community composition and functional guilds are more stable than those of the abundant taxa under certain conditions [85]. In other words, rare taxa contribute to maintaining the microbiome's function under environmental stress because some may be highly resistant to stress [86]. For example, a recent study found that rare microbial taxa might modulate the adverse effects of semi-arid grassland degradation drivers such as vegetation loss and eutrophication on soil organic matter decomposition [24]. Furthermore, rare taxa contributed more to soil C and N cycling and crop yield than the abundant taxa [25]. Thus, rare microbial taxa play a unique role in maintaining ecosystem functions [28,87]. However, our understanding of rare microbial taxa is still preliminary, and more attention should be given to rare microbial taxa in future studies of biodiversity and ecosystem multifunctionality.

As expected, in line with previous studies, high plant diversity increased soil multifunctionality as compared with low plant diversity [4,56,88]. However, our study elucidated the underlying mechanisms by which high plant diversity could enhance soil multifunctionality through fungal richness. According to current knowledge [89], the possible reason was that the competition for below-ground resources was less than that for above-ground resources in low diversity plant assemblages. Therefore, low diversity plant assemblages might increase the allocation of above-ground biomass to compete for light [89]. Although above-ground biomass was expected to be associated with litter production, below-ground C from rhizodeposition was vital for soil microbial communities [90]. C decomposes from above-ground plant litter to the soil surface; thus, it is not readily available to microorganisms inside the soil [36]. Therefore, a low diversity plant assemblage might intensify microbial C limitation [91], mainly by changes in plant-root-derived substrates [37]. Our results indicated that microbial C limitation significantly affected soil fungal richness but not bacterial richness (Figure 4a,b). A reasonable explanation might be that members of the fungal community are the first consumers of below-ground plant-derived C inputs to the soil [24,92]. Thus, intensifying microbial C limitation reduced soil fungal richness and ultimately decreased soil multifunctionality in the low diversity plant assemblage (Figure 9). However, high diversity plant assemblages increased below-ground biomass allocation (Figure 2a). A possible explanation for the increase in below-ground biomass may be that the interactions between different plant species roots may affect biomass allocation patterns, i.e., interspecific competition can increase plant below-ground biomass [93]. The total soil organic carbon content under high plant diversity was higher than that under low plant

diversity in this study area (Figure S3g). Although this study did not test the contents of soil organic matter mineralization, previous literature has shown that the quantity and quality of soil organic matter are the main driving factors affecting microorganisms [94]. The increasing plant below-ground biomass may improve the quantity and quality of soil organic matter in semi-arid zones [32], and the increased soil C resources provide rich substrates and energy for microorganism growth and basal metabolism [7], thus increasing fungal diversity. In this study, the decreased microbial C limitation favored fungal richness under high plant diversity, thus promoting soil multifunctionality (Figure 9). Furthermore, our results indicate that semi-arid grassland provides essential ecosystem services such as carbon sequestration, climate regulation, and biodiversity protection under high diversity plant assemblages. Our findings suggest that high diversity plant assemblages make semi-arid grassland ecosystems more efficient in regulating ecological processes, which local stakeholders or grassland conservation agencies want. Therefore, ecosystem services through the TESSA methodology (http://tessa.tools/) (accessed on 1 June 2021) [95] should contain future multifunctionality studies. The study of ecosystem services multifunctionality can provide important enlightenment for ecological protection and sustainable development under the increasing pressure of human activities and climate changes. In summary, this study revealed the mechanism of changes in soil multifunctionality in which plant-diversity-induced alterations of microbial C limitation regulate soil multifunctionality by affecting soil fungal richness in a semi-arid grassland.

**Figure 9.** A conceptual framework for understanding the effects of plant and fungal richness on soil multifunctionality in a semi-arid grassland. Changing plant biomass allocation patterns increased the ratio of plant below-ground biomass to above-ground biomass under high diversity plant assemblages, which can alleviate microbial carbon (C) limitation and thus enhance the fungal richness, finally promoting soil multifunctionality. The fungal richness was positively related to soil multifunctionality, but the bacterial richness was not. Saprotrophic fungi and rare fungal taxa were essential for maintaining the soil functions. The blue arrow represents an increase, and the red arrow represents a decrease. + and − describe promotion and inhibition effects. AGB: plant above-ground biomass; BGB: plant below-ground biomass; SMF: soil multifunctionality.

#### **5. Conclusions**

Our results demonstrated that the positive effect of plant diversity on soil multifunctionality was mainly due to the high plant diversity increasing plant below-ground biomass allocation, which alleviated microbial C limitation and favored fungal richness, finally promoting soil multifunctionality. Additionality, the high diversity of plant assemblages enhanced ecosystem services for semi-arid grasslands protection, and plant-fungus relationships were important to improve the assessment of ecosystem services. Instead of bacterial richness, the fungal richness was the crucial biotic predictor of soil multifunctionality in a semi-arid grassland. Furthermore, saprotrophic fungi and rare fungal taxa were major drivers of soil multifunctionality. In conclusion, this study provides a new perspective for evaluating the relative roles of fungal and bacterial diversity and biomass allocation patterns in maintaining soil functions in the context of global plant diversity loss and has important implications for biodiversity conservation and sustainable development in semiarid grasslands. Ecosystem multifunctionality and ecosystem services multifunctionality should be popularized and applied in future multifunctionality studies.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology11060870/s1. Table S1: The plant species composition of the selected low diversity with a single species and high diversity with multiple species assemblages in the Songnen grassland, China; Table S2: Variables used to estimate soil multifunctionality and their importance; Figure S1: Principal component analysis (PCA) to determine the difference of soil characteristics between low diversity plant assemblages and high diversity plant assemblages. The soil characteristics included soil water content, soil pH, soil electrical conductivity, soil available nitrogen, total soil nitrogen, total soil phosphorus, soil total organic carbon, and soil C/N; Figure S2: *Prior* structural equation models (SEM) of hypothetical relationships between plant diversity and soil multifunctionality (a). Models were accepted for the significant Fisher's C test (*p* > 0.05). Among significant models, the one with the lowest Akaike information criterion (AIC) was selected for the final SEM analysis (b–d); Figure S3: Soil water content (SWC) (a), soil pH (b), soil electrical conductivity (EC) (c), soil available nitrogen (AN) (d), soil total nitrogen (TN) (e), soil total phosphorus (TP) (f), soil total organic carbon (TOC) (g), and soil C/N (h) in response to low and high plant diversities. \* *p* < 0.05, \*\* *p* < 0.01, and \*\*\* *p* < 0.001 (*t*-test); Figure S4: Rarefaction curves for the soil bacterial (a) and fungal (b) communities at 97% sequence similarity in different plant diversity assemblages. Low diversity assemblages including single species *Carex duriuscula* C.A.Mey. (CD), *Lespedeza hedysaroides* (Pall.) Kitag. (LH), and *Calamagrostis rigidula* A.I.Baranov and Skvortsov (CR) assemblages. High diversity assemblages including multiple species assemblages with the dominant *Lespedeza daurica* (Laxm.) Schindl. (LD.ass), *Leymus chinensis* (Trin.) Tzvelev (LC.ass), and *Hierochloe glabra* Trin. (HG.ass); Figure S5: The relative abundance of bacterial and fungal communities at the phylum level. "Others" represents all phyla with relative abundance < 1%. \*\* *p* < 0.01, and \*\*\* *p* < 0.001 (*t*-test); Figure S6: Fungal (a) and bacterial richness (b) in response to plant diversity. \*\*\* *p* < 0.001 (*t*-test); Figure S7: The relationship between fungal (a–d) and bacterial richness (e–h) and soil multifunctionality, at four different thresholds 25%, 50%, 75%, and 90% of maximum. The black lines represent the fitted ordinary least squares (OLS) linear regressions. Red circles represent high plant diversity, while blue ones indicate low plant diversity.

**Author Contributions:** Conceptualization, Z.L. and F.X.; Data curation, Z.L.; Formal analysis, F.X.; Funding acquisition, F.X.; Investigation, Z.L. and F.X.; Methodology, Z.L.; Project administration, F.X.; Resources, Z.L.; Software, Z.L.; Supervision, Z.L. and M.Z.; Validation, Z.L.; Visualization, Z.L.; Writing—original draft, Z.L.; Writing—review and editing, Z.L., X.L. and F.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was financially supported by the National Natural Science Foundation of China (31670524).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data supporting this study's findings are available from the corresponding author upon reasonable request.

**Acknowledgments:** We gratefully acknowledge the staff of the Tongyu Observatory of semi-arid climate and environment for help with the logistics and permission to access the study site. Furthermore, we thank Yu Bi for her field and lab work support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Soil Microbial Diversity and Community Composition in Rice–Fish Co-Culture and Rice Monoculture Farming System**

**Noppol Arunrat 1,\*, Chakriya Sansupa 2, Praeploy Kongsurakan 3, Sukanya Sereenonchai <sup>1</sup> and Ryusuke Hatano <sup>4</sup>**


**Simple Summary:** The integration of fish in rice fields can influence the diversity and structural composition of soil microbial communities. Therefore, soil microorganisms between rice–fish co-culture (RF) and rice monoculture (MC) were compared. The key findings revealed that Actinobacteria, Chloroflexi, Proteobacteria, Acidobacteria, and Planctomycetes were the most dominant taxa across both paddy fields. The most abundant genus in MC belonged to *Anaeromyxobacter*, whereas that in RF was *Bacillus.* Nitrogen fixation, aromatic compound degradation, and hydrocarbon degradation were more abundant in RF. Phosphatase, β-glucosidase, cellulase, and urease enzymes were detected in both paddy fields. However, a 2-year conversion from organic rice to rice–fish co-culture may not be long enough to significantly alter alpha diversity indices.

**Abstract:** Soil microorganisms play an important role in determining nutrient cycling. The integration of fish into rice fields can influence the diversity and structural composition of soil microbial communities. However, regarding the rice–fish co-culture (RF) farming system in Thailand, the study of the diversity and composition of soil microbes is still limited. Here, we aim to compare the microbial diversity, community composition, and functional structure of the bacterial communities between RF and rice monoculture (MC) farming systems and identify the environmental factors shaping bacterial community composition. Bacterial taxonomy was observed using 16s rRNA gene amplicon sequencing, and the functional structures of the bacterial communities were predicted based on their taxonomy and sequences. The results showed that soil organic carbon, total nitrogen (TN), organic matter, available phosphorous, and clay content were significantly higher in RF than in MC. The most dominant taxa across both paddy rice fields belonged to Actinobacteria, Chloroflexi, Proteobacteria, Acidobacteria, and Planctomycetes. The taxa Nitrosporae, Rokubacteria, GAL15, and Elusimicrobia were significantly different between both rice fields. At the genus level, *Bacillus*, *Anaeromyxobacter*, and HSB OF53-F07 were the predominant genera in both rice fields. The most abundant genus in MC was *Anaeromyxobacter*, whereas RF belonged to *Bacillus*. The community composition in MC was positively correlated with magnesium and sand content, while in RF was positively correlated with pH, TN, and clay content. Nitrogen fixation, aromatic compound degradation, and hydrocarbon degradation were more abundant in RF, while cellulolysis, nitrification, ureolysis, and phototrophy functional groups were more abundant in MC. The enzymes involved in paddy soil ecosystems included phosphatase, β-glucosidase, cellulase, and urease. These results provide novel insights into integrated fish in the paddy field as an efficient agricultural development strategy for enhancing soil microorganisms that increase soil fertility.

**Keywords:** microbial diversity; microbial community composition; 16s rRNA gene; rice–fish co-culture; rice monoculture

**Citation:** Arunrat, N.; Sansupa, C.; Kongsurakan, P.; Sereenonchai, S.; Hatano, R. Soil Microbial Diversity and Community Composition in Rice–Fish Co-Culture and Rice Monoculture Farming System. *Biology* **2022**, *11*, 1242. https:// doi.org/10.3390/biology11081242

Academic Editors: Daniel Puppe, Panayiotis Dimitrakopoulos and Baorong Lu

Received: 3 August 2022 Accepted: 18 August 2022 Published: 20 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Microorganisms play important roles in soil and agricultural ecosystems [1,2]. They are responsible for several processes, such as biomass decomposition, nutrient circulation, and soil formation, which subsequently affect plant growth and production [1–3]. In recent years, soil microbial communities have been extensively investigated, as they may reflect soil fertility and ecosystem function [4,5]. Furthermore, the soil microbial community can be used as an indicator to track changes in various land management methods, such as tracking changes in restoration outcomes [6,7] or evaluating agricultural management methods [8].

Integrated rice and fish farming has been practiced in Thailand for more than 200 years by capturing wild seed fish in the rice fields [9]. The rice–fish co-culture (RF) is eulogized for improving ecosystems and alleviating poverty [10] and promoted as increasing biodiversity, reducing fertilizer and pesticide utilization, and contributing to system stability and sustainability [11,12]. Several studies (i.e., [13–15]) reported that RF generated the extra production of aquaculture, which increased farmers' income. Due to fish eating insects, pests, and weeds, the use of pesticides and herbicides can be reduced [16] while organic fertilizers and organic amendments are more applied. These promote the suitable condition for the abundant and diversified population of soil microorganisms, especially bacteria that play a crucial role in soil carbon and nitrogen mineralization. Proteobacteria, Bacteroidetes, Acidobacteria, and Chloroflexi were generally dominant phyla in the paddy soil [17–19], which play important roles in soil nutrient cycles [20,21]. Thus, soil microbial communities can be used as an indicator to explain soil health [22].

To date, the scientific knowledge on soil microbial taxonomic and functional composition, and their interactions with environmental factors of integrating fish in paddy fields remain unclear. Therefore, our study was carried out to fill this gap, aiming to (i) compare microbial diversity and community composition between rice monoculture (MC) and RF fields, (ii) identify the environmental factors shaping the bacterial community composition, and (iii) compare the functional structure of the bacterial communities between both study sites. This study can provide scientific knowledge for the development of a rice–fish co-culture farming system.

#### **2. Materials and Methods**

#### *2.1. Description of Study Sites*

The study sites were located in the Samnak Khun Nen Subdistrict, Dong Charoen District, Phichit Province, Lower North of Thailand. The maximum and minimum temperatures were 32.9 and 23.3 ◦C, respectively, while the average precipitation was 1264.8 mm year−1. A rice–fish co-culture farm (16◦04 04.1" N, 100◦32 31.1" E, Figure 1a) which has been producing organic rice for more than 10 years was selected. The International Federation of Organic Agriculture Movements (IFOAM) Standard was first certified in 2016, while the EU/USDA Organic Standard was approved in 2018. The "Riceberry" rice variety was usually grown once per year from August to December. Since 2019, this farm has been raising fish in the paddy field. The main species of fish were common snakehead (*Channa striata*), walking catfish (*Clarias batrachus* (L.)), and Nile tilapia (*Oreochromis niloticus*). The rice bran, vegetable, and fruit residues and cattle manure were applied in the paddy field as the food for the fish and nutrients for rice. The weeds were removed by hand, while biofermented juice was produced from lemongrass, neem leaves, fruits, and vegetables and then mixed with molasses and animal dung (poultry and cattle) to dispose of the insects. The type of rice–fish co-culture field was the canal refuge (Figure 1b). The transplanting method was used for rice planting, which was performed by hand. One-month-old fish were released into the paddy field 30 days after rice planting. The water level in the field was maintained at around 20–30 cm during rice growing. Rice harvesting was performed by hand, and all rice residues were left in the paddy field. Rice yield was approximately 3.6 tons ha<sup>−</sup>1, while the yield of fish was around 300 kg ha−1.

**Figure 1.** Study areas. (**a**) study sites, (**b**) canal refuge. The aerial image was taken from Google maps on 20 April 2022. The photo was taken on 27 December 2021 by Noppol Arunrat.

For a fair comparison, an adjacent conventional rice farm (16◦04 04.6" N, 100◦32 31.9" E) was selected as the comparison site (Figure 1a). The "RD41" (105 days) or "RD57" (110 days) rice varieties were chosen for planting once a year (August to November). The pregerminated rice seeds were sown by hand (broadcasting method). Then, N, P2O5, and K2O chemical fertilizers were applied, namely 46-0-0 (125 kg ha<sup>−</sup>1) and 16-20-0 (156.3 kg ha−1). Glyphosate (48% *w*/*v* SL) and alachlor (48% *w*/*v* EC) were used to kill the weeds, while acephate (75% S) and chlorpyrifos (40% EC) were applied to eliminate the diseases and insects. A harvesting machine was usually used for rice harvesting, then all rice residues were left in the paddy field. The rice yield was approximately 4.7 ton ha<sup>−</sup>1.

#### *2.2. Soil Sampling and Measurements*

In December 2021, soil samples were collected from the top layer (0–10 cm) of the rice–fish co-culture field and conventional rice field. Five duplicates of soil samples were collected for each field, and 10 soil samples were collected in total. The 50 g soil samples were kept in a cold storage box and brought to the laboratory for the extraction of soil microbial DNA.

Moreover, the soil cores (5.0 cm width × 5.5 cm length) were used to collect the soil samples for soil bulk density measurement and were then measured after 24 h of drying in an oven at 105 ◦C. The 1 kg soil samples were taken for soil physical and chemical properties analysis.

Soil texture was measured using a hydrometer. Electrical conductivity (ECe) was measured using an EC meter in saturation paste extracts (1:5) [23]. Soil pH was determined using a pH meter in a 1:1 soil-to-water mixture [24]. Available phosphorus (Avail. P) was determined following the molybdate blue method (Bray II extraction) [25]. Available potassium (Avail. K), available calcium (Avail. Ca), and available magnesium (Avail. Mg) were measured using atomic absorption spectrometry (NH4OAc extraction pH 7.0) [26]. Total nitrogen (TN) was measured using the micro-Kjeldahl method. Organic carbon (OC) was analyzed using the Walkley and Black [27] method and converted to organic matter (OM) by multiplying by 1.724. The SOC stock was calculated using the following equation:

$$\text{COC stock} = \text{OC} \times BD \times L,\tag{1}$$

where *SOCstock* is the soil organic carbon stock (Mg C ha<sup>−</sup>1), *OC* is organic carbon (%), *BD* is soil bulk density (Mg m<sup>−</sup>3), and *L* is soil depth (cm).

#### *2.3. DNA Extraction and Amplicon Sequencing of 16s rRNA Gene*

DNA was extracted from 0.25 g of soil using DNeasy PowerSoil Pro DNA Kit (Qiagen, Germantown, MD, USA) following the manufacturer's instructions. The extracted DNA was subjected to amplicon library preparation and sequencing. Briefly, PCR amplification, targeting the V3–V4 variable of the 16s rRNA gene, was performed using the universal primers 341F (5 -CCTAYGGGDBGCWSCAG) and 805R (5 -GGA CTACNVGGGTHTC-TAAT) (Klindworth et al., 2013). The amplicons were then sequenced on the Illumina Miseq platform (2 × 250 bp). The amplification and sequencing steps were run by Omics Sciences and Bioinformatics Center (Chulalongkorn University, Bangkok, Thailand).

#### *2.4. Sequencing Analysis and Microbial Taxonomic Identification*

The raw sequence dataset was analyzed with QIIME 2 v. 2021.8 [28]. The 16s rRNA primers were trimmed from forward and reverse reads using cutadapt [29]. The trimmed sequences were quality-filtered (MaxEE = 2; no ambiguous nucleotide) and merged (minimum overlap = 12 nucleotides), and chimeras were removed using the DADA2 plugin [30]. The high-quality sequence was clustered at 97% sequence identity into operational taxonomic units (OTUs) using the VSEARCH plugin [31,32]. Representative sequences of each OTU were taxonomically identified against the Silva v.138 database [33,34]. To eliminate potential sequencing error, rare OTUs (singletons, doubletons, and tripletons) were removed from the dataset. After that, the number of reads that remained in each sample

was randomly subsampled and normalized to the smallest number of reads per sample (24,676 reads/sample), to avoid sequencing depth bias, using rarefy was implemented in QIIME 2. These normalized datasets were used for further analysis.

#### *2.5. Functional Prediction*

Microbial associated functions were predicted using FAPROTAX [35,36] and PI-CRUSt2 [37]. The FAPROTAX predicted the ecologically relevant function of each taxon based on data of the cultured taxa. For example, if all cultured taxa of a bacterial genus were identified as nitrogen-fixing bacteria, all uncultured members of that genus will also be identified as nitrogen-fixing bacteria. On the other hand, PICRUSt2 predicted potential functions based on gene sequences presented in each taxon. In this study, the PICRUSt function was emphasized as enzyme activities that were potentially performed by the detected taxon. These functional analyses were performed following the instructions on the FAPRO-TAX (http://www.loucalab.com/archive/FAPROTAX/lib/php/index.php?section=Home accessed on 27 March 2022) and PICRUSt (https://github.com/picrust/picrust2/wiki accessed on 27 March 2022) webpages.

#### *2.6. Statistical Analysis*

Statistical analysis was performed on PAST [38] and R statistical software [39]. Independent *t*-tests were employed for comparison of soil physicochemical properties between monoculture rice fields and rice–fish co-culture fields. The correlations among the physicochemical variables were observed via Pearson's correlation matrixes. Differences in the relative abundance of microbial taxa detected in the two study sites were indicated using the linear discriminant analysis (LDA) effect size (LEfSe) [40]. Taxa with significant *p*-values (*p* < 0.05) and LDA score ≥ 2 were considered differentially abundant taxa. The LEfSe analysis was performed on an online interface of the Huttenhower lab Galaxy server (http://huttenhower.sph.harvard.edu/galaxy accessed on 27 March 2022). Alpha diversity, including observed OTU richness, Shannon, and Simpson indices were estimated using the diversity indices function in PAST. Differences in the alpha diversity indices between the two study sites were tested via t-test. Beta diversity, representing community composition, was analyzed using non-metric multidimensional scaling (NMDS) ordination based on Bray–Curtis dissimilarity, which was computed using the metaMDS function in the vegan R package. Permutational MANOVA (PERMANOVA) was used to the calculated difference between the two community compositions using the adonis function. The influence of soil properties on soil bacterial community composition was estimated by correlation analysis. The correlations were calculated using the envfit function in the vegan R package, and the *p*-values were corrected by Bonferroni's correction using the p.adjust function in the stat R package. The NMDS ordination with significantly correlated soil parameters was plotted using the ggplot2 R package. Bacterial-associated functions, predicted by both FAPROTAX and PICRUSt2, were visualized as extended bar plots in STAMP software. Statistical differences between each function were tested via *t*-test, and all *p*-values were corrected using Bonferroni's correction. Functional compositions were analyzed following the community composition (beta diversity) analysis as described above.

#### **3. Results**

#### *3.1. Soil Physicochemical Properties in Rice Monoculture and Rice–Fish Co-Culture Fields*

The soil samples, both from the rice monoculture (MC) and the rice–fish co-culture fields (RF), were silty clay. However, significant variances in soil physiochemical properties were found (Table 1). Lower acidity (6.0 ± 0.2, *<sup>p</sup>* < 0.01), bulk density (1.4 ± 0.02 Mg m<sup>−</sup>3, *<sup>p</sup>* < 0.05), ECe (0.4 ± 0.01 dS m<sup>−</sup>1, *<sup>p</sup>* < 0.01), available Ca (2279.0 ± 90.0 mg kg<sup>−</sup>1, *<sup>p</sup>* < 0.01), available Mg (175.1 ± 3.6 mg kg−1, *<sup>p</sup>* < 0.01), and %Sand (10.1 ± 0.8, *<sup>p</sup>* < 0.01) were easily observed in the RF field. Meanwhile, the RF soils also contained significantly higher contents of OM fraction (3.4% ± 0.2, *<sup>p</sup>* < 0.01), SOC (80.9 ± 3.5 Mg C ha−1, *<sup>p</sup>* < 0.01), TN (0.5% ± 0.02, *<sup>p</sup>* < 0.01), available P (20.0 ± 0.9 mg kg−1, *<sup>p</sup>* < 0.01), and %Clay content

(46.3 ± 0.9 mg kg<sup>−</sup>1, *<sup>p</sup>* < 0.01). However, there was no significant difference in available K content or %Silt content between the two sampling sites. Negative correlations between SOC stock and %Silt content were found in MC (r = -0.887, *p* < 0.05). In RF, negative correlations were found between SOC and ECe (r = –0.904, *p* < 0.05) as well as between available K and %Silt content (r = –0.992, *p* < 0.05). Moreover, total nitrogen positively correlated with %Sand content (r = 0.917, *p* < 0.05) (Table 2).

**Table 1.** Comparison of soil physicochemical properties of rice monoculture and rice–fish co-culture fields.


\*, \*\* indicate statistically significant with *p*-value < 0.05 and < 0.01, respectively. NS: No significant, BD = bulk density, OM = organic matter, TN = total nitrogen, ECe = electrical conductivity, CEC = cation exchange capacity, Avail. P = available P, Avail. K = available K, Avail. Ca = available Ca, Avail. Mg = available Mg.

**Table 2.** Pearson's correlation matrixes of soil properties in rice–fish co-culture (n = 5, white area) and rice monoculture fields (n = 5, grey area).


All values in bold print are significant (*p* < 0.05). BD = bulk density, OM = organic matter, TN = total nitrogen, ECe = electrical conductivity, CEC = cation exchange capacity, P = available P, K = available K, Ca = available Ca, Mg = available Mg.

#### *3.2. General Overview of the Sequencing Analysis*

A total of 337,778 high quality and abundance readings, representing 4597 OTUs, were derived from this study. After normalization, 4582 OTUs remained. Rarefaction curves of the detected OTUs derived from both MC and RF samples were flattened at the analysis

sequencing depth (24,676 reads/sample), indicating that the detected OTUs were sufficient to represent the microbial community in each sample (Figure 2).

**Figure 2.** Rarefaction curve of observed microbial OTUs detected in monoculture (MC) and rice–fish (RC) fields.

#### *3.3. Taxonomic Distribution and Differential Abundance of Soil Bacteria in Rice Monoculture and Rice–Fish Co-Culture Fields*

According to microbial analysis based on the 16s rRNA gene, microbial taxa detected in this study were classified into 47 phyla, 128 classes, 235 orders, 318 families, 479 genera and 4582 OTUs. Taxonomic distribution of the abundant bacteria (total relative abundance > 0.1%) is presented in Figure 2. The most dominant taxa across all samples belonged to Actinobacteria (MC = 22.78%, RF= 24.17%, on average), followed by Chloroflexi (MC = 18.44%, RF= 17.77%), Proteobacteria (MC = 18.25%, RF= 17.28%), Acidobacteria (MC = 11.16%, RF= 11.88%), and Planctomycetes (MC = 10.44%, RF= 8.76%) (Figure 3a). Taxa that were significantly different between the two sites were Nitrosporae, Rokubacteria, GAL15, and Elusimicrobia. The latter was enriched in RF, whereas the other three were enriched in MC (Figure 3b).

The difference in the microbial community was more noticeable at a deeper taxonomic level, as shown in Figure 4. More than 50% of all detected OTUs were found uniquely in one of the two study sites (Figure 4a). According to the LEfSe analysis at the phylum to genus level, a total of 135 differentially abundant taxa (*p* < 0.05; LDA score > 2) were detected, 111 of which were more abundant in MC than RF and 24 of which were more abundant in RF than MC (Figure 4b, Figure S1). At the genus level, *Bacillus*, *Anaeromyxobacter*, and HSB OF53-F07 were the predominant genera in both rice fields. The most abundant genus in MC was *Anaeromyxobacter*, whereas that in RF was *Bacillus* (Figure 4c). Seven out of the top 30 most prevalent genera, for example, *Bradyrhizobium*, *Bryobacter*, *Conexibacter*, *Nocadiodides,* and *Solirubrobacter,* were significantly more abundant in MC than in RF (Figure 4c). However, some low-abundance genera, such as *Chlorobaculum*, *Niastella,* and *Vicinamibacter*, were more abundant in RF than in MC (Figure S1).

**Figure 3.** Taxonomic distribution at phylum level. (**a**) The relative abundance of microbial phyla in each sample. (**b**) LEfSe analysis of differential taxonomy abundance at phylum level between the two study sites. The orange and blue horizontal bars represent the taxa enriched in monoculture and rice–fish fields, respectively.

#### *3.4. Richness, Diversity, Community Composition, and Their Correlation to Soil Properties*

Alpha diversity, reflecting richness and diversity, was indicated by observed OTU richness, Shannon, and Simpson indices. As shown in Figure 5, all alpha diversity indices were slightly higher in MC than in RF, but no significant difference was found between the two sites (*p* > 0.05, Figure 5a–c). On the contrary, beta diversity, presented by NMDS ordination with Bray–Curtis distance, showed a separated community between MC and RF (Figure 5d). This indicated that the community composition of bacteria in MC was different from that in RF. The difference was confirmed by PERMANOVA test (*F* = 0.251, *p* = 0.008).

The correlations between soil properties and bacterial community composition are shown in Table 3; 5 out of 12 measured parameters were significantly correlated with bacterial community composition in both study sites. Whilst the community composition in MC was positively correlated with Mg and sand, the same in RF was positively correlated with pH, TN, and clay (Figure 5). Mg was the most correlated factor (r<sup>2</sup> = 0.880), following closely by Sand (r2 = 0.866), pH (r2 = 0.857) and TN (r2 = 0.834) (Table 3).

**Figure 4.** Taxonomic differences of bacteria in monoculture and rice–fish co-culture fields. (**a**) The Venn diagram shows the number of OTUs found in monoculture (red) and rice–fish (green) fields. (**b**) The cladogram shows differential abundance taxa among the two rice fields. More information on the differential taxa is provided in Supplementary Figures S1 and S2. (**c**) Heatmap shows the relative abundance of the thirty most abundant microbial genera detected in each sample. The bar plot presents the mean relative abundance of the microbial genera detected in monoculture (orange) and rice–fish fields (green). Genus names with an asterisk are statistically different between the two sites (*p* < 0.05). MC: Monoculture, RF: rice–fish.


**Table 3.** Pearson's correlation between microbial community composition and soil properties.

\* Indicate significant parameters (*p* < 0.05). BD = bulk density, OM = organic matter, TN = total nitrogen, ECe = electrical conductivity, CEC = cation exchange capacity, Avail. P = available P, Avail. K = available K, Avail. Ca = available Ca, Avail. Mg = available M.

**Figure 5.** Bacterial diversity and community composition. Bar plot shows (**a**) OTU richness, (**b**) Shannon's, and (**c**) Simpson's indices measured in monoculture and rice–fish fields. (**d**) NMDS ordination, based on the Bray–Curtis dissimilarity measure, presents the community composition of soil microbes in the two study sites. MC: Monoculture, RF: rice–fish.

#### *3.5. Predictive Function*

Totals of 807 (17.61%) and 4532 (98.90%) OTUs were assigned to at least one function of the 63 ecologically relevant functions and 2238 enzymes, respectively. For FAPROTAX analysis, Cellulolysis, nitrification, ureolysis, phototrophy, nitrogen fixation, aromatic compound degradation, and hydrocarbon degradation, were found among the top 20 most abundant functions (Figure 6a). While the first four functional groups were more abundant in MC, the last three groups were more abundant in RF. However, no significant change was detected in any of those functions (*p* > 0.05), which is consistent with the results from the PICRUSt analysis. Enzymes involved in soil systems, such as phosphatase, βglucosidase, cellulase, and urease, were presented in this study. As shown in Figure 6b, no significant difference was found between the enzymes detected in MC and RF (*p* > 0.05). The functional potential structures, created from all detected ecologically relevant functions and enzymes, were shown in Figures 6c and 6d, respectively. Based on NMDS ordination and PERMANOVA analysis, no significant difference was detected (*p* > 0.05) between the functional structures in MC and RF.

**Figure 6.** Predictive functions. Extended bar plots show mean and difference in mean proportions of microbial function predicted by (**a**) FAPROTAX and (**b**) PICRUST2. Functions overrepresented in the monoculture (red) have a positive difference in the proportion, whereas functions overrepresented in the rice–fish (green) have a negative difference. NMDS ordinations based on the Bray–Curtis dissimilarity measure show the composition of bacterial functions, predicted by (**c**) FAPROTAX and (**d**) PICRUST2, in monoculture (red) and rice–fish field (green). MC: Monoculture, RF: rice–fish.

#### **4. Discussion**

#### *4.1. RF Farming System Can Increase Soil Nutrients*

Our results (Table 1) reveal that the rice–fish co-culture system significantly increased SOC stock and TN similarly to other studies [41–44], which found that the rice–fish coculture system can potentially increase the content of organic carbon and nitrogen in the soil through its mineralization of organic matter. This is due to the uneaten and excess feed as well as excreta produced during fish growth increasing SOC content and TN, which is consistent with the findings in the rice–crayfish–turtle co-culture by Li et al. [45] and rice– crayfish farming system by Si et al. [46]. They also revealed that crayfish and turtle feeding, molting, and excretion could increase the SOC and TN, while these aquatic animals help increase soil permeability by penetrating the soil surface in the paddy field. The supplies of some available elements (i.e., Ca and Mg) in the RF farming system were significantly lower than in MC (Table 1), which were possibly shaped by its divalent charge fraction on the ECe [47], following the lower ECe value in RF compared to MC (Table 1). The rice–fish co-culture system, however, shows a higher content of available P (Table 1), while in rice monoculture, it generally decreased due to the long-term cultivation. Slowly cycling P and labile P in soil increasing with time was also previously reported in the rice-frog-fish culture soil [48]. Higher clay content and organic matter were found in RF compared with MC due to a higher amount of organic amendments being added. Hassink [49] reported that clay particles could absorb organic matter and protect it from microbial decomposition.

#### *4.2. Microbial Diversity and Community Composition under Rice–Fish Co-Culture and Rice Monoculture*

To our knowledge, this study was among a few studies investigating the differentiation of both community and functional potential structures of soil bacteria in rice monoculture fields, compared with rice–fish co-culture fields. Here, we found that although the alpha diversity of bacteria did not differ significantly between MC and RF, the community composition did. Even though the alpha diversity indices, which present observed richness, Shannon's, and Simpson's indices, differed between the two sites, the differences were minor and did not reveal a significant difference. This result is consistent with previous studies on other rice co-culture fields, including rice–frog (Yi et al., 2019), rice–crab [19], and rice–fish fields [7,45]. Whilst Zhao et al. [7] reported that microbial diversity in rice–fish co-culture fields was similar to those in monoculture rice fields in the first year of cultivation and substantially changed after 5 years, Li et al. [45] found that it took up to 12 years to see the differences in rice–fish–turtle fields. Here, we showed that changes were still not found after three years of rice–fish cultivation. However, when it comes to community composition, all studies [7,45,50], including this one, found significant differences between co-culture and monoculture rice fields. A total of 135 differentially abundant taxa were found in both the MC and RF sites, four of which were among the abundant phyla and seven of which were among the top 30 most abundant genera, respectively. These results implied that three years of rice–fish co-culture agriculture, which was implemented after long-term monoculture rice cultivation, had a significant impact on microbial community composition, while diversity remained unchanged.

As shown in Figure 3a, the most dominant taxa in paddy fields (both RF and MC farming systems) belonged to Actinobacteria, Chloroflexi, Proteobacteria, Acidobacteria, and Planctomycetes. Meanwhile, all alpha diversity indices showed no significant difference between the two sites. This is consistent with the study of Jiang et al. [19], who found that the most abundant phyla in paddy soil and ditch sediment under the rice– crab co-culture system were Proteobacteria, Bacteroidetes, and Chloroflexi, while alpha diversity of bacterial diversity in paddy soil and ditch sediment was similar. The phylum Actinobacteria, which contributes to the turnover of the complex biopolymers and plant residue decomposition [51,52] by producing various carbon cycling enzymes [53], was the most dominant taxa in both RF and MC sites. The phylum Chloroflexi, which involves nitrification was the second-most dominant that was detected in RF and MC sites [54]. The phylum Proteobacteria is usually classified as "copiotrophs" (R-strategists), which indicates that its members are more abundant and have high growth rates under nutrientrich conditions [55]. It plays a key role in OM decomposition, produces many types of glycosyl hydrolases, and is involved in nitrogen fixation, which promotes plant growth [1]. The phylum Acidobacteria is involved in the carbon cycle of humus decomposition [56] and adjusts soil pH [57]. Members of the bacterial phylum Planctomycetes can act as slow-acting degraders of various biopolymers, cellulose, and chitin. They can degrade exopolysaccharides produced by other bacteria [58]. García-Orenes et al. [59] reported that some members of the phylum Planctomycete, such as *Blastopirellula,* are good indicators of soil fertility because they respond to the applications of manure or fertilizers faster than they react to soil physicochemical properties.

Nitrosporae, Rokubacteria, GAL15, and Elusimicrobia taxa were significantly different between the RF and MC fields (Figure 3b). This result indicated that the structure of the soil bacterial community was significantly changed after integrating fish into the paddy field. The phylum Nitrospirae plays a crucial role of decomposing soil mineral nitrogen and improving nutrient availability [60], and its functions can enhance crop productivity [61]. At the genus level, *Bacillus*, *Anaeromyxobacter*, and HSB OF53-F07 were found in both rice fields. The most abundant genus in MC was *Anaeromyxobacter*, whereas *Bacillus* was the most dominant in RF (Figure 4c). *Bacillus* plays multiple functions in the soil ecosystem for nutrient cycling, which is involved in nitrogen fixation, phosphorus nutrition, and potassium solubilization that promotes plant growth [62]. Similarly, *Anaeromyxobacter* can

fix and assimilate N2 gas via its nitrogenase [63]. This is why soil nutrients (TN, Avail. P, and Avail. K) were higher in RF than in MC soil (Table 1).

Soil physicochemical properties have a great influence on bacterial community composition [7,50,64,65]. Our findings indicated that the different soil management practices of RF and MC farming systems cause a significant difference in soil physicochemical properties, resulting in the different composition of soil bacterial communities and their metabolic functions (Figure 5 and Table 3). This is in line with the study of Viruel et al. [66], who reported that changes in SOC stock, TN, and pH cause the changes in soil bacterial communities and metabolic functions in farming systems of the semi-arid Chaco region, Argentina. Hartmann et al. [67] found that soil pH, SOC, and TN were the main predictors of bacterial community structure in long-term organic and conventional farming. As presented in Figure 5 and Table 3, we found that Mg, TN, pH, and soil texture (percent of sand and clay) were the main factors determining the community composition of bacteria in MC and RF. TN and pH have previously been identified as the key factors that influence the microbial community in rice co-culture fields [7,50]. While Mg was rarely measured in earlier works, this study found that it was one of the most important elements affecting community composition in rice fields. Hou et al. [68] found a strong effect of Mg on bacterial nitrification, which is in line with the study of Zhang et al. [69], who reported that the appropriate concentrations of Mg2+ could promote nitrification activity in the soil.

In addition, this study also showed the functional potential structure of the bacterial community, based on FAPROTAX and PICRUSt2, in MC and RF (Figure 6). Interestingly, despite differences in community composition between MC and RF, no difference in functional potential structure was found. Furthermore, there was no difference in the abundance of the individual function predicted by either tool. Recently, Chen et al. [70] demonstrated on global soil metagenomic data that microbial functional structure could remain stable while taxonomic diversity and composition changed, which is consistent with our findings. We reveal that, whilst bacterial community composition changed during the transition from monoculture to rice–fish co-culture, the functional structure could be maintained, and this was the first study to report such results on monoculture and rice–fish fields. However, it should be noted that the limitation of the prediction tools was that only taxa or genes from the databases were functionally assigned [35–37,71]. Regardless, using both tools to estimate the availability and abundance of genes or functions within the community could provide more insight into information on the complex community in soil.

#### **5. Conclusions**

This study shows that soil physicochemical properties (SOC, TN, OM, Avail. P, and clay content) were significantly higher in RF than in MC sites. The most dominant taxa across both paddy rice-farming systems belonged to Actinobacteria, Chloroflexi, Proteobacteria, Acidobacteria, and Planctomycetes. The taxa Nitrosporae, Rokubacteria, GAL15, and Elusimicrobia were significantly different between the RF and MC sites. At the genus level, *Bacillus*, *Anaeromyxobacter*, and HSB OF53-F07 were the predominant genera in both rice fields. The most abundant genus in MC was *Anaeromyxobacter*, that in RF was *Bacillus.* All alpha diversity indices between the RF and MC sites were not significantly different. The community composition in MC was positively correlated with Mg and sand, while in RF was positively correlated with pH, TN, and clay. Nitrogen fixation, aromatic compound degradation, and hydrocarbon degradation were more abundant in RF, while cellulolysis, nitrification, ureolysis, and phototrophy functional groups were more abundant in MC. The enzymes involved in paddy soil ecosystems included phosphatase, β-glucosidase, cellulase, and urease, but no significant difference was detected between MC and RF fields. However, more observation fields and long-term monitoring are necessary to conclusively confirm our findings in this study.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology11081242/s1, Figure S1: LEfSE analysis from phylum to genus level. The cladogram shows differential abundance in taxa (*p* > 0.05, LDA score ≥ 2) between monoculture rice fields and rice–fish co-culture fields. Taxa in red are significantly enriched in the monoculture, whereas taxa in green are significantly enriched in the rice–fish field. Figure S2: Bar plot shows differential abundance taxa with effect size (LDA score). Taxa in red are significantly enriched in the monoculture, whereas taxa in green are significantly enriched in the rice–fish field.

**Author Contributions:** Conceptualization, N.A.; methodology, N.A., C.S. and P.K.; investigation, N.A. and S.S.; writing—original draft preparation, N.A., C.S., P.K., S.S. and R.H.; writing—review and editing, N.A., C.S., P.K., S.S. and R.H.; supervision, R.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research project is supported by Mahidol University (contract number: A14/2564).

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of Institute for Population and Social Research, Mahidol University (IPSR-IRB) (COA. No. 2021/01-001).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Raw sequence data generated for this study are available in the Sequence Read Archives (SRA) of the National Center for Biotechnology Information (NCBI) under BioProject accession number: PRJNA819169.

**Acknowledgments:** The authors extend their appreciation to Mahidol University (contract number: A14/2564).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

