Elevational Diversity Patterns of Green Lacewings (Neuroptera: Chrysopidae) Uncovered With DNA Barcoding in a Biodiversity Hotspot of Southwest China

Elevational diversity patterns can reflect the responses of biodiversity to climate change spatially. We investigate the species diversity patterns of green lacewings (an important predatory group of insects) along the gradient of elevation from the Shaluli Mountains (Mts. Shaluli), which belong to the Hengduan Mountains in southwestern China, one of the important hotspots of global biodiversity. We combined multiple approaches, including Automatic Barcode Gap Discovery (ABGD), Assemble Species by Automatic Partitioning analysis (ASAP), General Mixed Yule Coalescent (GMYC), Poisson tree processes (bPTP), multi-rate Poisson tree processes (mPTP), to delimit the green lacewings species based on the standard barcoding region of cytochrome c oxidase subunit I (COI). The α-diversity and β-diversity patterns of green lacewings from the Mts. Shaluli along the gradient of elevation were analyzed, with further exploration on how the temperature effect elevational-diversity pattern on broad-scale (county scale) elevational gradients. The DNA barcoding reference library consisted of 40 green lacewing species from the Mts. Shaluli. The α-diversity of green lacewings decreased with the increasing elevation. The temperature was found to have a significant effect on the abundance and Shannon-Wiener diversity index but not on the species richness. Nestedness replaced turnover as the main component of Sørensen’s dissimilarity with the increasing elevation, and greater nestedness occurred at low temperature areas. The combination of a reliable DNA barcoding database could improve the accuracy and efficiency to investigate the species diversity patterns of green lacewings. Temperature, resource, and resultant interspecific competitions may have important roles in explaining the species diversity patterns of green lacewings from the Mts. Shaluli. Priority of conservation should be given to the species at low elevation, middle elevation, and relatively high temperature regions under the background of global climate warming.


INTRODUCTION
Elucidating the factors that drive global biodiversity is the fundamental goal in ecology and conservation (McCain, 2007) considering the complexity of biotic (i.e., interspecific interactions, food availability, and thermal tolerance) and abiotic (i.e., temperature, precipitation, atmospheric pressure, and wind velocity) factors in nature's ecosystem (Sundqvist et al., 2013). Temperature is commonly considered to be a main causal factor of species diversity dynamics (McCain, 2007;Deutsch et al., 2008;Diamond et al., 2016;Birkett et al., 2018;Barbarossa et al., 2021). The elevational gradients, which may reflect their responses to climate change spatially, also play an important role in providing powerful information on how the ranges of biodiversity are restricted by the environmental conditions (Grinnell, 1914;Sundqvist et al., 2013;Supriya et al., 2019). Although a variety of elevational diversity patterns are found among the various taxa and from different regions, exploration of the unknown patterns of diversity along the altitude gradients and their formation mechanisms are still ongoing projects. So far, the elevation patterns of overall community variation (α-diversity) are known as follows: (1) α-diversity declines with increasing elevation, (2) α-diversity shows a hump-shaped or bimodal-shaped distributions with increasing elevation, (3) α-diversity increases with increasing elevation, and (4) α-diversity has no linear relation with elevation (Rahbek, 2005;Peters et al., 2016;Noriega and Realpe, 2018;Supriya et al., 2019;Uhey et al., 2021). The community composition variation (β-diversity) has been interpreted as the following two different features. First, the nestedness (species are lost or gained) is prevalent at high elevation and low elevation, or only present at one of them, while the turnover (a species replaced by other species) is common at mid-elevation (da Silva et al., 2018;Uhey et al., 2021). Second, only turnover is present in all the elevational gradients (Paknia and Sh, 2015;Young et al., 2019). Combining both the α-diversity and β-diversity gives a more comprehensive understanding of the biotic heterogeneity of a certain area (Wang et al., 2000;Baselga, 2010).
The insects represent a major component of global biodiversity and are widely used as an indicator of biodiversity in ecological studies. The previous studies on the insect biodiversity along the elevational gradients are confined to a small number of groups with ecological significance, such as geometrid moths (Axmacher et al., 2004;Beck and Chey, 2008), butterflies (Despland et al., 2012), pollinating hymenopterans (Jackson et al., 2018;Perillo et al., 2021), ants (Reymond et al., 2013;Bishop et al., 2014), dung beetles (da Silva et al., 2018;Stanbrook et al., 2021), and ground beetles (Zou et al., 2014). Moreover, the studied taxa have distinctly differed distribution patterns along elevation (Rahbek, 2005). The green lacewings (Neuroptera: Chrysopidae), a group of important predatory biocontrol agents with a positive role in the management of aphids, coccids, thrips, whiteflies, and mites in the ecosystem (Alford, 2019;Lai and Liu, 2020), had rarely been investigated. The previous studies on the diversity of green lacewings mainly focused on the species richness and/or Shannon-Wiener diversity index in agro-ecosystem Thierry et al., 2005;Martins et al., 2019) or forest ecosystem (Yi et al., 2018), but rarely investigated on the pattern of β-diversity. For example, Bozdogan (2020) preliminarily explored the richness patterns of a green lacewing community along the elevational gradients from the East Mediterranean but did not document the β-diversity of green lacewing community along the elevational gradients and correlated temperature variation up to now. On the other hand, understanding the taxonomic diversity, distribution pattern, and driving factors of the green lacewings are all critical for the exploration of the dominant species as the natural enemies in the different agricultural ecosystems  and developing conservation strategies (Zou et al., 2014). However, the diversity of green lacewings had been rarely taxonomically known, and their distribution patterns and responses to the environmental factors were rarely investigated, particularly, in the remote mountainous region of Southwestern China (e.g., the Hengduan Mountains), one of the most important global diversity hotspots, where the complex topography and climate were thought to be in favor of high level of biodiversity (Wu et al., 2013;Wen et al., 2016).
In the recent two decades, owing to the molecular species delimitation accelerated by DNA barcoding (Hebert et al., 2003;Ruiter et al., 2013;Hendrich et al., 2015), the largescale biodiversity investigation has broken through the limit for analyzing the diversity at the family level (Negi and Gadgil, 2002;Zou et al., 2020), but is turning to focus further on the species richness (Smith et al., 2005;Woodcock et al., 2013;Young et al., 2019). Nevertheless, accurate molecular species identification relies on the establishment of the DNA barcoding library of the specific groups and the development of an appropriate analytical tool for species delimitation (Ruiter et al., 2013;Puillandre et al., 2021). So far, various approaches [e.g., General Mixed Yule Coalescent (GMYC), Automatic Barcode Gap Discovery (ABGD), Poisson tree processes (bPTP), multirate Poisson tree processes (mPTP), and Assemble Species by Automatic Partitioning analysis (ASAP)] are mainly applied in the species delimitation based on DNA barcoding. GMYC and bPTP estimate the absolute time and mutational time, respectively, at different nodes of the phylogenetic tree (Pons et al., 2006;Zhang et al., 2013;Puillandre et al., 2021). ABGD automatically identifies the limit between the smaller intraspecific distances and the larger interspecific distances based on the distribution of pairwise genetic distance within priori maximal genetic intraspecific divergence P (Puillandre et al., 2012(Puillandre et al., , 2021. mPTP alleviates the theoretical and technical shortcomings of bPTP and takes account of divergent intraspecific variation. ASAP proposes the species partitions are ranked by a scoring system that does not need the biological prior insight of intraspecific diversity (Puillandre et al., 2021). However, the resulted divisions of potential species inferred by these methods are sometimes inconsistent, which thus calls an integrative approach combining these methods and morphological evidence as well (Ducasse et al., 2020;Puillandre et al., 2021).
Based on a 2-year survey on the lacewing diversity from the Shaluli Mountains (Mts. Shaluli), we aimed to (i) explore the species diversity of green lacewings from the Mts. Shaluli by using multiple approaches of species delimitation (GMYC, ABGD, bPTP, mPTP, and ASAP) based on the DNA barcodes, and to (ii) reveal the community diversity patterns of green lacewings along elevation and the effect of temperature on their broadscale (county scale) community diversity. We predict that (i) α-diversity may reach the peak from the mid-elevation areas, as the green lacewings are apt to suffer human disturbance from the low-elevation areas (Bozdogan, 2020;Dupont and Strohm, 2020;Lai and Liu, 2020), and the fitness to the habitats from the high-elevation areas is low; (ii) α-diversity may increase with the increasing temperature; (iii) greater nestedness component of β-diversity may occur at sites with a high elevation and low temperature as many thermophilic organisms are limited at higher elevations by their inability to withstand cold temperature (Fu et al., 2007).

Study Area and Sampling
The Mts. Shaluli is one of the major parts of the Hengduan Mountains, located at the west of Sichuan, southwestern China, and is featured in a complicated topography, with striking vertical zonation, ranging from a dry-hot valley canopy to alpine permafrost. The climates vary dramatically, with the difference of mean annual temperature reaching up to 20.2 • C (Wu et al., 2013;Wen et al., 2016). The high floral diversity and wide range of habitat types make the Hengduan Mountains form a hotspot of biodiversity (Yao et al., 2015). Additionally, the Mts. Shaluli is a key area producing apple and pear in Sichuan Province (Xie et al., 2011). With conspicuous elevation gradients and rich species diversity of the green lacewings, the Mts. Shaluli is an ideal candidate for the regional ecosystem for exploring the responses of green lacewings to environmental change, and these findings will contribute to the bio-control application using the local green lacewings.
All the green lacewing specimens herein studied were collected from July to August of 2019 and 2020 from 48 sampling sites of 113 sampling sites (27.3094-30.2623 N, 99.0108-101.7959 E) in eight counties within the Mts. Shaluli (Figure 1 and Supplementary Table 1), containing six important reserves (Zhubalong Nature Reserve, Batang; Xiayong Nature Reserve, Derong; Fozhuxia Nature Reserve, Xiangcheng; Gexigou Nature Reserve, Yajiang; Haizishan Nature Reserve, Litang; and Qialangduoji Nature Reserve, Muli) and its adjacent regions. The number of sampling sites of every county varied from 11 to 22 depending on the size of the area of these nature reserves. The specimens were collected by sweep-nets with a diameter of 50 cm from 9 am to 6 pm within a day, and also by light trap with 450 W high-pressure mercury lamp from 8 pm to midnight. No green lacewings were collected from the remaining sites. The habitat types of the sampling sites comprise (1) croplands, (2) broadleaved deciduous forests, (3) residential areas, (4) mixed coniferous broad leaved forests, (5) mixed coniferous broad leaved forests and residential areas, (6) coniferous forests, (7) alpine shrubs, (8) grasslands, (9) wetland, and (10) rivers, ranging from 1,338 to 4,525 m in altitude. All the specimens are preserved in 95% ethanol or pinned, and preserved at −20 • C or room temperature, respectively. All the specimens are deposited in

Preliminary Species Identification Based on Morphology
The green lacewing specimens were identified to species or morphospecies based on the morphological characters by using the published keys to genera and species (Brooks and Barnard, 1990;Yang et al., 2005). The habitus photographs of the identified specimens are deposited and can be accessed in the BOLD database (dx.doi.org/10.5883/DS-CXB1).

DNA Barcoding
We selected 3-5 individuals per species or morphospecies to obtain the DNA barcodes. For the species that are difficult to be distinguished or have distinct potential intraspecific variations, over five individuals were selected for the analysis. The genomic DNA was extracted from the hind leg of each specimen using the TIANamp Genomic DNA Kit (TIANGEN Inc., Beijing, China). Because the universal DNA barcoding primers usually applied for the insects (i.e., LCO1490 and HCO2198; as shown in the reference Folmer et al., 1994) did not work well for the green lacewings in our study, the barcoding region of the mitochondrial COI gene (Hebert et al., 2003) was amplified using primers COI1-F (5 -ATTCAACCAATCATAAAGATATTGG-3 ) and COI1-R (5 -TAAACTTCTGGATGT-CCAAAAAATCA-3 ), which were accessed from the data named as "mitochondrion Chrysopa oculata (golden-eye lacewing)" in NCBI 1 . In case of unsuccessful amplification, an alternative forward primer COI4-F (5 -TGTAAAACGACG-GCCAGTAAACTAATARCCT TCAAAG-3 ) 2 (reference to Park et al., 2010) was applied to amplify the selected region. Each 25 µl volume contained 12.5 µl 2 × EasyTaq PCR SuperMix (TransGen Biotech, Beijing, China), 9.5 µl ddH 2 O, 1 µl of each primer, and 1 µl DNA template. The thermal cycling conditions followed Yi et al. (2018), but the annealing temperature was 52 • C for COI4-F/COI1-R and 50 • C for COI1-F/R. The PCR products were sequenced on ABI3730XL Genetic Analyzer by Beijing Ribio BioTech Co., Ltd., Beijing, China.
The sequence reads of COI in forward and reverse directions were assembled by ContigExpress (NY, United States), and we trimmed the tRNA-W region that owned stop codon at <200 bp of the 5 end sequence based on the translation of amino acid. The sequence of standard DNA barcoding region of Chrysoperla nipponensis (trimmed from the complete mitochondrial genome DNA of this species, GenBank accession number: AP011623.1) was selected to join in the multiple sequence alignment in MEGA v.7.0 (Kumar et al., 2016) to confirm the tRNA-W region completely trimmed.
To assess the efficiency of our DNA barcoding database, the DNA barcoding gap was explored by calculating the nearest neighbor (NN) distance and maximum intraspecific (MI) distance based on a Kimura 2-parameter (K2P) model in MEGA v.7.0. The sequences with at least 500 bp were retained for the species delimitation.
The species delimitation based on the DNA barcodes was conducted by using the five methods as follows. First, based on the pairwise genetic distances, we performed the ABGD analysis (Puillandre et al., 2012) using the webserver 3 . The settings were as follows: relative gap width X = 1.0, K2P distance, and the prior maximum divergence of intraspecific diversity (P) values ranged from 0.001 to 0.100, other parameter values employed defaults. Second, we performed the Assemble species by ASAP (available at) 4 , which is also a K2P distance-based approach. Partition with the smallest score is considered as the final result of the species delimitation (Puillandre et al., 2021). The remaining three analyses, i.e., GMYC (Pons et al., 2006), PTP , and mPTP (Kapli et al., 2017), are phylogeny-based approaches. The maximum-likelihood (ML) tree was generated by using an IQ-TREE webserver (Nguyen et al., 2015;Trifinopoulos et al., 2016) based on the COI sequences under the GTR model 5 . The trees were saved as Newick format in the software FigTree v.1.4.3 (Rambaut and Drummond, 2016), then imported into the PTP 6 and mPTP 7 web server to delimit the species. For the GMYC analysis, the ultrametric trees were reconstructed with BEAST v.2.6.3 (Bouckaert et al., 2014). The relaxed lognormal clock model, Coalascant Constant Population tree priors were used to estimate the relative divergence times (Monaghan et al., 2009;Michonneau, 2015). The Markov Chain Monte Carlo (MCMC) generations were set to 200,000,000. TRACER v.1.7.0 (Rambaut et al., 2018) was utilized to check whether all the effective sample size (ESS) values have exceeded 200 for assessing the convergence of runs. The consensus trees were generated after discarding 25% of trees as burn-in (Puillandre et al., 2021).

Diversity Analysis
For the analysis of the diversity patterns of green lacewings along elevation, to avoid environmental heterogeneity among the different counties, 14 sampling sites within the same county (i.e., Derong County) with obvious elevational gradient (2093-4118 m) were selected to analyze (Supplementary Table 2). The elevational gradients of Derong County (range from 2,093 to 4,188 m) were partitioned into lowelevation (below 2,500 m), mid-elevation (2,501-3,500 m), and high-elevation (more than 3,501 m). For the analysis of the temperature effect on broad-scale (county scale) community diversity of the green lacewings, as it was considered as a key driver of an elevational diversity pattern (Bishop et al., 2014), we analyzed the data from the 33 sampling sites among the six sampling counties, because their historical daily maximum temperature, daily minimum temperature, and daily mean temperature were available in the China Meteorological Data Service Centre 8 (Table 1), and the mean elevation of these The correlation between α-diversity (abundance, species richness, and Shannon-Wiener diversity index) (Kaltsas et al., 2018) of green lacewings community and the variate we analyzed (elevation and temperature) were assessed in SPSS v.19.0 (IBM, NY, United States) with Pearson's correlation coefficient. If there is a correlation, the linear regression model (y = bx + a) or quadratic regression model [(y = cx 2 + a) or (y = cx 2 + bx + a)] will be used to fit (Werenkraut and Ruggiero, 2014). The goodness of fit was evaluated based on the adjusted R 2 (Yang et al., 2020). The species component difference between pairs of sites was assessed by β-diversity. The Sørensen dissimilarity index (βsor) based on the species richness (Baselga, 2010) was used to measure the overall β-diversity (Fu et al., 2019). Following Baselga (2010), βsor was partitioned into Simpson dissimilarity index (βsim) and nestedness-resultant dissimilarity (βnes). To determine the most important component, the comparisons between βnes and βsor were performed for pairs of sites (Dobrovolski et al., 2012). The formulas are as follows: where a is the number of species common to both sites, b is the number of species that only occur in the first site, and c is the number of species that only occur in the second site (Baselga, 2010).
To evaluate whether the α-diversity and β-diversity between the different elevational gradients are significantly different, the ANOVA of the Shannon-Wiener diversity index and βsor among the different elevational gradients would be conducted in SPSS v.19.0 with the generalized linear models, a paired comparison was performed by Tukey's significant difference post-hoc test (Yang et al., 2020).

Species Delimitation
A total of 1,168 specimens of Chrysopidae, belonging to 12 genera and 40 morphospecies (Supplementary Table 1), were collected from the Mts. Shaluli. In total, 246 specimens of 40 morphospecies were selected and their DNA barcodes (547-705 bp) were successfully obtained ( Table 2; GenBank accession number: MZ557095-MZ557340). There are 181 (73.58%) barcodes from 22 morphospecies identified as the described species. The remaining 18 morphospecies were identified only to the genus.
After alignment and trimming, there are 245 barcodes with 557 bp used for the calculation of the genetic distances. A distinct barcode gap was presented in 93.10% of all the analyzed species (Figure 3). The minimal barcode gaps were presented in Ch. annae and Ch. nipponensis. The MI distance was close to the interspecific distance between NN species (MI vs. NN: 0.016 vs. 0.013 for Ch. annae; 0.016 vs. 0.009 for Ch. nipponensis) (Supplementary Tables 3, 4).

Species Diversity
Among the 40 species herein delimitated, the four species are newly recorded in China, and 11 species are newly recorded in Sichuan Province (Supplementary Table 5). Apertochrysa owns the highest richness (22.50% of species), followed by Chrysopidia (17.50% of species). Additionally, Apertochrysa owns the highest abundance (33.73% of individuals), followed by Chrysopa (17.51% of individuals), and Cunctochrysa (16.35% of individuals). Three dominant species from the Mts. Shaluli was recognized and listed as follows (Figure 4). Apertochrysa barkamana (Yang et al.), accounting for 20.98% of the total number of individuals, is a dominant species in Yajiang County with elevation

Elevational Diversity Pattern
Given the obvious elevational gradient (2,093-4,118 m) and sufficient sampling, Derong County was selected for the analysis of the elevational diversity pattern of the green lacewings. Based on the above species delimitation, there are 17 species (a total of 249 specimens) from Derong (Supplementary Table 2). The abundance, richness, and Shannon-Wiener diversity index of the green lacewing communities from Derong County were significantly correlated with elevation (R 2 = 0.66, df = 11, P < 0.01; R 2 = 0.53, df = 12, P < 0.01; and R 2 = 0.34, df = 12, P < 0.05, respectively), and have the same elevational-diversity pattern, being decreased with the increasing elevation, as shown by the better linear regression model (Figure 5A), but unlike richness and Shannon-Wiener diversity index, the abundance of the green lacewing communities was not monotonically decreasing, and its rate of decrease gradually flattened out. The ANOVA result (Table 3) showed that the Shannon-Wiener diversity index of the green lacewing communities in high-elevation was significantly different from that in the low-elevation and mid-elevation, and there was no significant difference between the low-elevation and mid-elevation. There was a difference in the Sørensen dissimilarity index (βsor) based on the species richness of the paired sites in the different elevational gradients, but their differences among the elevational gradients were not significant ( Table 3). Nestedness-resultant dissimilarity (βnes) became a major component of β-diversity with the increasing elevation (Figures 5B,C).

Species Diversity and Temperature Variables
The abundance of the green lacewing communities from six counties (each herein defined as a temperature block; Table 1) located at the Mts. Shaluli significantly increased with the increasing temperature ( Figure 6). However, the richness was not significantly related to temperature in the present results. Shannon-Wiener diversity index of the green lacewing communities was positively correlated with daily minimum temperature and daily mean temperature, but not daily maximum temperature (Figure 6). Compared with the other temperature blocks (Table 1 and  Supplementary Table 6), the ratio of βnes to βsor from the lowtemperature blocks (i.e., Daocheng County and Litang County; daily maximum temperature below 21.40 • C, daily minimum temperature below 7.10 • C, and daily mean temperature below 13.50 • C) was 1 or 0, indicating that low temperature could lead to the nesting of communities or the absence of common species shared among the communities. Besides, the results show that when the difference of daily minimum temperature or daily mean temperature between the paired temperature blocks exceeds 6 • C, the difference of species components or nestedness between the communities reaches the greatest level.

Barcoding of Green Lacewings
The green lacewings represent one of the most diverse lineages of Neuroptera. Although they are commonly found in various ecosystems, such as forests, orchards, and croplands, the identifications of many green lacewing species, even of some frequently used natural enemy species, are difficult to practice solely based on the morphological characters because their interspecific differences are indistinct (Pappas et al., 2011;Lai and Liu, 2020). So far, there have been few studies that attempted to identify the green lacewings using the DNA barcodes (Morinière et al., 2014;Yi et al., 2018). Our analysis of the green lacewing barcodes from the Mts. Shaluli shows a distinct barcode gap presented in 93.10% of all the analyzed species (Figure 3), allowing unambiguous identification of 92.5% species. It suggests that DNA barcoding is promising for the accurate identification of the green lacewing species. Nevertheless, a well-established library of the DNA barcodes is an important precondition, which requires species delimitation facilitated by integrative data from morphology, molecule, and distribution.
Our results (Figure 2) suggest that mPTP is too conservative to be implemented in species delimitation of the green lacewings based on the DNA barcodes. In addition, ABGD and ASAP similarly performed as mPTP, being relatively conservative to delimitate closely related species. Here, GMYC and bPTP performed better because they could correctly distinguish the most closely related species or even cryptic species. A similar case is found for the species delimitation of Psylloidea (Hemiptera: Sternorrhyncha) (Martoni et al., 2018). Accordingly, for the common and species-rich genera Chrysoperla, Apertochrysa, and Mallada, GMYC and bPTP are recommended for molecular species delimitation even though they are time-consuming (7 s for ABGD and ASAP, 7 h for bPTP, and 33 h for GMYC in this study). When delimitating closely related species, however, approaches combining morphological and ecological evidence are needed (Puillandre et al., 2012;Huang et al., 2020).

Patterns of Elevational Diversity
The previous studies on some organisms from the Hengduan Mountains (Gong et al., 2005;Rahbek, 2005;Fu et al., 2006Fu et al., , 2007Li et al., 2009;Wu et al., 2013;da Silva et al., 2018;Uhey et al., 2021) found that the hump-shaped elevational α-diversity patterns are most prevalent, while a monotonic FIGURE 5 | The elevational diversity patterns of the green lacewing community from Derong County. (A) Relationships between elevation and α-diversity. The linear and quadratic regression were used to fit the elevation pattern. The goodness of fit was assessed based on the adjusted R 2 . **, *, and ns indicate significance at P ≤ 0.01, P ≤ 0.05, and non-significant, respectively. The solid lines indicate linear regression and the dashed lines indicate quadratic regression. Dissimilarity analysis (Sørensen dissimilarity index, βsor) of the green lacewings from Derong Country partitioned into Simpson dissimilarity index (βsim) and nestedness-resultant dissimilarity (βnes) of (B) elevational site pairs (i.e., each site and its immediately adjacent uphill site), and (C) elevational gradients.
decreasing pattern with increasing elevation is also frequently reported (Rahbek, 2005;Sam et al., 2015;Noriega and Realpe, 2018;Finnie et al., 2021). The richness and Shannon-Wiener diversity index of the green lacewings from Derong County are consistent with the latter pattern mentioned above, the relationship of the abundance of green lacewings and elevation fit with the right half of hump-shaped (Figure 5A), which contradict our hypothesis of highest α-diversity from the midelevation areas. This may be explained in the following two aspects. First, incomplete sampling from the areas with different elevations could result in biased patterns of elevational diversity. Second, the definition of low-or mid-elevation may vary among the different regions. Here, low-elevation (2,000-2,500 m) is frequently defined as mid-elevation (Peters et al., 2016;Supriya et al., 2019). However, Bozdogan (2020) found the monotonic decreasing pattern of green lacewing species-richness along the increasing elevation from the East Mediterranean area of Turkey, which has a relatively lower elevation ranging from 400 to 1,400 m. Moreover, the elevational diversity is found to shift from a unimodal pattern toward a monotonic decline with increasing taxonomic coverage of the plant and animal communities (Peters et al., 2016). Thus, the patterns of elevational diversity appear to be variable rather than uniform among taxa (Rahbek, 1995). Different taxa, seasonality, climate, resource availability, and interaction between species (Beck et al., 2010;Sundqvist et al., 2013;Supriya et al., 2019) may contribute to the differentiation of diversity patterns. Concerning the presently revealed decline of α-diversity of the green lacewings along increasing elevation, the decline of temperature and the decrease of resources with higher elevation may have major effects on driving this pattern, because the lower temperature could hamper the growth and induce diapause of the green lacewings (Kostál, 2006;Wang et al., 2019), and the cannibalism among green lacewing larvae could be provoked when the preys are scarce (Mochizuki et al., 2006;Ye and Li, 2020). Being components of β-diversity, turnover and nestedness may better reflect the extent of changes in community composition among the sites (Baselga, 2010). We found that species turnover decreased and nestedness increased with the increasing elevation ( Figure 5). Herein, turnover is the main component of β-diversity from low and middle elevation. It is supported by the distributions of most green lacewing species, which are not overlapped on a small scale at low and middle elevations. However, the difference of the Sørensen dissimilarity index (βsor) based on the green lacewing species richness among the elevational gradients is not significant ( Table 3). In other words, the turnover from low-and mid-elevation has a similar FIGURE 6 | Relationships between temperature and α-diversity of the green lacewing community from six counties (Batang County, Derong County, Litang County, Muli County, Yanyuan County, and Daocheng County). The historical temperature data of the above counties are from China Meteorological Data Service Centre. The linear and quadratic regression were used to fit the temperature pattern. The goodness of fit was assessed based on the adjusted R 2 . **, *, and ns indicate significances at P ≤ 0.01, P ≤ 0.05, and non-significant, respectively. The solid lines indicate linear regression model and the dashed line indicates quadratic regression.
effect on the species composition as nestedness from a high elevation. Turnover is commonly due to species competitive exclusion and/or niche specialization (Ulrich et al., 2009;Baselga, 2010). Because there are a relatively high richness and abundance of the green lacewing communities at low-and middle-elevation, species competitive exclusion may be common in limited space at low-and mid-elevation for the predatory green lacewings, eventually leading to higher turnover in the two elevations. Other studies also found that the nestedness component is greater in higher latitude areas (Dobrovolski et al., 2012;Paknia and Sh, 2015), which may be relevant to physiologic limitation as it is too cold to survive for most green lacewings at high elevation.

Temperature Shaping Species Diversity
Temperature is considered to be an important factor to explain population dynamics (Geurts et al., 2012) and an important predictor of species richness along with the elevational gradients (Fu et al., 2007;Peters et al., 2016). In this study, the species diversity of the green lacewings is strongly affected by the temperature on a broad-scale elevational gradient. The relationship between the temperature and α-diversity of the green lacewings community was incompletely consistent with our hypothesis of α-diversity increases with the increasing temperature. We found that the abundance of green lacewings significantly increased with the increasing temperature (Figure 6). This finding further confirms that the terrestrial ectotherms may experience an initial increase in the population growth rates at mid-to high-latitudes with climate warming (Deutsch et al., 2008). Considering the richness of the green lacewing communities and temperature, there was no significant correlation. It indicates that other factors may independently or together with temperature change the richness of the green lacewings communities in different temperature blocks, such as water availability and geographic distance (Qian et al., 2005;Peters et al., 2016;Liu et al., 2020). The Shannon-Wiener diversity index of the green lacewing communities was positively correlated with daily minimum temperature and daily mean temperature, but not daily maximum temperature (Figure 6). The previous studies (Araújo et al., 2013;Liu et al., 2020) suggested that physiological tolerance to cold evolves more quickly than tolerance to heat in plants and animals. As such, the organisms are expected to have limited plasticity in the increase of their upper thermal limits (Bai et al., 2019), explaining more variability of the green lacewing species diversity among the cooler areas than that among the warmer areas.
The higher nestedness component of β-diversity occurred at the low-temperature blocks (Supplementary Table 6), which is consistent with our hypothesis. We found only three species from the low-temperature blocks, even though we made extensive collecting from 22 sampling sites within the two low-temperature blocks (Supplementary Table 1). The low-temperature blocks are located around Haizishan Mountain and its adjacent areas and their elevational range are above 4,000 m, where the glaciation was particularly active during the Pleistocene and massive extinction of thermophilic organisms, such as green lacewings might have taken place (Cao and Zhao, 2007;Dobrovolski et al., 2012), as many thermophilic organisms are limited at higher elevations by their inability to withstand cold temperature (Fu et al., 2007).

Conservation Implications
As the increase in temperature above optima has negative effects on the organisms (Bozinovic et al., 2011), the species with limited thermal tolerance usually undergo uphill shifts or extinct with climate warming (Chen et al., 2009;Birkett et al., 2018;Stanbrook et al., 2021). Particularly, the tropical species are considered more sensitive to temperature than their counterparts from the temperate regions, because they currently live under nearly optimal temperature (Addo-Bediako et al., 2000;Deutsch et al., 2008;Chen et al., 2009), although our sampling region can be classified as temperate-zone, the climatic conditions vary dramatically among the different sites (Wu et al., 2013;Wen et al., 2016). The species of the community at low-elevation and from relatively hightemperature blocks with the greater turnover component are expected to shift uphill due to climate warming, which may alter the monotonic decline pattern to the unimodal pattern along increasing elevation. Besides, along with the increasing of species diversity at mid-elevation, extinction could be triggered by more severe interspecific competitions. Accordingly, the priority of green lacewing conservation should be given to the species from the regions with low-or mid-elevation, for example, protecting or artificially offering favorable habitats (forests for most species), and assessing the threatened species in these regions.
Cunctochrysa albolineatoides, A. barkamana, and Chrysopa sp. 1 are the dominant green lacewing species in Mts. Shaluli, and could be exploited as the potential biological control agents for the local region. Thus, more attention to conservation should be paid to these three species, especially for A. barkamana and Chrysopa sp. 1, because they may be more susceptible to climate warming as they mainly live in low-elevation and hightemperature regions.

CONCLUSION
We first analyzed the α-diversity and β-diversity of the green lacewings along elevation in a biodiversity hotspot of southwest China, combining DNA barcoding and the five molecular species delimitation methods. The distinct barcode gaps are present in 93.10% of species, which highlights the availability of DNA barcoding in regional species diversity assessment of the green lacewings. The elevational α-diversity of the green lacewings communities did not show a hump-shaped pattern similar to most of the previous studies but decreased with the increasing elevation. Nestedness replaced turnover as the main component of Sørensen dissimilarity index of the green lacewing communities with the increasing elevation, and greater nestedness occurred at low temperature blocks.
Species diversity of the green lacewings from Mts. Shaluli is strongly affected by the temperature on a broad-scale elevational gradient. Additionally, the resource and resultant interspecific competitions may also be the potential drivers on the elevational diversity pattern of the green lacewings. Priority of conservation should be given to the species at low-and mid-elevation and relatively high-temperature regions facing the challenge of climate warming.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
XL, YhL, and YL conceived the idea for this study, designed the research, and wrote the manuscript. YL provided and analyzed the data. All authors contributed to the article and approved the submitted version.