The Effect of Environmental Factors on the Diversity of Crane Flies (Tipulidae) in Mountainous and Non-Mountainous Regions of the Qinghai-Tibet Plateau and Surrounding Areas

Simple Summary Understanding the macro pattern and underlying mechanism of species diversity is critical for research on biodiversity and biological conservation. However, there have been few reports about the different effects of the same environmental factor on biodiversity between mountainous and non-mountainous regions. This study revealed the diversity patterns of Tipulidae in the Qinghai-Tibet Plateau and its surrounding areas, investigated the influence of environmental factors on its diversity in mountainous and non-mountainous regions, and deduced the richness model of Tipulidae in mountainous regions. Our results revealed three highly endemic regions and provided a species richness model of a group of decomposer insects. The warmest quarter precipitation and topographic heterogeneity have main effects on the diversity of Tipulidae in mountainous regions. These findings provide a reference for the diversity model of decomposers to aid in the protection of biodiversity in Asian mountains. Abstract Tipulidae, one of the most diverse families of Diptera, is widely distributed in the world. The adults have weak flight ability, making it an ideal model for studying the formation of insect diversity. This study aims to explore the species diversity and endemism of Tipulidae in the Qinghai-Tibet Plateau and the surrounding areas, as well as analyze the relationships between the diversity pattern and 25 environmental factors in mountainous and non-mountainous regions. To this end, we collected 2589 datasets for the distribution of 1219 Tipulidae species, and found three areas with high diversities of Tipulidae around the QTP, including the Sikkim-Yadong area, Kamen River Basin, and Gongga Mountain. Further R, generalized additive model (GAM), and stepwise multiple regression analysis indicated that the richness and endemism of Tipulidae is mainly influenced by the warmest quarter precipitation and topographic heterogeneity in mountainous regions, but in non-mountainous regions, the richness is mostly affected by the precipitation seasonality, while there is no regularity in the relationship between endemism and environmental factors. In addition, the richness model in mountainous regions was in conformity with the results of GAM.


Introduction
Exploration of biodiversity patterns is of great significance for biogeographical research and biodiversity conservation [1][2][3]. It is also crucial to clarify the patterns of species richness and endemism and the driving forces behind these phenomena. Elucidation of the mechanisms underlying the macro patterns of species richness can help to explain species differentiation and conservation [4,5].
taset mainly included the latest checklist and distribution information of Tipulidae in the QTP and its surrounding areas. The taxonomic information was further validated using the Catalogue of the Crane Flies of the World (https://ccw.naturalis.nl/ (accessed on 25 December 2021)). The duplicated and erroneous data were excluded, and the close distribution records of the same species were removed to avoid sampling bias. For the distribution data providing longitudes and latitudes in the original sources, the georeferenced data were retained directly. For those data only providing the names of the distribution sites, the geographic coordinates were determined by using the Google map. All distribution data were kept at the county or below level. Finally, 2589 distribution datasets of 1219 Tipulidae species were obtained. The distribution data were divided into mountainous and non-mountainous groups according to the Mountain Regions of Earth provided by the Center for Macroecology, Evolution and Climate GLOBE Institute (https://macroecology.ku.dk/ (accessed on 11 February 2022)) ( Figure 1b). Distribution data of Tipulidae in the Qinghai-Tibet Plateau and its surrounding areas (Some data from Mongolia, Sri Lanka, South Korea and Russia were also included, but outside the main study areas). Green and red points represent distribution in mountainous and non-mountainous regions, respectively.

Mapping of Species Richness and Endemism
For evaluation of species richness pattern, 1° grid scale was chosen for the final analysis. The fishnet tool was employed to divide the study area into 1° × 1° (~110 × 110 km)

Figure 1. (a) Study area (malachite green part). (b) Distribution data of Tipulidae in the Qinghai-Tibet
Plateau and its surrounding areas (Some data from Mongolia, Sri Lanka, South Korea and Russia were also included, but outside the main study areas). Green and red points represent distribution in mountainous and non-mountainous regions, respectively.
The distribution data of Tipulidae were collected from three sources: (1) literature and collection sites for each species of Tipulidae based on the Catalogue of the Crane Flies of the World (https://ccw.naturalis.nl/ (accessed on 25 December 2021)); (2) specimen collection records from the Entomological Museum of China Agricultural University; and (3) distribution data from the Global Biodiversity Information Facility (https://www.gbif. org/ (accessed on 6 June 2022)) crawled by using R 4.1.2 and RStudio (package: "dismo") (some data from Mongolia, Sri Lanka, South Korea and Russia were also included) [29]. The dataset mainly included the latest checklist and distribution information of Tipulidae in the QTP and its surrounding areas. The taxonomic information was further validated using the Catalogue of the Crane Flies of the World (https://ccw.naturalis.nl/ (accessed on 25 December 2021)). The duplicated and erroneous data were excluded, and the close distribution records of the same species were removed to avoid sampling bias. For the distribution data providing longitudes and latitudes in the original sources, the georeferenced data were retained directly. For those data only providing the names of the distribution sites, the geographic coordinates were determined by using the Google map. All distribution data were kept at the county or below level. Finally, 2589 distribution datasets of 1219 Tipulidae species were obtained. The distribution data were divided into mountainous and non-mountainous groups according to the Mountain Regions of Earth provided by the Center for Macroecology, Evolution and Climate GLOBE Institute (https://macroecology.ku.dk/ (accessed on 11 February 2022)) ( Figure 1b).

Mapping of Species Richness and Endemism
For evaluation of species richness pattern, 1 • grid scale was chosen for the final analysis. The fishnet tool was employed to divide the study area into 1 • × 1 • (~110 × 110 km) grids with the furthest distribution points as the boundary based on ArcGIS 10.2 (ESRI, Inc., Redlands, CA, USA). The number of species in each grid was counted, and the endemism of each grid was determined by weighted analysis method [30]. Species richness pattern and endemic pattern were displayed in different colors.

Environmental Data
As the larvae of most Tipulidae species feed on humus and only a few species feed on living plants, they are important decomposers [23,24]. Soil organic matter content is an important driver of decomposer abundance [31]. Compared with the diversity of strictly herbivorous insects, that of Tipulidae may be less correlated with vegetation. Hence, we chose the global soil organic matter content (GSOC) instead of the vegetation factor as the contributing factor [26]. The data of GSOC were obtained from ISRIC-World (https://www.isric.org/ (accessed on 11 January 2022)). Larger-scale patterns of species diversity are generally considered to be closely related to abiotic factors [32]. In order to explore how the environmental factors drive the diversity pattern of Tipulidae in the QTP and its surrounding areas, we selected 19 climatic factors (Abbreviations: Bio1 = Annual Mean Temperature; Bio2 = Mean Diurnal Range (max temp-min temp); Bio3 = Isothermality (Bio2/Bio7) (×100); Bio4 = Temperature Seasonality (standard deviation ×100);  [33], with the spatial resolution of 30 s. For environmental heterogeneity, we selected latitude, longitude, above mean sea level (AMSL) and topographic heterogeneity as the factors. For topographic heterogeneity, we selected elevation range (elev-range) and elevation standard deviation (elev-stdev) [14]. The elevation data were obtained from the World Meteorological Database. Environmental data of all grids with a spatial resolution of 1 • × 1 • were extracted by calculating the average of each grid cell. Data extraction was performed using the software of ArcGIS 10.2 (ESRI, Inc., Redlands, CA, USA).

Data Analysis
We used Mountain Regions of Earth (https://macroecology.ku.dk/ (accessed on 11 January 2022)) as the template and crop tool in ArcGIS to distinguish the 247 grids for the mountainous group and 111 grids for the non-mountainous group. The richness index and center coordinates of each grid were exported by ArcGIS, and the endemic index was calculated by Excel's COUNTIF and SUMIF functions matched with the grid, where each distribution point was located to center coordinates by the ROUNDDOWN and VLOOKUP functions. To determine whether there are significant differences in the environmental factors driving the diversity pattern of Tipulidae in mountainous regions and non-mountainous regions, the data of both groups were analyzed using the paired normality test. If the data followed a normal distribution, the "MASS" package of "R" was used to conduct a variance homogeneity test and independent sample t-test to judge the difference [34]; otherwise, a Wilcoxon rank sum test was performed, and the significance level was set at p < 0.05.
In order to test the correlation of all environmental variables with the richness and endemism of Tipulidae in mountainous regions and non-mountainous regions, the "ggcorrplot" package of "R" was used to screen the factors highly correlated with richness and endemism [35]. The significant level of correlation was set as p < 0.05.
The screened environmental factors related to species richness and endemism were fitted by the generalized additive model (GAM) using the "mgcv" (Wood) package in "R", respectively. When there was an ambiguous relationship between the explanatory variables and response variables, the GAM was used to determine whether there was a non-linear relationship between the variables [36,37]. In this study, when the value of K = 4, and the goodness of fit reached p < 0.05, the model was taken as optimal, and an adjusted R 2 was used to quantify the variance interpretation rate of the environmental variable.
Finally, the Pearson correlation coefficients of all factors were tested to establish a prediction model for the richness pattern of Tipulidae in the QTP and its surrounding areas. Due to the involvement of many independent variables, stepwise multiple regression analysis was used to explore the relationships between multiple variables, reveal the collinearity among factors with significant influence, and verify the quality of the model according to the coefficients of the screened factors in the fitted regression equation. The screened environmental factors related to species richness and endemism were fitted by the generalized additive model (GAM) using the "mgcv" (Wood) package in "R", respectively. When there was an ambiguous relationship between the explanatory variables and response variables, the GAM was used to determine whether there was a non-linear relationship between the variables [36,37]. In this study, when the value of K = 4, and the goodness of fit reached p < 0.05, the model was taken as optimal, and an adjusted R 2 was used to quantify the variance interpretation rate of the environmental variable.

Richness and Endemism Patterns of Tipulidae
Finally, the Pearson correlation coefficients of all factors were tested to establish a prediction model for the richness pattern of Tipulidae in the QTP and its surrounding areas. Due to the involvement of many independent variables, stepwise multiple regression analysis was used to explore the relationships between multiple variables, reveal the collinearity among factors with significant influence, and verify the quality of the model according to the coefficients of the screened factors in the fitted regression equation.

Correlation of Environmental Factors with Tipulidae Richness and Endemism
Based on the results of our Wilcoxon rank sum test, Bio1, Bio4, Bio5, Bio7, Bio8, and Bio10 in mountainous regions had significantly lower values than those in nonmountainous regions (p < 0.01), while it was opposite for Bio3, Bio18, GSOC, elev-range, elev-stdev, and AMSL (p < 0.01). The remaining factors showed no significant difference (p > 0.01) (Supplementary Materials Table S1).
The GAM analysis results revealed that the species richness patterns of Tipulidae are highly correlated with its endemic patterns in the QTP and its surrounding countries and area) (p < 0.001). However, they were affected by different factors, and each factor might have different relative importance for species richness and endemism. Rainfall factors (Bio [12][13][14][15][16][17][18][19] and topographic heterogeneity factors were the most important factors, while energy factors (Bio 1-11) were the least relevant factors (Figure 3). range, and elev-stdev (p < 0.001) and related to Bio7 (p < 0.01) (Figure 3a).
In non-mountainous regions, the species richness pattern was associated with Bio15 (p < 0.01) and no factor was significantly associated with the endemism (Figure 3b). Notably, some factors with close values in the two groups showed significant correlations in mountainous regions but extremely low correlations in non-mountainous regions with species richness and endemism.
In non-mountainous regions, the species richness pattern was associated with Bio15 (p < 0.01) and no factor was significantly associated with the endemism (Figure 3b). Notably, some factors with close values in the two groups showed significant correlations in mountainous regions but extremely low correlations in non-mountainous regions with species richness and endemism.

Effects of Environmental Factors in Mountainous and Non-Mountainous Regions
According to the results of R-sq.(adj) and p-value (Supplementary Materials Table S2), the species richness of Tipulidae in the QTP and its surrounding areas was significantly influenced by Bio2, Bio7, Bio12, Bio13, Bio14, Bio15, Bio16, Bio18, elev-range, and elev-stdev.

Prediction Model of Species Richness
In this study, p-values greater than 0.7 were considered as strong autocorrelation (Supplementary Materials Table S3) [38]. A total of 13 representative factors were screened out through the Pearson correlation coefficients, including Bio2, Bio3, Bio4, Bio12, Bio14, Bio15, Bio18, Bio19, GSOC, elev-range, AMSL, longitude and latitude. We further performed stepwise multiple regression analysis on these 13 factors, and finally screened four main contributing factors, including Bio18, elev-range, longitude and Bio19. Since the species richness of Tipulidae fluctuated widely in the QTP and its surrounding areas, a logarithmic transformation was carried out on species richness, and a multiple regression equation was obtained based on correlation coefficients: y = −2.52766 + 0.000835 Bio18 + 0.000332 elev-range + 0.023789 longitude + 0.0012 Bio19 (R 2 = 0.2183, F = 16.824) (Figure 6; Supplementary Materials Tables S4 and S5). Due to the strong autocorrelation of species richness and endemism and the same GAM results, we did not predict the model of the endemism in mountainous regions. The GAM had low explanatory power for the species richness and endemism in non-mountainous regions (Supplementary Materials Table S2). Hence, we did not predict the model in non-mountainous regions.

Prediction Model of Species Richness
In this study, p-values greater than 0.7 were considered as strong autocorrelation (Supplementary Materials Table S3) [38]. A total of 13 representative factors were screened out through the Pearson correlation coefficients, including Bio2, Bio3, Bio4, Bio12, Bio14, Bio15, Bio18, Bio19, GSOC, elev-range, AMSL, longitude and latitude. We further performed stepwise multiple regression analysis on these 13 factors, and finally screened four main contributing factors, including Bio18, elev-range, longitude and Bio19. Since the species richness of Tipulidae fluctuated widely in the QTP and its surrounding areas, a logarithmic transformation was carried out on species richness, and a multiple regression equation was obtained based on correlation coefficients: y = -2.52766 + 0.000835 Bio18 + 0.000332 elev-range + 0.023789 longitude + 0.0012 Bio19 (R 2 = 0.2183, F = 16.824) (Figure 6; Supplementary Materials Tables S4 and S5;). Due to the strong autocorrelation of species richness and endemism and the same GAM results, we did not predict the model of the endemism in mountainous regions. The GAM had low explanatory power for the species richness and endemism in non-mountainous regions (Supplementary Materials Table S2). Hence, we did not predict the model in non-mountainous regions.

Diversity Pattern
Higher species richness of Tipulidae was found in the eastern regions of the QTP. Such a distribution pattern has also been observed in Hemiptera and birds [3,39]. These results confirm the important role of the QTP as a glacial refuge in Asia for most biological taxa [39][40][41]. The species richness of Tipulidae was concentrated in seven areas, which are mostly located in broad continuous mountains or block mountains in China and India.
In central and eastern China, the highly endemic areas of Tipulidae include Wuyi mountain (Fujian), Tianmu mountain (Zhejiang), and Shennongjia mountain (Hubei), which are mostly noncontiguous areas. In western China, the highly endemic areas of Tipulidae form a large continuous area, from western Sichuan, Gaoligong Mountains to

Diversity Pattern
Higher species richness of Tipulidae was found in the eastern regions of the QTP. Such a distribution pattern has also been observed in Hemiptera and birds [3,39]. These results confirm the important role of the QTP as a glacial refuge in Asia for most biological taxa [39][40][41]. The species richness of Tipulidae was concentrated in seven areas, which are mostly located in broad continuous mountains or block mountains in China and India.
In central and eastern China, the highly endemic areas of Tipulidae include Wuyi mountain (Fujian), Tianmu mountain (Zhejiang), and Shennongjia mountain (Hubei), which are mostly noncontiguous areas. In western China, the highly endemic areas of Tipulidae form a large continuous area, from western Sichuan, Gaoligong Mountains to southeastern Tibet, and is even extended to northern India and Assam. Gongga Mountain, Kamen River, and Sikkim-Yadong were three areas with the highest endemism of Tipulidae in the QTP and its surrounding areas. In Gongga Mountain and Sikkim Yadong, Tipulidae can adapt to high elevation from 1500 m to 3500 m, and is distributed from low to high elevation (1000-2500 m) in the Kamen River Basin. We imagine that the high endemism of Tipulidae in the Gongga Mountains, Kamen River, and the Sikkim-Yadong areas could be attributed to a combination of environmental factors, historical events, and research depth.
In the early 20th century, the British troops entered Tibet from the Sikkim-Yadong areas [42], and some western entomologists collected numerous specimens from Tibet [43]. In the 1930s, many missionaries (such as Armand David) and explorers found numerous biological resources in the Gongga Mountains [44]. The species of Tipulidae in the Kamen River Basin were systematically collected and published from 1961 to 1971 by Alexander [45,46]. Therefore, the high species richness and endemism in small-scale regions can be ascribed to a combination of different factors such as climate, topographic heterogeneity, and research history. The high diversity of Tipulidae suggests that dramatic speciation events might have occurred in these regions and indicates the complex research history and diverse ecological patterns around the QTP.

Correlation between Environmental Factors and Tipulidae Species Diversity
Environmental factors have inconsistent effects on the species richness of Tipulidae between mountainous and non-mountainous regions in the QTP and its surrounding areas (Figures 4 and 5, Supplementary Materials Table S2), suggesting that the richness pattern of Tipulidae is also affected by environmental factors such as elevation, monsoon, and microclimate [47,48]. Precipitation and topographic heterogeneity are the most important predictors of the diversity of Tipulidae. Although soil organic matter content is significantly different between mountainous and non-mountainous regions, it is largely irrelevant to the species richness and endemism of Tipulidae.
Most energy factors significantly differ between mountainous and non-mountainous regions (Supplementary Materials Table S1). However, it has no significant effect on the diversity of Tipulidae. The diversity pattern of Tipulidae in QTP and its surrounding areas showed that low-latitude regions have higher species richness than high-latitude regions. Hydrodynamic theory has predicted that species richness is primarily determined by the availability of water and environmental energy: the former has a more significant effect on low-latitude regions and the latter has a greater effect on high-latitude regions of the northern hemisphere [32]. Besides, the diversity of Tipulidae is less affected by warmer climates, which confirms the previous finding that Tipulidae has strong freezing resistance and can use diverse habitats for adaptation to temperature changes [49].
Tipulidae is mostly terrestrial and prefers humid habitats with a weak tolerance to drought [49]. The rich precipitation in the warmest season results in a gradual increase in the richness of Tipulidae in mountainous regions, while in non-mountainous areas, precipitation in the warmest season has no explanatory power for the richness of Tipulidae. The rainfall in the warmest season of mountainous regions was significantly higher than that of non-mountainous regions (Supplementary Materials Table S1), which may be an important driving factor for the higher richness of Tipulidae in mountainous regions.
Precipitation seasonality has a greater impact on the diversity of Tipulidae than other factors in non-mountainous regions, while it shows no obvious effect in mountainous regions. A high variation coefficient of precipitation seasonality indirectly represents instability and drastic changes of weather, which will cause a significant reduction of species diversity. However, the special environmental conditions in mountainous regions can somewhat alleviate the negative influence of dramatic weather changes, thereby better protecting the species diversity.
Climatic factors, rather than topographic factors, have the strongest explanatory power for species richness of birds, mammals, amphibians and bats [50], but topographic heterogeneity was also found to be crucial for the richness of birds in some studies [51].
Our results showed that the richness of Tipulidae sharply increases with elevation range over 4000 m, which may be associated with the special geographical conditions on the edge of the QTP; these regions tend to have very large elevation ranges and diverse climates.
Topographic heterogeneity only showed a significant effect on the species richness of Tipulidae in mountainous regions (Figure 4). In non-mountainous regions, there are also some mountains as indicated by the elevation range, but these mountains are discrete and isolated and therefore have little mountain-mass effect [52]; they can hardly play a full role as a "bridge" or "shield" [53][54][55] to weaken the impact of extreme climates to create stable microclimates and diverse habitats [50,56]. Generally, climatic stability is a key factor for high species richness in an area [57], and precipitation seasonality is more influential than topographic heterogeneity in non-mountainous regions. However, the mountain-mass effect can be enhanced by massive mountains in mountainous regions, and the warmest quarter precipitation and topographic heterogeneity are the most important factors.
In terms of endemism, topographic heterogeneity also has less significant impacts on Tipulidae in non-mountainous regions, while in mountainous regions, the greater topographic heterogeneity may bring about more complex and diverse resource heterogeneity and endemic climates [58], resulting in an increase in the endemism of Tipulidae.
In this study, the prediction model only has an interpretability of 21.83%, and therefore more factors should be considered in future studies, such as human factors and local microclimate.

Conclusions
In this study, the diversity of Tipulidae was firstly analyzed in the QTP and its surrounding areas. Three areas had high diversities of Tipulidae around the QTP, including the Sikkim-Yadong area, Kamen River Basin, and Gongga Mountain. It was found that the warmest quarter precipitation and topographic heterogeneity had significant effects on the diversity of Tipulidae in mountainous regions, and precipitation seasonality negatively affected its richness in non-mountainous regions. We speculate that Tipulidae may have sensitive responses to local microclimate, which may be a key factor driving the differentiation of Tipulidae in mountainous regions.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/insects13111054/s1, Figure S1: Relationships between endemism and Bio14, Bio15, Bio17 based on the GAM; Table S1: Differences of each factor in mountainous and non-mountainous regions; Table S2: GAM interpretability; Table S3: Correlations of each factor; Table S4: Summary of the species richness prediction model; Table S5: Correlation coefficients for the richness prediction model.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author. The data are not publicly available due to the uncompleted subject.