Distribution and Assemblage Variation of Benthic Macroinvertebrates: A Uniform Elevational Biodiversity Pattern Among Different Groups?

Biodiversity patterns along the elevational gradient of vertebrates have been widely focused on in previous studies, but they are still insufficient on invertebrates in lakes to a wide elevational extent. Based on field samplings and literature, we compared biodiversity patterns among different taxonomic groups of benthic macroinvertebrates in 104 lakes of China and India along an elevational gradient of 2–5,010 m a.s.l. and revealed the key driving factors, and then, we discussed the key mechanisms underlying elevational biodiversity patterns. We found that elevational biodiversity patterns of different taxonomic groups were not uniform, e.g., an exponentially decreasing pattern of Bivalvia, a first horizontal and then decreasing pattern of Gastropoda, and a linear decreasing pattern of Oligochaeta and Insecta. Elevation and elevation-controlled variables (temperature and salinity) were the key driving factors to biodiversity patterns. Their effects were strongest on Bivalvia and less on Gastropoda, whereas they were relatively weak on Oligochaeta and Insecta. Finally, we discussed three important mechanisms that shaped elevational biodiversity patterns and assemblage variations of benthic macroinvertebrates by linking our results with the classic hypotheses about biodiversity patterns, including climate/productivity, environmental heterogeneity, and dispersal/history. These results could improve our understanding of biodiversity patterns and biodiversity conservation.


INTRODUCTION
Understanding the distribution patterns and their underlying mechanisms of biodiversity is one of the central objectives in biogeography and ecology (Gaston, 2000;Lomolino, 2001;Mittelbach et al., 2007;Laiolo et al., 2018). Large-scale spatial biodiversity patterns have recently received extensive attention, especially in the context of global climate change (Collen et al., 2014;Birrell et al., 2020). Along elevational gradients, organisms generally show a decreasing trend in species richness, and this pattern has been verified in several taxonomic groups (Bryant et al., 2008;Teittinen et al., 2017). Learning this elevational biodiversity pattern is helpful not only for determining the underlying drivers of biodiversity but also for predicting the biological response to global climate change, which is very important for biodiversity maintenance and protection (Lomolino, 2001;Rahbek, 2005;Liu et al., 2021).
In the previous studies, it was widely thought that biodiversity decreased monotonically with increasing elevation. However, recent studies showed that elevational biodiversity may show either a monotonically decreasing or a hump-shaped pattern (McCain, 2005;Rahbek, 2005;Fu et al., 2006;Beck et al., 2017;Teittinen et al., 2017). Whatever the distribution pattern has shown, climatic, area, and local environmental variables are the most cited driving factors (Lomolino, 2001;Hawkins et al., 2003;McCain, 2004McCain, , 2007McCain, , 2009Birrell et al., 2020). Besides those explanatory factors above, biodiversity distribution along elevational gradient also varies due to specific taxonomic groups (Bryant et al., 2008;Kessler et al., 2011;Nunes et al., 2020;Pandey et al., 2021). Therefore, it is necessary to study elevational biodiversity patterns on different taxonomic groups.
In this article, we studied benthic macroinvertebrates (i.e., Bivalvia, Gastropoda, Oligochaeta, and Insecta) in 104 lakes along an elevational gradient of 2∼5,010 m a.s.l. in China and India. Our objectives were to compare elevational biodiversity patterns among different benthic macroinvertebrate groups and to unravel whether the effects of elevation on the biodiversity patterns were uniform. Through explaining the relationship between elevation and elevation-related variables, we tried to reveal the mechanisms underlying these patterns. Our study could enrich the cognition about elevational biodiversity patterns of different benthic macroinvertebrate groups to a complete and wide extent and verify that those patterns were primarily determined by climate/productivity factors.

Data Source and Distribution of Studied Lakes
Study areas are located in China and India. A total of 104 lakes were studied with an elevational gradient of 2-5,010 m a.s.l. (Figure 1). In our study, benthic macroinvertebrates from 41 lakes in the Tibetan Plateau were collected from our field investigations in 2015 and 2016. To ensure the comparability of data, referring to our field investigations, we filtered the studies that used similar collecting methods and summarized the distribution data of benthic macroinvertebrates from 63 lakes, including 25 lakes from the Eastern Plain of China, 10 lakes from the Inner Mongolia-Xinjiang Plateau of China, 11 lakes from the Yunnan-Guizhou Plateau of China, three lakes from the Tibetan Plateau of China, and 14 lakes from the mainland of India (PRISMA flow diagram for literature searching and filtering is presented in Appendix 1). Locations, elevation and collecting methods, and data sources of studied lakes are presented in Appendix 2. To avoid the incomparability of taxonomic information from different sampling efforts among studies, we classified benthic macroinvertebrates groups into the uniform taxonomic level. Oligochaetes, gastropods, and bivalves were unified into genus. Insects were unified into family or genus.

Field Investigations and Source of Collected Environmental Variables
The macroinvertebrates sampling habitats included pelagic zones and littoral zones in each lake. In pelagic zones, six samples in two different stations were collected with a weighted Petersen grab (0.0625 m 2 ). In littoral zones, different habitats in two stations were randomly selected along the shoreline between depths of 15-30 cm, and two samples were collected using a D net (0.25 × 2 m, with 420 µm mesh). All the samples were sieved with a 420 µm sieve in the field. Specimens were manually picked from sediment with a white porcelain plate in the field and preserved in a 10% formalin solution. All specimens were brought back to the laboratory for microscopic examination, counting, and weighing. Benthic macroinvertebrates were identified to the uniform taxonomic level according to professional references (Liu et al., 1979;Morse et al., 1994;Epler, 2001;Wang, 2002).
In the field investigations, elevation was measured with a handheld GPS produced by GARMIN, USA. Water dissolved oxygen (DO), pH, and salinity were measured with a YSI ProPlus (Yellow Spring Inc., USA). Water transparency (transparency) was measured with a Secchi disk. Total nitrogen (TN), total phosphorus (TP), ammonium nitrogen (NH + 4 ), nitrate nitrogen (NO − 3 ), and chemical oxygen demand (COD) were measured according to the environmental quality standards for surface water of China (GB 3838-2002;Wei et al., 1989). The DO, pH, salinity, transparency, TN, TP, NH + 4 , NO − 3 , and COD of collected lakes were extracted from the literature.
Drainage area, lake area, lake length, lake width, recharge coefficient, average water depth, and precipitation were collected from in the book of Lakes of China (Wang and Dou, 1998). Annual mean temperature (Temp_mean), maximum temperature of the warmest month (Temp_max), minimum temperature of the coldest month (Temp_min), temperature annual range (Temp_range), and evaporation were downloaded from WorldClim website (http://www.worldclim.org/; Fick and Hijmans, 2017), and values were extracted at each lake scale in ArcGIS 10.1.

Statistical Analyses
First, we examined the relationship between elevation and other environmental variables with Spearman's rank correlation analysis to explore whether other variables were significantly correlated with elevation. Principal component analysis (PCA) was used to analyze the relationship among 22 environmental variables. According to the relationships above, we divided environmental variables into two groups, namely, elevationrelated variables (significant correlation with elevation), and elevation-unrelated variables (no significant correlation with elevation). Then, by exploring the relationships between species richness (uniform genus or family level) of different groups (Bivalvia, Gastropoda, Oligochaeta, and Insecta) and elevation by regression analysis and quantile regression with the packages of quantreg (Koenker, 2021) in R 3.4.3 (R Core Team, 2017), we aimed to understand the distribution variations of different groups along the elevational gradient. Canonical correspondence analysis (CCA) was used to examine the relationship between macroinvertebrates composition and the explanatory variables. Forward selection with Monte Carlo permutation tests (999 permutations) was used to select a parsimonious set of explanatory variables under the cutoff point of p < 0.05. Partial canonical correlation analysis (pCCA) was used to analyze the relative importance of elevation-related variables and elevation-unrelated variables. Then, the total variation in macroinvertebrate composition was partitioned into four independent fractions: elevation-unrelated variance, elevationrelated variance, shared variance, and unexplained variance. Some environmental variables were log-transformed before analyzing. PCA, CCA, and pCCA analyses were performed using the software CANOCO for Windows 5.0 version (Ter Braak and Šmilauer, 2012).

Environmental Variables Along the Elevational Gradient
A PCA of 22 environmental variables showed that the first ordination axis was driven by elevation and elevation-related variables (mainly climatic variables) and explained 28.6% of the variance. The second ordination axis was mainly related to the area and other environmental variables and explained 16.0% of the variance (Figure 2) TP showed negative correlations with elevation (p < 0.05). Conversely, transparency, evaporation, average water depth, pH, and salinity showed positive correlations with elevation (p < 0.05). Other environmental variables including lake length, lake width, latitude, recharge coefficient, lake area, drainage area, Temp_range, and NO − 3 showed no significant correlation with elevation (p > 0.05). Minimum, maximum, and mean value of environmental variables and their relationships with elevation are given in Appendix 3.

Taxa Composition and Distribution Patterns Along Elevational Gradient
A total of 638 taxa belonging to four classes and 106 families were recorded in our studied lakes. Insecta was the most diverse group with 353 taxa belonging to 76 families. Gastropoda was also a diverse group with 137 taxa belonging to 17 families. There were fewer taxa of Oligochaeta and Bivalvia, which were recorded as 78 taxa belonging to five families and 70 taxa belonging to eight families, respectively.
Taxa distributions along the elevational gradient were not uniform among different taxonomic groups of benthic macroinvertebrates (Figure 3). Insecta and Oligochaeta were the widely distributed groups, and some taxa only occurred in highelevational lakes. The distributions of Bivalvia and Gastropoda were relatively narrow.
The responses of benthic macroinvertebrate incidence to elevation were mostly negative, and only a few taxa showed a positive response (see Appendix 4 for details). Species richness of benthic macroinvertebrates showed decreasing tendency along elevation. However, these decreasing patterns were different among taxonomic groups: Bivalvia showed an exponentially decreasing pattern that showed a rapid decline in biodiversity at low elevation; Gastropoda richness showed a horizontal and then decreasing pattern and presented a significant decline over 2,000 m a.s.l.; and Oligochaeta and Insecta showed a linear decreasing pattern (Figure 4).

Explanatory Variables of Assemblage Variation
Climatic variables were the most important explanatory factors for assemblage variation of benthic macroinvertebrates. Variables controlled by climate such as temperature, precipitation, salinity, and pH have significant effects on taxa distribution ( Table 1). Area variables were also important factors in explaining assemblage variations, such as watershed area and water depth, which showed a significant impact on the distributions of Bivalvia and Gastropoda. Moreover, the watershed area showed a significant impact on the distribution of Insecta (Table 1). Compared with climatic variables and area variables, the impacts of other environmental variables on assemblage variation were relatively weak ( Table 1). In general, elevation and elevation-controlled variables were the key driving factors for assemblage variation of benthic macroinvertebrates ( Table 1).
FIGURE 4 | Biodiversity patterns along elevational gradient in different groups of benthic macroinvertebrates in lakes of China and India. Oligochaetes, gastropods, and bivalves were identified to genus or species. Insects were identified to family or genus. The solid lines represent the best regression curve using the least-squares method (statistical significance of r 2 values is shown). The dashed lines represent the quantile regression curve of 95th percentiles.
Frontiers in Ecology and Evolution | www.frontiersin.org TABLE 1 | Explanations (%) and order of entry of significant explanatory variables in the forward selection with Monte Carlo permutation tests to the assemblage variation of benthic macroinvertebrates (α = 0.05).

Relative Contributions of Elevation-Related and Elevation-Unrelated Variables on Assemblage Variance
Elevation-related variables explained much more variation than elevation-unrelated variables. However, the proportions of assemblage variation explained by elevation-related variables and elevation-unrelated variables were different among taxonomic groups (Figure 5). Bivalvia showed the largest proportion of variation explained by the pure effect of elevation-related variables, followed by Gastropoda, while Oligochaeta and Insecta showed a relatively small proportion of variation explained by the pure effect of elevation-related variables (Figure 5). In total, the effect of elevation on taxa distribution was strongest in Bivalvia and less strong in Gastropoda, whereas it was relatively weak in Oligochaeta and Insecta.

Elevational Biodiversity Patterns of Different Biological Groups
It is a comprehensive study that explored the biodiversity patterns of multiple taxonomic groups of benthic macroinvertebrates (including Bivalvia, Gastropoda, Oligochaeta, and Insecta) along such an extensive elevational gradient at a large geographical scale. Our results elucidated that species richness of benthic macroinvertebrates declined with increasing elevation, but the patterns varied among taxonomic groups, which showed an exponentially decreasing pattern of Bivalvia, a horizontal and then decreasing pattern of Gastropoda, and a linear decreasing pattern of Oligochaeta and Insecta. Rahbek (2005) summarized 204 data sets from 140 studies on elevational biodiversity patterns, including plants, vertebrates, and invertebrates, and assigned them to one of five patterns, namely, monotonic decreasing; horizontal, then decreasing; hump-shaped; increasing, and other. Among the patterns above, a hump-shaped elevational biodiversity pattern was more typical (c. 50%) than a monotonic decreasing pattern (c. 25%). However, the elevational range of previous studies was not wide enough (e.g., <500 m, or more than 2,000 m) to reflect the continuous biodiversity pattern. For example, a hump-shaped pattern would be concealed if elevations below 500 m were excluded (Neate-Clegg et al., 2021). And if the elevation was over 4,000 m, the monotonic increasing pattern would often hardly appear.
Our results showed a decreasing elevational biodiversity pattern in species richness, rather than the more typical humpshaped elevational biodiversity pattern in the literature, which may be caused by the studied groups. For example, nonvolant small mammals showed a clear hump-shaped elevational biodiversity pattern in species richness (McCain, 2005;Gebert et al., 2019). The hump-shaped pattern was also typical in vertebrates such as birds (McCain, 2009;Neate-Clegg et al., 2021), reptiles (McCain, 2010), and amphibians (Fu et al., 2006;Wang et al., 2019), and terrestrial invertebrate species, such as FIGURE 5 | Proportion of total explained variance related to elevation-related variance, shared variance, and elevation-unrelated variance in different groups of benthic macroinvertebrates. Oligochaetes, gastropods, and bivalves were identified to genus or species. Insects were identified to family or genus.

Key Mechanisms of Elevational Biodiversity Patterns
Three important hypotheses, namely, climate/productivity, environmental heterogeneity, and dispersal/history, have been most mentioned to explain the biodiversity patterns (Field et al., 2009;He et al., 2020). In our study, elevation-related variables (i.e., temperature and salinity) and elevation play a dominant role in biodiversity patterns, which could be corresponded with the three hypotheses above. Among them, temperature, which represents the climate/productivity hypothesis, is the most important explanatory variable in explaining assemblage variation for macroinvertebrates and different groups (Bivalvia, Gastropoda, Oligochaeta, and Insecta; Table 1). Temperature has been mostly recognized as a primary factor to shape biodiversity patterns for the groups to adjust quickly to changing climate (He et al., 2020). Temperature monotonically declines with elevation, so some tropical species find it hard to survive at high elevational areas, which leads to the decrease of species richness (Fu et al., 2004;Wang et al., 2019). In addition, low temperature could limit the species richness in two ways. One way is to reduce metabolic rates, thereby increasing generation times, and then hinder evolutionary processes (Rozen-Rechels et al., 2019). The other way is to lead to foraging restrictions in ectotherms and reduce their ability to utilize food resources (Sanders et al., 2003;Classen et al., 2015). All in all, the temperature has important impacts on almost all biological processes, such as interspecific relationships, gene mutation, speciation, extinction, and evolution. Through those ways, it affects the biodiversity patterns (Rozen-Rechels et al., 2019).
Salinity, which represents the environmental heterogeneity hypothesis, is also an important explanatory variable to explain assemblage variation in this study. Elevation generally influences the precipitation and evaporation and ultimately determines the salinity of the lakes and finally affects the biodiversity patterns (Williams et al., 1990). Several studies have shown that salinity is an ecological barrier to determining geographical distribution. With salinity increased, biodiversity declined exponentially. Lin et al. (2017) invested zooplankton communities in 45 lakes of the Tibetan Plateau and found that species richness declined significantly with salinity. In addition, our unpublished data of benthic macroinvertebrates in 40 lakes of the Tibetan Plateau showed that the taxa composition in saline lakes was significantly different from that in freshwater lakes. Oligochaeta, Gastropoda, and Bivalvia were generally distributed in freshwater lakes with salinity below 3 g/L. Chironomidae were distributed in saline lakes with salinity up to 35 g/L. Only a few taxa of Insecta (i.e., Ephydridae and Hydrophilidae) can survive in lakes with salinity over 35 g/L.
Besides the influences of temperature and salinity, elevation, which represents the dispersal/history hypothesis, also plays an important role in determining biodiversity patterns (Currie et al., 2004;Hagen et al., 2021). Elevational patterns in China were FIGURE 6 | Mechanisms of elevational effects on benthic macroinvertebrates biodiversity in lakes of China and India. "±" indicate whether there is a positive or a negative tendency between two variables. characterized by complex geological tectonic histories since the Cenozoic (Favre et al., 2015). Dispersal ability is a key condition for species distributions (e.g., Heino, 2013;Siefert et al., 2015). Island biogeography studies have shown that weak dispersers are often absent from distant islands or habitat fragments, leading to a decrease in species richness with isolation (Prugh et al., 2008;Sandel et al., 2020). Our results confirmed that this phenomenon also occurs in lake ecosystems. Most bivalves with passive movement ability are absent from high altitude areas, leading to the poor biodiversity of bivalves in Tibetan Plateau lakes. Furthermore, we have also investigated two genera of small bivalves (Pisidium and Sphaerium) in Tibetan Plateau lakes, as small bivalves may be easier to be transferred from the lowland to high plateau areas by waterbirds (Cai et al., 2018). Compared with strong aerial dispersers with flying adults (i.e., Insecta), weak and passive dispersers with aquatic adults (i.e., Gastropoda and Bivalvia) are more sensitive to elevation restrictions (Heino, 2013). This is consistent with our study. In addition, McCulloch et al. (2017) studied New Zealand's entire plecopteran fauna (100 species), observed a positive correlation between geographical range and relative wing length, and found that even small reductions in wing size could significantly influence geographical range.
In summary, the impacts of elevation on biodiversity patterns of benthic macroinvertebrates may be controlled through three important mechanisms (Figure 6): (i) with the elevation increasing, the air temperature drops sharply, leading to decreasing of the water temperature and ultimately reducing the biodiversity of benthic macroinvertebrates; (ii) precipitation and evaporation are significantly affected by the increasing elevation, which eventually increases the salinity of lakes and thereby reduces the biodiversity of benthic macroinvertebrates; (iii) due to isolation, elevation hinders the exchange and diffusion of species and ultimately affects biodiversity.
Our study concluded different biodiversity patterns for different groups of benthic macroinvertebrates along the elevational gradient. Climatic factors played a key role in the elevational biodiversity patterns in all the groups. Insecta distribute along a wider elevational gradient than Gastropoda and Bivalvia. Finally, we revealed three mechanisms that shaped the biodiversity patterns of benthic macroinvertebrates.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
BW and YC put forward the scientific hypothesis and participated in the investigations in the fields. BW, YZ, and YH identified the macrozoobenthos. BW and YH collected the literature data and analyzed the data. BW wrote the manuscript. YC, YH, and YZ edited the manuscript. All authors read and approved the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the National Natural Science Foundation of China (No. 31471962) and the National Science and Technology Basic Special (2014FY210700).

ACKNOWLEDGMENTS
We are grateful to Dekui He and Xiong Xiong from the Institute of Hydrobiology, Chinese Academy of Sciences, for the acquisitions of lake environmental data from Tibet Plateau and comments on this manuscript.