Geographic patterns of Lucanus (Coleoptera: Lucanidae) species diversity and environmental determinants in China

Abstract Clarifying the geographic patterns of species diversity and the determinant factors can provide essential information for species conservation and management. Stag beetles (Coleoptera: Lucanidae) of Lucanus are important saproxylic insects and can be used for biomonitoring forests. Most of Lucanus species are facing conservation concerns due to their limited distribution and fragmented habitats, particularly in China, which has the richest species diversity of this genus. The distribution patterns of species diversity of Lucanus at large spatial scales remain portly understood. We studied the distribution patterns of Lucanus and its environmental and geographic determinants in China. Distribution data for 72 species and subspecies were examined. All these species are distributed in southern China except for Lucanus maculifemoratus dybowskyi, which is mainly distributed in north China. The hotspot for Lucanus in China is southeastern Tibet. Our study indicated that the species richness of Lucanus in China was shaped by the precipitation of the wettest and driest month, net primary productivity, digital elevation model, and latitude at a large scale. These variables collectively explained 56.2% of the variation in species richness; precipitation contributed the most (44.1%). Our results provide valuable insights to improve the conservation of Lucanus and can contribute to furthering our understanding of the biogeography of stag beetles in China.

Hypotheses on the factors shaping large-scale species diversity have been proposed (Brown et al., 2004;Hubbell, 2001;Palmer, 1994). The environmental heterogeneity hypothesis based on climatic factors suggests that higher environmental heterogeneity promotes a number of habitat types and enhances species diversity (David, 1991;Lin et al., 2013;Stein et al., 2014;Tews et al., 2004).
However, studies on leaf beetles (Coleoptera: Chrysomelidae) indicated that the effects of environment and geographic position were similar (Baselga & Jiménez-Valverde, 2007).
Stag beetles (Coleoptera: Lucanidae) are saproxylic insects that commonly inhabit the lowlands or mountains dominated by broadleaf woody plants. They play important roles in the carbon cycle because their larvae live in and feed on decaying wood, and adults feed on tree sap or decaying fruit (Songvorawit et al., 2017;Tanahashi et al., 2009). Lucanidae species have been used as forest biodiversity indicators in Europe because temperature and deadwood significantly affect the presence of European stag beetles (Lachat et al., 2012;Thomaes et al., 2008). Climate has been speculated as the main driving factor for Colophon (Lucanidae) evolution (Switala et al., 2014). Morphological evidence has suggested that the evolution of mandibles of stag beetles was related to environmental heterogeneity rather than genetic differentiation (Huang & Lin, 2010;Shingleton et al., 2011). These studies on a limited number of species showed that environmental factors are critical to stag beetles.
Nevertheless, the diversity distribution patterns of stag beetles at large scales remain poorly understood. More works are necessary to explore the factors on driving distribution patterns of stag beetles. Lucanus is widely accepted as the most typical representative of Lucanidae. Studies of molecular data indicated that Lucanus diverged circa 51.3 mya (95% HPD: 48.7-53.9 mya), and its sister group is Prismognathus (Hosoya & Araya, 2005;Kim & Farrell, 2015).
Species (including subspecies) of Lucanus are particularly abundant in the Oriental region, including south and southwest China, India, Laos, Vietnam, and Myanmar (Fujita, 2010;Huang & Chen, 2010;Lin, 2017;Wan, 2007). In the present study, we examined the distribution and diversity patterns of Lucanus species in China and quantified the relative contribution of environmental conditions and spatial factors. The results will help us to understand the species diversity and distribution of these stag beetles.

| Collection of species distribution data
Distribution data of Lucanus species in China were retrieved from three sources: volumes I-III of "Stag Beetles of China" (Huang & Chen, 2010, 2017; collecting records from 1982 to 2018 of Lucanus preserved in the ecological geology specimen room at Anhui University (collected insects were not listed as "protected"); collecting records from other museums and universities (Wan, 2007) bio-nica.info/home/index.html). After removing problematic and duplicate information, we obtained 1856 distribution records for 72 species of Lucanus. Although we attempted to obtain available records for different species, the sampling effort was inevitably uneven.

| Environmental and geographic data
To analyze the association between species diversity and environmental and geographic factors, we studied 24 variables (Table S1) Data for 19 climate variables were obtained from WorldClim version 2 with a spatial resolution of 5 min (Fick & Hijmans, 2017). NDVI (Xu, 2018), NPP, and DEM data were obtained from the Resource and Environment Data Cloud Platform (http://www.resdc.cn). LAT and LON were transformed from species distribution information.
We extracted raster data in ArcGIS 10.2, including reprojection using Lambert conformal conic projection, resampling, extracting the value by distribution point, and calculating the value in each cell (500 × 500 km 2 ). Among them, 19 bioclimatic variables were char-

| Statistical analysis
Species accumulation curves (SAC) were widely used to assess the completeness of sampling effort or sampling completeness (Colwell et al., 2004;Mongombe et al., 2019). To estimate the sampling completeness of Lucanus species in China at 500 × 500 km 2 , we constructed SACs based on the number of samples (each grid represented a sample). The function "specaccum" in the package vegan (Oksanen et al., 2013) in R 3.6.1 (R-Core-Team, 2019) was used to find the SAC with the classic method "random," which adds grids in a random order. The number of permutations was set to 100.
To eliminate the effect of the square of the study area on species richness, we used an equal-area grid cell size of 500 × 500 km 2 (Ding et al., 2006;Marcelo & Douglas, 2004). We overlapped all species distributions with the 500 × 500 km 2 fishnet and then counted the number of species present within each grid as the species richness values for each grid. Drawing fishnets and calculating species richness were carried out in ArcGIS 10.2.
We filtered the 24 variables through Pearson's correlation analysis (|r| > 0.8) to eliminate the strong correlation between environmentally heterogeneous variables (Garg & Tai, 2012). Pearson's correlation analysis was performed in IBM SPSS Statistics v24. To assess the relationship between species richness and environmental variables, we applied generalized linear modeling (GLM) with Poisson regression and negative binomial regression as those models are typically used for count data. We validated the models by applying the log-likelihoods test and the null hypothesis to compare the Poisson and negative binomial GLMs. Regression models and model validation were carried out in R; Poisson GLMs were constructed by the function "glm" in the vegan package; negative binomial GLMs were constructed by the function "glm.nb" in the MASS package (Venables & Ripley, 2002); model selection process was by the function "drop1," all variables being significant at the 5% level suggested that the model selection was finished; model validation was conducted by the function "odTest" in the pscl package (Zeileis et al., 2008).

| Distribution of species diversity
The SAC illustrates trends in species richness with changing number of samples. As the number increases, the SAC increases sharply at first and then slowly, indicating sufficient sampling for data analysis

| Environmental and spatial associations
Nine variables selected through the correlation analysis, namely, annual mean temperature, temperature seasonality, precipitation of wettest month (PWM), precipitation of driest month (PDM), NDVI, NPP, LAT, LON, and DEM (Table S2). Further model selection of Poisson GLMs only retained five variables ( Table 1). All of these terms in the model were significant (p < .01). The output of the negative binomial GLM was similar to the Poisson GLM output (Table 2), except that the corresponding significance level was lower, and the AIC was slightly higher. The chi-square test statistic was equal to 1.5972, and the p-value was equal to .1031 (p > .05), which did not support the null hypothesis. Hence, there was strong support for the Poisson GLM.
To further quantify the contribution of different variables to species richness, we classified these five variables into four matrixes: precipitation (PWM, PDM), NPP, LAT, and DEM for conducting variation partitioning. The variation partitioning results showed that the total contribution of these five variables was 56.2%; precipitation (PWM, PDM) independently explained 44.1% (Figure 3), which was much higher than that of several other variables; NPP independently explained 7.5%; DEM independently explained 0.5%; LAT independently explained less than zero (−0.18%).

| D ISCUSS I ON
China is the main distribution center of Lucanus species, but many species of this genus present a narrow distribution, which may be due to a weak dispersion ability or specific habitat needs, or even the different ecological effects caused by environmentally heterogeneous habitats (Tews et al., 2004;Thomaes, 2009;Tini et al., 2017).
For example, the species distributed in Taiwan are endemic. Only L. m. dybowskyi occurs in north China. This is related to the strong diffusion ability of this species and its unique adaptability to the cold F I G U R E 1 Sample-based species accumulation curves of Lucanus species in China environment and limited broad-leaved forests (Fang & Yoda, 1991;Qiao et al., 2015).
Therefore, moderate water content enhances the occurrence and abundance of larval stag beetles (Songvorawit et al., 2017). Huang (2014 reported that climate change may threaten the survival of Lucanus species due to their specific habitat preferences. However, our model did not include temperature, which is unexpected. Lucanus inhabit mature deciduous forests, especially oak woodland, and feed on dead wood (Bardiani et al., 2017). NPP was also one of the major factors in our model, which is an important manifestation of vegetation productivity. The species richness of Lucanus can be directly affected by vegetation. NPP was positively correlated with Lucanus species richness, which means that high vegetation productivity can provide more resources, indicating that more species can be accommodated. Furthermore, studies revealed that plant cover has a significant effect on the diversity of many other groups, such as ladybugs (Coleoptera: Coccinellidae) (Sushko, 2018), stink bugs (Heteroptera: Pentatomidae) (Reisig et al., 2015), and small mammals (Keller & Schradin, 2008). However, plant diversity is under increasing threat due to human activities (Bardiani et al., 2017;Sharrock et al., 2011). Therefore, we suggest that strategies for conserving Lucanus should focus on protecting habitats. An interesting example included the use of artificial habitats comprised of rotten woodpiles and buckets with rotten deciduous leaves for larvae provides us with references (New, 2005). Besides, studies on the distribution patterns of ants (Hymenoptera: Formicidae) and ground beetles (Coleoptera: Carabidae) in China reported that changes in altitude could be a key factor affecting the distribution of species diversity (Shen et al., 2016;Yang et al., 2017). This result is consistent with our results, but the independent contribution of DEM was low.
With global climate change and large-scale deforestation, many species of Lucanus are facing problems from habitat fragmentation and loss. These species must be urgently protected. L. datunensis is the only species that is currently protected in China, and there are no protection measures or threat level assessments for the others. The distribution information for Lucanus in China is lacking. Therefore, our research can attract the attention of relevant researchers to undertake more effective monitoring and protection plans. Our analysis only found part of the relevant factors, and more factors (e.g., percentage of broad-leaved forest) and smaller scales should be explored in the future.

ACK N OWLED G M ENTS
We appreciate the laboratory colleagues Wen-Wen Chen, Qiang-Qiang Wu, and Jiao-Jiao Yuan for their kind help on data analyses.
We also thank Dr. Luca Bartolozzi, Dr. Feng Zhang, and Dr. Tian-Long Cai for providing suggestions on revision. The research was supported by the National Natural Science Foundation of China (NO. 31872276 and NO. 31572311).

CO N FLI C T O F I NTE R E S T
There have no competing interests exist.