*Article* **Geology Can Drive the Diversity–Ecosystem Functioning Relationship in River Benthic Diatoms by Selecting for Species Functional Traits**

**Evangelia Smeti 1,\*, George Tsirtsis <sup>2</sup> and Nikolaos Theodor Skoulikidis <sup>1</sup>**


**Simple Summary:** The way that diversity affects ecosystem functioning is of great importance, as it helps us understand the health state of an ecosystem. Primary producers contribute to ecosystem functioning through biomass production, which is considered to be a proxy of ecosystem functioning. In rivers, the primary producers of the biofilm are diatoms, unicellular algae with cell walls of silica. In this study, we tested the way diatom species affect biomass production across nine rivers in Greece. Nutrient concentrations that drive primary production are linked to river geology. We found that the geological substrate of a river could be responsible for the diversity–biomass relationship: in rivers with a siliceous substrate, more diatom species increased biomass, whereas in rivers with a calcareous substrate, a change in diatom species number did not change biomass. By using model simulations, we found that this difference could be attributed to the different stages of the biofilm in time. Our results show the importance of different factors that affect diatom species, their functional traits and biomass production and what we should consider when testing for ecosystem functioning.

**Abstract:** The biodiversity–ecosystem functioning (BEF) relationship has been studied extensively for the past 30 years, mainly in terrestrial plant ecosystems using experimental approaches. Field studies in aquatic systems are scarce, and considering primary producers, they mainly focus on phytoplankton assemblages, whereas benthic diatoms in rivers are considerably understudied in this regard. We performed a field study across nine rivers in Greece, and we coupled the observed field results with model simulations. We tested the hypothesis that the diversity–biomass (as a surrogate of ecosystem functioning) relationship in benthic diatoms would be affected by abiotic factors and would be time-dependent due to the highly dynamic nature of rivers. Indeed, geology played an important role in the form of the BEF relationship that was positive in siliceous and absent in calcareous substrates. Geology was responsible for nutrient concentrations, which, in turn, were responsible for the dominance of specific functional traits. Furthermore, model simulations showed the time dependence of the BEF form, as less mature assemblages tend to present a positive BEF. This was the first large-scale field study on the BEF relationship of benthic diatom assemblages, offering useful insights into the function and diversity of these overlooked ecosystems and assemblages.

**Keywords:** biofilm; diatoms; rivers; Greece; model simulations

### **1. Introduction**

Ecosystem functioning comprises multiple processes that account for ecosystem health and sustain ecosystem services. For the past 30 years, research has been focusing on proving the pivotal role of diversity in driving ecosystem functioning [1]. In particular, the study of the form of the biodiversity–ecosystem functioning (BEF) relationship is considered important in light of global change and species extinctions [2]. It could be further used as a proxy of ecosystem health, resilience and species interactions, providing great insight

**Citation:** Smeti, E.; Tsirtsis, G.; Skoulikidis, N.T. Geology Can Drive the Diversity–Ecosystem Functioning Relationship in River Benthic Diatoms by Selecting for Species Functional Traits. *Biology* **2023**, *12*, 81. https:// doi.org/10.3390/biology12010081

Academic Editor: José Carlos Hernández

Received: 13 December 2022 Revised: 29 December 2022 Accepted: 30 December 2022 Published: 4 January 2023

**Copyright:** © 2023 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/).

into an ecosystem's need for conservation [3,4]. BEF studies, especially early ones, have focused on experimental work, mainly on terrestrial plants, whereas aquatic environments and especially freshwater, remain understudied [5]. Field studies are rare, and regarding microalgae, they have focused on phytoplankton in Scandinavian and USA lakes [6,7] and the Baltic Sea [6,8]. A few studies on biofilm are limited to estuaries [9,10], but studies on the productivity of the river biofilm, especially on benthic diatoms, are almost missing (see [11]).

Diatoms, a major component of phytobenthos in rivers and the most diverse group of protists, are unicellular algae with silica cell walls, responsible for 20% of O2 production, and are important indicators of water quality [12]. As primary producers, their growth depends on nutrient concentrations and light, contributing immensely to primary biofilm productivity, an important ecosystem function. The importance of benthic diatoms on biofilm biomass production, along with their high diversity, renders them an ideal group of organisms for studying the BEF relationship in river biofilms. Furthermore, in recent years, functional traits related to cell size, to adherence to substrates and life forms are increasingly used in describing benthic diatom assemblages [13,14]. Despite the growing evidence that functional richness could be more important in driving ecosystem functions than taxonomic richness [15], functional diversity metrics are not widely used in BEF studies.

Although initial research focused on positive BEF relationships (i.e., an increase in ecosystem functions with increasing diversity), recent research and meta-analysis suggest different relationships (e.g., negative or hump-shaped relationships), depending on the type and duration of the study (i.e., observational field or experiment), the ecosystem type and the taxa studied [5]. Furthermore, the form of the BEF relationship can be affected by multiple diversity levels in space and time, such as the regional and local species pools or the initial and realized diversity [16]. Although they have the advantage of the large scale and the natural world, observational field studies can give ambiguous results. Natural systems are very complex and dynamic, and patterns could be masked by different diversity scales that are part of different assembly phases. On the other hand, modeling studies could provide mechanistic interpretations of the form of the BEF relationship [17,18]. For example, modeling studies on phytoplankton have suggested that when species utilize all the available resource space (e.g., when the system is at a mature steady-state), there is no BEF relationship, whereas when a part of the resource space stays unutilized (e.g., after a species extinction), then a strong positive BEF is apparent [17]. Therefore, coupling field observations with numerical modeling could give better insights into the drivers of the BEF relationship.

The aim of this study was to test the form of the BEF relationship in benthic diatoms in rivers and the drivers of this relationship. Toward this aim, we collected samples along nine Greek rivers, varying in their geographical location, geology, drainage area and nutrient concentrations. As rivers are highly dynamic ecosystems and Greece is a country with diverse landscapes and geology, we hypothesized that a general BEF relationship would be hard to observe and that it would be driven by additional, possibly abiotic, factors. We further investigated species traits (size, attachment to the substrate) that can be responsible for the observed relationship. To further understand species coexistence and the consecutive biomass production during a succession of the biofilm in time, we ran model simulations and checked the BEF relationship at different maturity levels of the biofilm. We hypothesized that a less mature biofilm (at the beginning of time succession) would present a positive BEF, as resources are still unutilized, whereas a mature biofilm, where all available resources are used, would not present a BEF. To the best of our knowledge, this is the first large-scale observational study of the BEF relationship in benthic river diatoms, and although incomplete, it can give important insights into the function and diversity of these overlooked systems.

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

This study combined field observations with numerical modeling. Field sampling was conducted in river biofilm, comprising samples for both microscopic observation and biomass. Water physico-chemical parameters and nutrient concentrations were measured at each site. Diversity (i.e., species richness and evenness) of diatoms in the biofilm was defined using microscopic counts of species abundances, and functional traits were assigned to species to account for functional diversity. Biomass was measured in the lab as chlorophyll a concentration and was used as a proxy of ecosystem functioning. The shape of the BEF relationship was tested at different spatial scales, and nutrient concentrations and functional traits were investigated as possible drivers of the observed shape. Model simulations, using a well-known numerical model on species competition for available resources, were run, and the BEF relationship was observed at different time points to test for the dependence of the relationship to the maturity of the assemblage.

#### *2.1. Field Sampling*

Nine Greek rivers (Nestos, Lissos, Fonias, Spercheios, Mornos, Alfeios, Arkadikos, Neda and Evrotas) were sampled in the summer of 2020, at a low flow period, when no major disturbances would cause shifts in the assemblages and their biomass (Figure A1, Table A1). These rivers were selected based on accessibility and appropriate sampling substrate (stones) as well as due to their differences in terms of size, geology and environmental conditions. In each river, five sampling sites were sampled from upstream to downstream, apart from Arkadikos and Lissos, where only four samples were taken. In order to ensure replication, in each site, three spots were sampled, each comprising three stones. From each stone, two surfaces of the defined area were scraped, the first used for chlorophyll analysis (immediately put in a dark bag and frozen) and the other for species identification and counting (preserved with 70% ethanol). This ensured the direct comparison between species diversity and biomass production. At each site, physico-chemical parameters (Temperature, DO, pH, Conductivity, Turbidity) were also measured in situ using a Portable multiparameter Aquaprobe, and water samples were collected for the determination of nutrients (NO3, NO2, NH4, TN, PO4, TP and SiO2).

#### *2.2. Analysis of Samples*

In the laboratory, after filtration through 0.45 μm pore size membrane filters, nutrients were determined by a Skalar San++ Continuous Flow Analyzer [19]. For the determination of chlorophyll, the trichromatic equations were applied [20], where all three main chlorophylls were measured (Chl-a, Chl-b, Chl-c), and their concentrations were in mg/cm2. Chl-a is a measure of the whole phytobenthos biomass production, whereas Chl-c is more indicative of the biomass produced by benthic diatoms.

Diatom species samples were treated with hot hydrogen peroxide to remove organic matter and obtain clean frustules, to be used for diatom species identification [21]. Clean frustules were mounted with Naphrax®, identified to species level with a light microscope (Nikon Eclipse Ci-L, Nikon Microscope Solutions, Europe) at 1000× magnification and counted until no more new species were detected in each sample. As the surface scraped out of the stone was defined, and the volume at each step of the procedure was also measured, the counting reflected the absolute abundance of cells per cm2. For the taxonomy, the work of [22] was mainly used.

#### *2.3. Data Analysis*

In order to ensure that the sampling effort was adequate for all rivers examined, species accumulation curves (SACs) were constructed, showing that, indeed, most species were observed under the specific sampling and analysis procedures (Figure A2). Furthermore, to check that species richness counts were not biased due to macroecological patterns and large differences in the drainage areas, species–area relationships (SAR) for the selected rivers were performed, demonstrating the absence of a relationship between the area and the observed species richness (Figure A3). Another potentially confounding factor in natural systems is pollution. In the present study, pollution levels slightly differed, even between sites of the same river, based on a biological quality diatom index, but quality classes did not play an important role in the BEF relationship (interaction term *p*-value = 0.07).

Taxonomic diversity was calculated using both species richness (S) and evenness (J), to account for the abundance distribution of individuals among species (i.e., assemblage structure). For calculating functional diversity, the functional richness index was used, defined as the total branch length of a functional dendrogram based on species' functional traits [23]. Functional traits used were cell size (L/W ratio, biovolume), substrate adherence (high profile, low profile, motile and planktonic guilds), life forms (colonial, singular) and nitrogen fixation [14]. For the calculation of biovolume, equations of geometric shapes were used [24], and dimensions that could not be measured in our samples (e.g., cell height) were defined based on the literature [14]. The total biovolume of each sample was divided by the total abundance of the sample for the calculation of the average cell size of each assemblage, aiming to compare cell size between different groups of rivers [25].

Biomass metrics tested were Chl-a, expressing biomass production of the entire biofilm, Chl-c and Total biovolume, linked to benthic diatoms. Chl-a and Chl-c were highly correlated (Spearman r = 0.86, *p*-value < 0.001) and presented the same trends. Therefore, only Chl-a was used as a surrogate of biomass production. In order to show more clearly linear trends, Chl-a concentrations were ln-transformed.

The form of the relationship between the different diversity metrics and Chl-a was determined for the whole dataset, searching for a general pattern in the examined Greek rivers, as well as for each river separately to test for possible differences in the BEF relationship between rivers. Rivers were further grouped in two previously defined hydrochemical zones in Greece [26] with distinct silicate and phosphorus concentration ranges, affected by geology; in zone A, siliceous substrates are more prominent, and silicate and phosphorus concentrations in water are higher, whereas, in zone B, calcareous substrates dominate and silicate and phosphorus concentrations in water are lower, the latter due to adsorption on carbonate-rich particles and sediments [27,28]. Substrate geology (i.e., siliceous vs. calcareous) is known to select for species diatom species [22]. Therefore, as phosphorus and silica are important nutrients for diatom growth, their different concentrations in these two zones could affect assemblage characteristics and, thus, the corresponding BEF relationships. For the rest of the manuscript, when we refer to substrates (siliceous or calcareous), we refer to the geologic substrate of a river basin.

For testing the significance of the BEF relationship when used in different groups in the dataset, generalized linear mixed-effects models were used, with the river as a random factor. Data analyses and illustrations were performed in R (v. 4.0.3) [29], using packages vegan v. 2.5-7 [30], BAT v. 2.7.1 [31], lme4 v. 1.1-28 [32], ggplot2 v. 3.3.5 [33] and plotly 4.10.1 [34].

#### *2.4. Model Simulations*

Model simulations were performed in an effort to understand the importance of temporal succession and the maturity of the biofilm on the BEF relationship. Applied models were based on well-known models for phytoplankton competition for resources, assuming a continuous inflow of nutrients [35]. This model describes the population dynamics of 400 diatom species (Ni) competing for two nutrients (*Rj*), namely nitrogen and phosphorus. The initial species number (*n* = 400) is based on the total diatom species observed in all field samples.

$$\frac{dN\_i}{dt} = N\_i \left( \min\left(\frac{\mu\_{\max\_i i} \times R\_1}{K\_{1i} + R\_1}, \frac{\mu\_{\max\_i i} \times R\_2}{K\_{2i} + R\_2}\right) - m\_i \right), \quad i = 1 - n \tag{1}$$

$$\frac{dR\_i}{dt} = D\left(S\_{\bar{j}} - R\_{\bar{j}}\right) - \sum\_{i=1}^{n} c\_{\bar{j}i} \times \min\left(\frac{\mu\_{\max\_i i} \times R\_1}{K\_{1i} + R\_1}, \frac{\mu\_{\max\_i i} \times R\_2}{K\_{2i} + R\_2}\right) \\ \text{N}\_i, \quad j = 1, 2 \tag{2}$$

*Ni* is biomass of species i, and *Rj* is the concentration of nutrient j; *μmaxi* is the specific maximum growth rate of species i, and Kji is the half-saturation constant of resource j for species i, based on the Monod model of growth limitation; mi is the mortality induced by flushing, and it was calculated as the flushing rate (D) divided by the maximum growth rate of each species in the model; *D* is the nutrient flushing rate; *Sj* is the input nutrient j concentration; and *cji* is the intracellular content of nutrient j in species i. In diatoms, the maximum growth rate is linked to species size, with smaller species presenting a higher growth rate [36]. As larger species in a biofilm tend to be more affected by flushing than smaller species that tend to adhere to the substrate stronger, we assumed that they are more affected by flushing, which increases larger species mortality.

The two nutrients used in the model are phosphorus and nitrogen, as they are both essential nutrients for growth. Based on field observations, phosphorus was mainly the limiting nutrient, whereas nitrogen limitation was also observed in some cases. The two nutrients in the model are added synchronously and in a continuous manner during the simulations at concentrations following the Redfield ratio. This synchronous and continuous flow simulates an ideal river environment, from upstream (nutrients entering the system) to downstream (nutrients flushing). Following the N:P:Si ratio in field observations, Si was never found to be limiting; therefore, even though an important nutrient for diatom growth, it was not considered in model simulations.

Life history traits were assigned to species based on the literature values and on species functional traits that we observed in the field samples. The three main life history traits we focused on were the specific maximum growth rate (μmax), the competitive ability for Phosphorus (KP) and the competitive ability for Nitrogen (KN). Based on field data, smaller species tended to be at low nutrient concentrations; therefore, we assigned three groups of species, with each group being superior for two life history traits: one group consisted of fast-growing species with the increased competitive ability for phosphorus but not for nitrogen (high μmax and low KP but high KN), one group consisted of fast-growing species with the increased competitive ability for nitrogen (high μmax and low KN but high KP) and one group consisted of slow-growing species with the increased competitive ability both for phosphorus and nitrogen (low μmax and low KP and KN). Keeping a trade-off was important as the presence of a "superspecies", superior for all traits, would exclude all other species, and thus, species richness in an assemblage would be extremely low. When assigning traits to virtual species, we made sure that there was a trade-off between R\*P and R\*N, with a level of complementarity equal to 0.49 [37]. R\* is the minimum concentration of a resource at which a species could keep its population stable, and it is a summary value of both growth rate and Kj. Life history traits were assigned to species using R (v. 4.0.3).

The mathematical equations were solved numerically using a specially developed Fortran code following [37] and adapted to meet the characteristics of the studied systems. The BEF relationship was tested at each time step using the species richness and evenness against the log-transformed abundance of the 100 replicates. For each replicate, the initial biomass of each species and the total initial abundance varied randomly. The model parameters values, ranges and initial conditions are detailed in Table A2.

#### **3. Results**

#### *3.1. Field Observations*

The general BEF relationship (when all samples were pooled together) when using species richness (S) as the diversity predictor of biomass was positive, albeit rather weak (*p*-value < 0.01, Figure 1a). A seemingly similar but not significant trend was apparent when functional richness was used as a biomass predictor (*p*-value = 0.235, Figure 1c). The lack of a significant relationship was also present when evenness (J) was used as a diversity predictor of biomass (*p*-value = 0.676, Figure 1b).

**Figure 1.** Diversity (expressed as (**a**,**d**) species richness S, (**b**,**e**) evenness J and (**c**,**f**) functional richness) and ecosystem functioning (presented as ln(Chl-a) on y-axis) relationships in each river on the whole dataset (**a**–**c**) and at different geological substrate levels (**d**–**f**).

When each river was considered separately, the relationship between species richness and biomass production was variable between the different rivers sampled (Figure A4a). Indeed, Chl-a was best explained when the interaction between species richness and the river was also considered (adjusted R2 = 0.55, *p*-value < 0.001). This variation is also apparent when considering other diversity metrics (evenness J and functional richness, Figure A4b,c). Variability among rivers was also evident in environmental conditions, as depicted in the physico-chemical parameters and nutrient concentrations measured (Figure A5).

A strong interaction effect is apparent when testing for the substrate geology (interaction term *p*-value < 0.001). In siliceous substrate, an increase in species richness resulted in an increase in biomass production (positive BEF-slope = 0.097, *p*-value < 0.001), whereas, in the calcareous substrate, an increase in species richness did not have any effect on biomass production (no BEF relationship-slope = −0.0014, *p*-value = 0.9—Figure 1d). Functional richness presented the same trend (interaction term *p*-value < 0.05, Figure 1f), but evenness (J) had no effect on predicting biomass ((interaction term *p*-value = 0.204, Figure 1e). There was no significant difference between species richness or biomass for the two groups of rivers.

In rivers with a siliceous substrate, all tested nutrients (TinN (i.e., sum of NO2, NO3, NH4), PO4, SiO2) presented higher concentrations than in rivers with a calcareous substrate (*p*-value < 0.05–Figure 2a–c). Regarding species traits, rivers in siliceous substrates have diatom assemblages comprised of bigger and motile species, whereas rivers in calcareous substrates have diatom assemblages comprised of smaller, low-profile species (*p* < 0.05–Figure 2d–f). Overall, motile species tended to increase with increased phosphorus concentrations, whereas low-profile species tended to decrease with increased phosphorus concentrations (Figure 3a,b). The other guilds (high-profile and planktonic species) did not present any consistent relationship between geology or nutrient concentrations. Furthermore, a higher relative abundance of low-profile species resulted in high dominance assemblages (Figure 3d), whereas higher evenness was observed when more motile species were present (Figure 3c).

**Figure 2.** Difference between nutrient concentrations: (**a**) total inorganic nitrogen; (**b**) orthophosphates; (**c**) silicates and species functional traits; (**d**) average size in the assemblage; (**e**) relative abundance of motile species; (**f**) relative abundance of low-profile species between calcareous and siliceous substrates. Y-axis in (**a**–**d**) is ln-transformed.

**Figure 3.** Relationship between PO4 concentration and diatom guilds: (**a**) relative abundance of motile species; (**b**) relative abundance of low profile species; and (**c**,**d**) the same diatom guilds and evenness J.

#### *3.2. Model Results*

The above field results on nutrient concentrations and species guild and size indicate that small, fast-growing cells are also good competitors for phosphorus. This was the assumption we used in the model parameterization regarding life history traits of the initial species pool (explained in the methods above).

Model results varied with time during the simulation. At the very beginning of the simulation period, species richness started to increase, along with total biomass, and the BEF relationship was positive for these initial time steps (Figure 4a, Table 1). During succession, once all the species presented detectable biomass, no significant relationship was apparent between species richness and total biomass production (Figure 4a, Table 1). Even later in succession, when species started to go extinct and the total biomass started to reach the maximum carrying capacity of the system, there was still no significant relationship, or a negative one, between species richness and total biomass production (Figure 4a, Table 1). However, species richness and total biomass did not vary a lot between replicates at later stages of succession. On the other hand, as species started to go extinct and the system reached its maximum biomass, evenness presented a higher variability between replicates and a negative relationship with total biomass, whereas assemblages with higher dominance also presented higher biomass (Figure 4b, Table 1).

**Figure 4.** Model results between (**a**) species richness S and (**b**) evenness J (on x-axis) and biomass (y-axis) for different days (0–1000) during simulations. For each day, the 100 replicates are plotted.

**Table 1.** Slopes and their statistical significance for the equations. ln(biomass) = aS + b and ln(biomass) = aJ + b.


Note: \*\*\* < 0.001, \*\* < 0.01, \* < 0.05.

Nutrient concentrations started to decline and reached their minimum fast, with the system being phosphorus-limited early in succession, whereas nitrogen concentrations took longer to decline (Figure 5). It was during this period of nitrogen depletion that total biomass increased further and reached its maximum when both nutrients reached their minimum values (Figure 5). During succession, the species that first went extinct were the slow-growing species (Figure 6), whereas, at the end of the simulation period, the species that survived and contributed the most to the total biomass were the ones with high growth rate and high competitive ability for phosphorus (i.e., low KP) (Figure 6).

**Figure 5.** (**a**) Total biomass and nutrient concentrations; (**b**) nitrogen and (**c**) phosphorus during time succession in model simulations.

**Figure 6.** Three-dimensional diagram with the species life-history traits (μmax—maximum specific growth rate, KN-half saturation constant for nitrogen, KP-half saturation constant for phosphorus) in the model at day 0 (initial species pool), day 100 and day 1000 (end of simulation). Different colors are the three groups of species based on their competitive abilities (see methods for details).

#### **4. Discussion**

Overall, our study suggests that the BEF relationship in river benthic diatoms, although it could be regarded as positive, it seems to be a function of different factors. The main driver of the BEF seems to be geology, directly linked to nutrient availability, which, in turn, selects for specific functional traits. Furthermore, the maturity of the assemblage (i.e., time point during the succession of the biofilm) seems to be an important factor in the observed relationship, as suggested by the model simulations. This is the first attempt to generalize the BEF relationship in river benthic diatoms, using large-scale field observations and numerical modeling, and although the conclusion should be driven with caution, it offers valuable insight into these ecologically important assemblages.

The high variability among rivers regarding their environmental conditions (physicochemical parameters, nutrient concentrations, drainage area) led to variable diatom assemblages and biomass production. Therefore, it was not surprising that the BEF relationship would also vary among rivers, greatly masking the effect of diversity on biomass in the whole dataset. However, when rivers were split into two groups based on the geology of the substrate, the two patterns were very clear: positive BEF in a siliceous substrate and no BEF in a calcareous substrate. The two substrate groups differed in nutrient concentrations, which were higher in sites with a siliceous substrate. This was expected for silica, as it originates from silicate rock weathering [38]. Regarding phosphorus, this pattern has already been shown in calcareous substrates, as phosphorus is being removed from the water column due to adsorption mechanisms on carbonate material [27,39,40]. Nitrogen was also lower in calcareous substrates, although this trend was not so pronounced. Although most of the sites were phosphorus-limited, there were some sites that were nitrogen limited, belonging, though, to both geological groups. This is consistent with previous studies in Greek rivers, suggesting that the limiting nutrient was site-dependent [26].

Nutrient concentrations largely affected diatom guilds. More specifically, motile species were more abundant in increased phosphate concentrations (and thus siliceous substrates), whereas low-profile species were more abundant in low phosphate concentrations (and thus calcareous substrates). This is in agreement with previous studies, where low-profile species showed a preference for low nutrient concentrations, whereas motile species abundance started to increase with increased nutrient concentration [13]. The fact that low-profile species (i.e., species that adhere strongly to the substrate) were more abundant in calcareous substrates could indicate that species in this functional group take advantage of the precipitated phosphorus. Furthermore, most of the low-profile species found in the study (especially *Achnanthidium* spp.) have a small size, are fast growers and tend to present high populations, increasing dominance [11].

The difference in the BEF relationship between the two substrates could be explained by the combination of nutrient concentrations and traits predominance and by the maturity of the biofilm. According to field data, at higher nutrient concentrations (siliceous substrate), the addition of species could increase biomass, suggesting that species do not occupy all available niches and new arriving species make use of available space, increasing biomass [41]. On the other hand, a stable BEF relationship suggests that the species present occupy all the available niches, consuming all the available resources, and the system has reached a saturated state, even from a few species [17,18]. Model results suggest that assemblages with a positive BEF could be at an early assembly process, whereas a stable relationship could be an indication of a later in succession, more mature assemblage. This is in agreement with previous modeling studies [17] for phytoplankton, using similar models but with different parameters regarding species' life history traits. This could be an indication of a general trend in microalgae assemblages.

Model outcomes suggest that the predominant species traits are related to fast growth and strong phosphorus competition. This is related to our field results, where phosphorus limitation was more predominant, and it would select species with low phosphorus requirements [35]. The selection for fast-growing species was also highly enforced by the penalty induced in slow-growing species, a rather simplistic function that selects for

specific traits. The model applied in the present study followed many assumptions and generalizations and could not capture the complexity of a natural system. For example, nutrient inputs follow similar ratios as observed field nutrient concentrations but could be different in many cases, such as, for example, in highly polluted systems or when pointsource pollution increases the concentration of a particular nutrient [42]. However, when comparing observational and model results, it is important to remember that it is not the absolute values that are being compared but rather the trends that could give indications on mechanisms underlying observed patterns. Therefore, we believe that our model results reinforce our field findings and assumptions on the time-dependent BEF observations.

The lack of studies on the BEF relationship in benthic diatoms can be explained by a number of challenges and restrictions that it entails, some limitations of which were also apparent in the present study. Specifically, in rivers that are highly dynamic environments, biofilm assemblages can be highly affected by incidents such as heavy rains and floods and point source pollution that could make results evaluation harder. This was one of the reasons that sampling took place during summer, at low flow conditions, when there was a lower probability of heavy rain events, and we expected to collect a more mature biofilm. However, other stressors, such as pollution and desiccation, could be affecting our results [11]. Another limitation of benthic studies is the quantification of benthic concentrations and abundances and the overall sampling effort. In our study, we tried to eliminate this by scraping the biofilm of a defined surface and by using the same stone for both biomass measurement and diversity quantification. The use of chlorophyll a as a surrogate of ecosystem function is widespread in the literature, and it focuses on the biomass of primary producers of the biofilm and the general ecosystem state [43]. On the other hand, photosynthetic biofilm (i.e., phytobenthos) is a complex formation comprising many different groups of photosynthetic organisms apart from diatoms, including cyanobacteria. Therefore, different groups of species and pigments should be carefully considered in order to cover the full spectrum of the BEF relationship of the biofilm. Moreover, as water samples for the quantification of nutrients were from the water above the biofilm, and this differed from nutrient concentrations on the biofilm [44], the use of other ecosystem function metrics, such as the resource use efficiency (accounting for both biomass production and nutrient assimilation in cells, [6]), could not be directly related to our study.

#### **5. Conclusions**

This was the first large-scale field study searching for a BEF relationship in benthic diatoms in rivers. Despite the limitations recognized in a field study on benthic microorganisms, it offers important insights into species' contribution to biomass production. It highlights the importance of geology and nutrient concentrations on the form of BEF relationship and indicates species functional traits that could be responsible. The coupled modeling approach demonstrates the time-dependence of the BEF relationship during the succession of the biofilm formation and agrees with field observations on species functional traits. Further experimental work and application of different model scenarios could expand our knowledge and understanding of the ecosystem function of this ecologically important group of organisms.

**Author Contributions:** Conceptualization, E.S.; methodology, E.S.; software, G.T.; formal analysis, E.S.; investigation, E.S.; resources, E.S.; data curation, E.S.; writing—original draft preparation, E.S.; writing—review and editing, G.T. and N.T.S.; visualization, E.S.; supervision, N.T.S.; funding acquisition, E.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is co-financed by Greece and the European Union (European Social Fund-ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project "Reinforcement of Postdoctoral Researchers—2nd Cycle" (MIS-5033021), implemented by the State Scholarships Foundation (IKΥ).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are available on request.

**Acknowledgments:** We thank K. Gritzalis and I. Karaouzas for assisting in river and site selection, A. Lampou, A. Masouras and Yiorgos Maskalidis for assisting in field sampling, S. Laschou, N. Kapetanaki and G. Filippi for performing nutrient analyses, A. Abonyi and R. Ptacnik for fruitful discussions and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **Appendix A**

In Appendix A are figures and tables related to the methodology and detailed river results. All the available Figures and Tables are mentioned in the text.

**Figure A1.** Map with the hydrographic network of Greece and the sampling points for the nine rivers. Hydrographic network is from https://geodata.gov.gr/en/dataset/udrographiko-diktuo (accessed on 6 December 2022). Map was created in R (v. 4.0.3), using packages rnaturalearth (v. 0.1.0), sf (v. 1.0-6) and ggplot2 (v. 3.5.3).

**Figure A2.** Species–accumulation plots for each river. When the curve reaches a plateau, we assume that sampling effort was enough to detect most of the species present. Analysis was performed in R (v. 4.0.3), using package vegan (v. 2.5-7).

**Figure A3.** Species–area relationship for the sampled rivers. There is no clear evidence that species richness increases with increased drainage area.

**Figure A4.** Diversity (expressed as (**a**) species richness-S, (**b**) evenness-J, (**c**) functional richness-FRichness in x-axis) and biomass (ln Chl-a in y-axis) relationship for each river.

**Figure A5.** Physico-chemical (temperature, turbidity, conductivity), nutrient concentration (Si, TN, TP), diversity (species richness, functional richness, evenness) and biomass (Chl-a, Chla-c, Total biovolume) variation between the nine rivers of the study.


**Table A1.** Coordinates and sampling dates of each site. Site numbers 1–5 correspond to sites from upstream to downstream.

**Table A2.** Model parameters.


#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **Giant Panda Microhabitat Study in the Daxiangling Niba Mountain Corridor**

**Wei Jia 1,2, Shasha Yan 3, Qingqing He 3, Ping Li 4, Mingxia Fu <sup>4</sup> and Jiang Zhou 1,3,\***

<sup>1</sup> School of Life Sciences, Guizhou Normal University, Guiyang 550001, China


**Simple Summary:** The giant panda is an endemic species in China and the flagship species of global wildlife conservation. Habitat studies of the giant panda corridor can reveal their survival status, habitat environment, and the threats they face, which is crucial for giant panda population recovery and habitat conservation. The results of this study show that due to the opening of the National Road 5 (G5) Niba Mountain tunnel and the completion of the Niba Mountain giant panda corridor plan, the recovery of vegetation within the Niba Mountain giant panda corridor has led to the emergence of giant panda activity in the area, which may have spread to the central part of the reserve through the corridor. The habitat selection characteristics of the giant pandas in the corridor were clarified by investigating the microhabitats of the giant panda corridor in Niba Mountain. These findings can provide a reference for scientists to formulate practical habitat conservation and management measures for giant pandas in the study area.

**Abstract:** Habitat reduction and increased fragmentation are urgent issues for the survival and recovery of the giant panda (*Ailuropoda melanoleuca*). However, changes in the distribution and microhabitat selection of giant panda habitats in different seasons in the same region have rarely been assessed. To further understand giant panda habitat requirements, this study analyzed the giant panda habitat selection characteristics and differences using the sample data of the giant panda occurrence sites collected during 2020–2022. The results showed that the giant panda in both seasons selected medium altitudes (2000–2400 m), southeastern slopes, slopes less than 15◦, taller tree layers (8–15 m) with a larger diameter at breast height (17–25 cm) and medium density (25–55%), shorter shrub layers (<4 m) with sparse density (<30%), and taller bamboo (>2 m) with high density (>35%). The giant panda microhabitat survey in the Niba Mountain corridor clarified the characteristics of suitable habitat selection for the giant panda in the corridor. The findings of the study can provide scientific references for the development of practical habitat conservation and management measures for giant pandas in the study area.

**Keywords:** *Ailuropoda melanoleuca*; Niba Mountain corridor; suitable habitat; habitat selection; principal component analysis

#### **1. Introduction**

The giant panda is an endemic species in China and the flagship species of global wildlife conservation. Currently, the species is distributed in six mountains: Qinling [1], Minshan [2], Qionglai [3], Daxiangling [4], Xiaoxiangling [5], and Liangshan [6,7]. In recent decades, the giant panda population has decreased, with suitable habitat areas becoming increasingly shrunken and fragmented due to an imbalance between economic development and ecological conservation [8,9]. Habitat reduction and increased fragmentation are urgent challenges for the survival and recovery of the giant panda today [10].

**Citation:** Jia, W.; Yan, S.; He, Q.; Li, P.; Fu, M.; Zhou, J. Giant Panda Microhabitat Study in the Daxiangling Niba Mountain Corridor. *Biology* **2023**, *12*, 165. https://doi.org/10.3390/ biology12020165

Academic Editors: Daniel Puppe, Panayiotis Dimitrakopoulos and Baorong Lu

Received: 11 December 2022 Revised: 17 January 2023 Accepted: 18 January 2023 Published: 20 January 2023

**Copyright:** © 2023 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/).

To effectively protect the giant panda's habitats, the establishment of habitat corridors for the giant panda has been widely considered [11]. In 2007, with the support of the World Wildlife Fund, the earliest exploration of giant panda habitat corridor construction began in the Tudeling area, connecting the A and B populations of the giant panda in Minshan [12]. According to the Fourth National Giant Panda Survey (2011–2014), four giant panda corridor zones were identified in the distribution area of the giant panda. The Niba Mountain giant panda corridor zone, located in the Daxiangling Mountains, was identified as a priority ecological corridor for construction [13]. There has been no historical distribution of the giant panda in Niba Mountain, and no traces of giant panda activity were found in the area during the Fourth National Giant Panda Survey. Niba Mountain was classified as a key corridor area connecting the Daxiangling giant panda population and the Qionglai population. The construction of giant panda corridors is important for alleviating the declining genetic diversity of giant panda populations among different mountains [14]. In recent years, Gong et al. found that the habitat pattern of the giant panda is an important basis for corridor site selection. Additionally, they suggested that the study of all habitat-related microhabitat factors should be focused on [15]. Habitat studies of the giant panda corridor can more fully reveal the giant panda's survival status, habitat environment, and threats they face. This is crucial for the population recovery and habitat conservation of wild giant pandas [16].

Various constraints affect the habitat selection of the giant panda, including topography (elevation, slope, and aspect) and community structure (tree size and bamboo cover) [17]. These constraints may affect giant panda mobility, shelter availability, and the palatability of edible bamboo for the panda [18]. Previous studies showed that between 1999–2000 and 2011–2014, the giant panda in several mountains in Sichuan experienced a shift in habitat use. The giant panda increasingly utilized secondary forests, which had recovered due to protective measures. The giant panda migrated to higher elevations, despite the availability of bamboo food sources at lower elevations [10]. The implementation of natural forest conservation programs, infrastructure construction, livestock encroachment, and a range of emerging threats may have affected giant panda habitat selection [10].

Previous studies have also investigated the giant panda's habitat selection in the Daxiangling Mountains and reported changes in the habitat selection of this species between 2001 and 2020 [4,19,20]. However, these studies only described the overall habitat selection of the giant panda throughout the year. The habitat selection and utilization of the giant panda under the influence of different environmental conditions in different seasons were not further investigated. The habitat selection of the giant panda varies over time. Furthermore, they found that highly edible bamboo and good shelter sites had an important influence on the habitat selection of the giant panda [18]. The giant pandas feed mainly on *Chimonobambusa szechuanensis* and *Qiongzhuea multigemmia* at lower elevations (1500–2200 m) and *Arundinaria faberi* at higher elevations (2000–2600 m) in the Daxiangling Mountains [21,22]. However, changes in the distribution and microhabitat selection of giant panda habitats in different seasons in the same region have rarely been assessed. Determining a suitable habitat distribution pattern and selection characteristics within the protected area is essential for formulating a scientific and reasonable giant panda conservation and management plan [8]. Therefore, based on the sample data of giant panda occurrence sites collected during 2020–2022, this study aimed to apply MaxEnt and other methods to evaluate the habitat selection characteristics, including the spatial distribution pattern of suitable habitats for the giant panda in the Daxiangling Niba Mountain corridor [23]. The research conclusions can provide a scientific reference for giant panda habitat conservation and management in the area.

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

#### *2.1. Study Area*

The Daxiangling Mountains are located in the transition zone between the Sichuan Basin and the Qinghai–Tibet Plateau in the eastern part of the Hengduan Mountains, covering an area of about 6440 km2. The vegetation in the region has obvious vertical distribution zones, and the distribution order is evergreen broad-leaved forest (below 1400 m), deciduous broadleaf forest (1400–1800 m), mixed coniferous forest (1800–2600 m), coniferous forest (2600–3100 m), and alpine scrub meadow (above 3100 m) [24].

According to the results of the Fourth National Giant Panda Survey, there are 14 wild giant pandas distributed in the Daxiangling Reserve and Longcang Valley [18]. Daxiangling is the main conservation area of the Yingjing Area of Giant Panda National Park and one of the key areas of Giant Panda National Park [25,26]. The Daxiangling Niba Mountain giant panda ecological corridor (102◦29 –102◦52 E, 29◦28 –29◦43 N) is located in the southwestern part of the reserve in Yingjing County and Hanyuan County, Sichuan Province, covering an area of about 38.5 km2. The National Road 108 (G108) Ya'an section and the G5 Yaxi Expressway Niba Mountain tunnel were completed and opened to traffic in 2000 and 2012, respectively. Both cross the corridor from north to south. For the construction of these highways, many bamboo trees were cut down, which are major food sources for pandas. Furthermore, frequent human activities have caused the separation of the giant panda habitat in Daxiangling [24]. The Niba Mountain corridor has also become a key area for the exchange of giant panda populations in the Daxiangling Mountains [13].

#### *2.2. Giant Panda Habitat Selection Analysis*

From October 2020 to April 2022, the data of giant panda traces were recorded using the sample line method and infrared camera monitoring method, and a sample survey was conducted in March–April and October–November (from 2020 to 2022) in the areas where giant panda traces were recorded to determine the preferred environmental factors of the giant panda habitat. A total of 34 sample lines ≥3 km were set at intervals of ≥500 m. The sample lines covered as many vegetation types and as many potential giant panda distribution areas as possible. Combining data from 158 infrared cameras placed in the study area, the entire Daxiangling Reserve was divided into 145 square grids of 2 km2 each, with each camera spaced at least 500 m apart to ensure uniform camera coverage (Figure 1). Ten microhabitat variables were recorded in a 10 × 10 m sample square centered on the site of giant panda traces. The classification criteria for different environmental variables are shown in Table 1. A control sample was randomly set up along the sample line for every 500 m of walking or 100 m of elevation climb without traces of giant panda activity to reflect the environmental background information, and the setting and habitat variables of the control sample were recorded in the same way as the utilization sample [27]. A total of 348 samples were set up [23].

**Figure 1.** Map of infrared camera sites and giant panda occurrence sites in the study area.


**Table 1.** Classification criteria for microhabitat variables.

These habitat selection and ecological niche data were input into Excel for the relevant conversions. Following the conversion, the data were entered into SPSS13.0 for normality testing via the one-sample K-S test. Data that conformed to a normal distribution were tested through one-way analysis of variance (ANOVA), and data that did not conform to a normal distribution were tested using the Mann–Whitney *U* test.

#### *2.3. Construction of a Suitable Habitat Model for the Giant Panda in the Niba Mountain Corridor*

The estimation of suitable habitats for giant pandas in the study area was performed using the MaxEnt model. A total of 62 giant panda occurrence sites were obtained in the field, with 44 and 18 occurrence sites in the rainy and snow seasons, respectively. To reduce autocorrelation, an 1125 m radius buffer was generated in ArcGIS 10.2 with giant panda occurrence sites in the rainy and snow seasons. When the occurrence site buffers overlapped with each other, one of them was randomly retained, and the rest were eliminated, resulting in 20 and 10 occurrence sites retained in the rainy and snow seasons, respectively [27]. The giant panda has rigorous requirements for habitat, usually choosing primary forests with low human interference [28,29]. Climate and land-use types are also important factors that influence the spatial distribution of the giant panda [9,30]. Therefore, climate, topography, vegetation, and human disturbance are important factors affecting the spatial distribution of the giant panda. In the construction of the model, the above variables were selected to evaluate the habitat. Access to each variable is shown in Table 2 [23]. Because the prediction accuracy of the model was affected by the correlation between environmental factors, the "caret" function in R 4.2.1 was used to remove the highly correlated variables, and the factors with Pearson correlation coefficients greater than 0.8 were removed. Finally, nine factors were retained (Table 2).

**Table 2.** List of environmental factors for habitat suitability evaluation of the giant panda in the study area.



**Table 2.** *Cont.*

In total, 75% of the occurrence sites of the giant panda were selected for modeling, and 25% of the occurrence sites were retained for validation. The importance of each environmental factor was assessed using the Jackknife method, and the output was in the logistic format. The model prediction results were tested using the receiver operating characteristic (ROC) curve [31]. The evaluation criteria were as follows: the area enclosed by the ROC curve and the area under the curve (AUC value) was 0.5–0.6 for failure, 0.6–0.7 for poor, 0.7–0.8 for fair, 0.8–0.9 for good, and 0.9–1.0 for excellent [32]. The means of 10 calculation results were averaged to gain the habitat suitability index (HSI). The suitable habitat range in the study area was divided using Youden's index as the threshold.

#### *2.4. PCA of Microhabitat Factors for Giant Panda Habitat Selection in the Niba Mountain Corridor*

PCA can project the high-dimensional original data onto the low-dimensional mutually orthogonal principal components. This process can maximize the information content of the original data while reducing the dimensionality of the data and can effectively overcome the correlation or multicollinearity problem between variables through mutually orthogonal principal components [27]. The raw data of giant panda microhabitat variables were analyzed using PCA; the mean and covariance matrices of the sample data matrix were calculated, and the eigenvalue of the correlation matrix was found. The principal components with eigenvalues >1 were extracted, and each principal component and its contribution rate were derived from the eigenvalue to determine the factors that play a major role in giant panda habitat selection [33]. PCA was performed in IBM SPSS Statistics 26.0. In the statistical analysis, *p* < 0.05 was considered statistically significant.

#### **3. Results**

#### *3.1. Results of Giant Panda Habitat Selection Analysis*

The giant panda occurrence sites included 18 sites recorded using the sample line method and 26 sites photographed by infrared cameras in the rainy season. All occurrence sites were photographed by infrared cameras during the snow season. The total number of control samples was 155 in the rainy season and 131 in the snow season (Table 3).


**Table 3.** Comparison of microhabitat characteristics of the giant panda during the rainy and snow seasons (mean ± standard deviation).

The one-sample K–S test was performed for 10 ecological factors in the habitat selection and control plots. The results showed that six variables—namely, the altitude, average height of trees, diameter at breast height of trees, tree coverage, the average height of shrub, and bamboo height—were normally distributed, while the other four variables, including aspect, were not. The plot type was used as the control variable for the ANOVA. In normally distributed variables, there were significant differences between altitude, the average height of trees, and diameter at the breast height of trees (*p* < 0.05), but no significant differences between tree coverage, average height of trees, and the average height of bamboo (*p* > 0.05). In the control and habitat selection plots, the Mann–Whitney *U* test revealed significant differences between two variables, namely shrub cover and bamboo cover (*p* < 0.05). However, there was no significant difference between the aspect and slope (*p* > 0.05; Table 4).

**Table 4.** Comparison of variables in habitat selection plots and control plots in the study area by analysis of variance (ANOVA) and Mann–Whitney *U* tests.


\* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001.

In both seasons, giant pandas preferred medium altitudes (2000–2400 m), southeastern slopes, and areas with slopes of 15◦. The average elevation preferred by the giant panda in the snow season (2209 m) was higher than that in the rainy season (2147 m), i.e., a difference of 62 m. The giant panda was captured at an average aspect of 158.98◦ during the rainy season and 141.39◦ during the snow season. Furthermore, the giant panda was captured at an average slope of 13.41◦ during the rainy season and 8.89◦ during the snow season.

The preferred community structure of the giant panda habitat was characterized by a preference for tall (8–15 m), large (17–25 cm) diameter at breast height, and moderately dense (25–55%) trees; short (<4 m) and sparse (<30%) shrub; and tall (>2 m) and dense (>35%) bamboo trees in both seasons. The mean height of trees during the rainy season (10.98 m) and the snow season (12.28 m) were both greater than 8 m. The mean diameter at the breast height of trees during the rainy season (24.27 cm) and the snow season (17.72 cm) were both greater than 17 cm. Furthermore, the depression of trees during the rainy season (38.30%), and the snow season (40.56%), were both greater than 25%. The mean height of the shrub during the rainy season (2.47 m) and the snow season (3.58 m) were both less than 4 m. Additionally, the shrub cover during the rainy season (18.75%) and the snow season (27.22%) were both less than 30%. The mean height of the bamboo during the rainy season (2.12 m) and the snow season (2.39 m) were both greater than 2 m. Furthermore, the bamboo cover during the rainy season (59.32%) and the snow season (63.06%) were both greater than 30%.

#### *3.2. Assessment of Suitable Giant Panda Habitat in the Niba Mountain Corridor*

According to the ROC test results of the Maxent model, the AUC values were 0.964 and 0.967 for the rainy and the snow seasons, respectively, the prediction accuracy reached the level of "excellent," and the maximum Youden's index values were 0.316 and 0.527, respectively. The prediction results were transformed in ArcGIS 10.2. The results showed that the suitable habitat areas for the giant panda in the rainy and snow seasons were about 95.61 km<sup>2</sup> and 41.56 km2, respectively, accounting for 7.30% and 3.17% of the total area of the study area, respectively. In addition, most of the suitable habitat areas for the giant panda in the rainy and snow seasons were located between 1800–2600 m (mixed coniferous and broad forest). The area of suitable habitat for the giant panda in the snow season overlapped with that in the rainy season, while the overlapping area of suitable habitat in the rainy and snow seasons was about 26.85 km2, with an area overlap rate of 24.34%. Suitable habitats in the rainy season were located in the western part of the Daxiangling Reserve, in Niba Mountain, and the Longcang Valley (Figure 2), while suitable habitats in the snow season were mainly located in Niba Mountain (Figure 3). The area of the west side of G108 in the rainy and snow seasons was about 55.88 km2 and 19.86 km2, respectively, while the area of the east side of G108 was approximately 22.23 km2 and 14.12 km2, respectively. The suitable habitat during the rainy season was located in the eastern periphery of the Daxiangling Reserve in the Longcang Valley, with an area of about 17.14 km2; the suitable habitat during the snow season was also scattered in the eastern part of the Reserve, and the periphery of the Longcang Valley, with an area of about 7.07 km2. The distribution of suitable habitats in both the rainy and snow seasons showed obvious fragmented characteristics.

**Figure 2.** Suitable giant panda habitat during the rainy season.

**Figure 3.** Suitable giant panda habitat during the snow season.

#### *3.3. Main Factors Affecting Microhabitat Selection in the Giant Panda*

Among the 10 principal components of the rainy season and snow season, the eigenvalues of the first four principal components were >1, and the cumulative contributions reached 62.71% and 58.64% during the rainy and snow seasons, respectively. The first four principal components were extracted, and their corresponding eigenvectors were calculated. The variables with higher loadings in PC1 during the rainy and snow seasons were the mean height of the tree layer (10.98 m), the mean height of bamboo (2.12 m), the mean height of the tree layer (12.28 m), and the mean diameter at breast height of the tree layer (17.72%). The variables with higher loadings in PC2 were the mean height of the shrub layer (2.47 m), the shrub layer cover (18.75%), the mean height of bamboo (2.39 m), and the bamboo cover (63.06%). The variables with higher loadings in PC3 were the mean height of bamboo (59.32%) and the shrub layer cover (27.22%). Furthermore, the variables with higher loadings in PC4 for both the rainy and snow seasons were altitudes of 2147.84 m and 2209.17 m, respectively (Tables 5 and 6).


**Table 5.** Principal component analysis loadings of habitat factors at the microhabitat scale of the giant panda during the rainy season.

**Table 6.** Principal component analysis loadings of habitat factors at the microhabitat scale of the giant panda during the snowy season.


#### **4. Discussion**

*4.1. Main Activity Area of the Giant panda in the Niba Mountain Corridor*

According to the data from previous giant panda surveys, there were no traces of giant panda activity in Niba Mountain before 2015 [18]. A trail of giant panda activity was discovered in Niba Mountain in 2016. Infrared cameras began to be placed in 2017, and giant panda activity was photographed in 2018. The traces of giant panda activity were recorded during different seasons in the Niba Mountain corridor (Figure 1). Furthermore, the giant pandas were spread out in the middle of the corridor. Therefore, the giant panda may have been able to use the corridor to migrate, and the planning of the giant panda corridor in Niba Mountain is reasonable.

In 2012, the G5 Niba Mountain tunnel was completed and put into use, and the traffic flow of the G108 dropped sharply, which greatly reduced the transportation disruption in Niba Mountain. However, the vegetation in the area is still in a slow recovery stage, and no giant panda activity was found during the Fourth Survey [34]. In 2015, the planning of the Niba Mountain giant panda corridor was completed, and vegetation restoration work began in Niba Mountain [13]. Traces of giant panda activity were recorded in 2016, which indicates that it will take at least 4 years for the vegetation within an 1125 m radius of giant panda activity on the left and right sides of G108 to recover to a stage that can be utilized by the giant panda. We found that some human activities, such as bamboo shoot collection, grazing, and medicinal plant collection, occurred within an 1125 m radius of giant panda activity. However, human activities are controlled on a limited scale. Therefore, no major disturbance impacting the activities of the giant panda in the area has occurred [35].

#### *4.2. Giant Panda Microhabitat Characteristics in the Niba Mountain Corridor*

Our study showed that there were distinctive features of giant panda habitat selection in the study area. The giant panda chose to move in southeastern, gently sloping areas at medium altitudes with tall trees in both seasons, which is consistent with the results reported by Fu et al. and Bai et al. [17,18]. The giant panda always chose habitats with lower energy consumption requirements as well as higher nutritional value and net energy gain [36]. A survey of the microhabitats of the giant panda in Niba Mountain clarified that the suitable habitat characteristics of the giant panda in this area were similar to those in other areas of Sichuan Province [17,23,37].

The giant panda preferred to move at moderate altitudes in both seasons, but the average distribution elevation in the snow season was about 62 m higher than that in the rainy season. We believe that this may be related to their preference to feed on bamboo species in different seasons. The giant panda prefers to feed on *Qiongzhuea multigemmia* and *Chimonobambusa szechuanensis* at relatively low altitudes (1500–2200 m) in the rainy season and *Arundinaria faberi* at relatively high altitudes (2000–2600 m) in the snow season. This indicates that the giant panda constantly migrates according to the seasons to meet their energy needs [38], similar to the findings of Chen et al. and Liu et al. on the existence of "seasonal vertical movement" in the giant panda [39,40]. Human activities in the Daxiangling Mountains are mainly concentrated in valley areas at low altitudes (<1400 m). There are few human activities in the high-altitude areas (>1800 m) where the giant panda is active, and the giant panda can avoid anthropogenic disturbance by choosing middleand high-altitude areas as their habitat [10,11]. The results of PCA showed that altitude and bamboo growth status were significantly associated with the habitat preference of the giant panda. This is because different altitudes can provide different cumulative temperatures for the bamboo. For example, *Qiongzhuea multigemmia* only reaches the conditions for shoot development at an altitude of 1500–2200 m above sea level and when its accumulated temperature reaches a certain standard, coupled with the high ambient precipitation in March [41]. Our results indicated that altitude affects the bamboo shoot collection, the staple food of the giant panda, and thus influences the habitat selection of the giant panda in different seasons.

The Maxent model was used to predict the exact distribution and area of suitable giant panda habitats in the region during different seasons. The distribution of suitable giant panda habitats was more fragmented in snow than in rainy seasons. With the construction of the corridor, the suitable habitat area for the giant panda in the Daxiangling Reserve and Longcang Valley reached about 95.61 km<sup>2</sup> in the rainy season and 41.56 km2 in the snow season, which can promote further rejuvenation of the giant panda population. The suitable habitat area for giant pandas in the snow season overlapped with that during the rainy season. The overlapping area was about 26.85 km2, which accounts for 24.34% of the total suitable habitat area. The overlapping area, which was primarily distributed in Niba Mountain, was also the core area for the giant panda's year-round activities [40]. This is because the giant panda prefers to select different bamboo species distributed at different elevations and in different seasons. Furthermore, this verifies the research result that the average distribution elevation of the giant panda in the snow season is about 62 m higher than in the rainy season. The next step will be to introduce a finer level (four levels: most suitable, suitable, less suitable, and unsuitable) of assessment criteria for giant panda habitat assessment in the study area to explore the potential habitat and dispersal pathways of giant pandas in the study area.

The results showed significant seasonal differences in the selection of microhabitats by the giant panda, reflecting that the species has different resource requirements in different seasons [42].

#### *4.3. Newly Recorded Site for Giant Panda Activity in the Middle of the Niba Mountain Corridor*

A giant panda activity trail was discovered in the central part of the Niba Mountain corridor in 2019. Two infrared cameras captured giant panda activity in February 2022. Both sites are located within the Niba Mountain giant panda corridor, and it is possible that the giant panda in Niba Mountain is spreading westward along the corridor or the giant panda in the Longcang Valley is spreading eastward. The Niba Mountain corridor is a dispersal channel for the giant panda populations on the left and right sides, creating conditions for the exchange of giant panda populations between the two sides. However, at present, the Niba Mountain corridor is only planned to be about 5 km east of the G5 in the central part of the reserve and does not reach the Longcang Valley. Therefore, it is suggested that another ecological corridor for the giant panda can be planned from the central part of the reserve to the Longcang Valley. Alternatively, the eastern part of the Niba Mountain corridor can be extended to the Longcang Valley to strengthen the exchange of giant panda populations [18]. In 2016, a new point of giant panda activity was discovered in Niba Mountain, and in 2019, it was discovered in the middle of the Niba Mountain corridor. This indicates that it takes about 3 years for the giant panda to move from the west to the east of the reserve and from the north to the south [35].

We are aiming to better verify the actual role of the Niba Mountain giant panda corridor in the dispersal of the giant panda. In the next step, biological samples such as fresh feces will be collected from the giant panda in different areas. Then, we will explore the genetic diversity of the giant panda population by combining microsatellite and mitochondrial control regions, using feces as the main experimental material. Furthermore, we will explore the origin of newly recorded giant panda occurrence sites in the central part of the country to determine the dispersal path of the giant panda population in the Daxiangling Mountains [43].

As one of the important methods used to connect local populations of the giant panda and restore their habitat, giant panda corridors have been established in key areas under the requirements of China's giant panda habitat management plan. The establishment of these corridors has restored vegetation in the relevant areas, connected the more severely fragmented giant panda populations, and promoted exchanges between various giant panda populations [44]. The construction of the G5 Niba Mountain tunnel also provides a reference for the conservation of giant panda habitat in other mountain systems in Sichuan Province, demonstrating that it is possible to minimize the impact on the connectivity of the original giant panda habitat using elevated roadways or tunnels. This strategy not only meets the transportation development needs of human society but also maintains the normal distribution and dispersal of giant panda populations to a certain extent. We recommend monitoring the movement of the giant panda near roads and their use of corridors (e.g., road tunnels) to evaluate the impact of corridors on the giant panda. Further research should assess the effects of bamboo cover on the foraging and movement of the giant panda. In addition, it is suggested that field patrols be strengthened to decrease human activities in the giant panda habitat, including grazing, herb collection, bamboo shoot collection, and hunting. The future production and living modes of residents can be driven by ecological tourism development and the exploration of agricultural and sideline products with regional characteristics, which may change the negative attitudes of residents toward ecological conservation. This method can promote the industrial transformation to realize the balanced development of economic construction and ecological conservation [18]. Additionally, the Daxiangling Mountains should focus on the continuous conservation of major habitats for the giant panda, decrease habitat fragmentation caused by anthropogenic interference, and strengthen habitat restoration for the giant panda [23].

#### **5. Conclusions**

The present study shows that the giant panda prefers to move to medium-altitude areas; southeastern slopes; and areas with slopes less than 15◦, with an overlap of suitable habitats in both seasons. Due to the opening of the G5 Niba Mountain tunnel and the completion of the Niba Mountain giant panda corridor plan, the recovery of vegetation within the Niba Mountain giant panda corridor has led to the emergence of giant panda activity in the area, which may have spread to the central part of the reserve through the corridor. The findings can provide a reference for the future planning and construction of giant panda corridors.

**Author Contributions:** Data curation, W.J. and M.F.; formal analysis, W.J.; investigation, W.J. and P.L.; methodology, W.J.; supervision, J.Z.; validation, W.J., S.Y. and Q.H.; visualization, W.J.; writing—original draft preparation, W.J.; writing—review and editing, S.Y., Q.H. and J.Z.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA23080000).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are available on request.

**Acknowledgments:** We thank Hong Yang, Daxiangling Nature Reserve Administration, for his guidance on the fieldwork and the members of the Daxiangling Patrol and Monitoring Team for their support of the fieldwork and their valuable suggestions on data integration and processing.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
