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 the Entomological Museum of China Agricultural University (CAU), Beijing, China.

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).

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¹. In case of unsuccessful amplification, an alternative forward primer COI4-F (5′-TGTAAAACGACG-GCCAGTAAACTAATARCCT TCAAAG-3′)² (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₂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³. 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)⁴, 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 (Zhang et al., 2013), 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⁵. The trees were saved as Newick format in the software FigTree v.1.4.3 (Rambaut and Drummond, 2016), then imported into the PTP⁶ and mPTP⁷ 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).

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 low-elevation (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⁸ (Table 1), and the mean elevation of these counties form an elevational gradient on broad-scale. Then, we calculated the mean of these temperature data, respectively, for the final analysis. The temperature data from the sampling site 7–6 (27.7333 N, 100.7473 E) belonging to the Ninglang County was used as same as that from the sampling site 7–9 (27.7418 N, 100.7539 E) belonging to the Yanyuan County, because they are geographically near to each other (1,151 m), and the Ninglang County lack climate station.

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² + a) or (y = cx² + bx + a)] will be used to fit (Werenkraut and Ruggiero, 2014). The goodness of fit was evaluated based on the adjusted R² (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).

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.

The results from the different species delimitation methods generally correspond to the identifications of the morphospecies, but with some discrepancies (Figure 2). mPTP was the most conservative method on the species delimitation, with the fewest mOTUs (=36) recovered. Both the genetic distance-based approaches, i.e., ABGD and ASAP, resulted in more mOTUs (=37) than that from mPTP. However, GMYC and bPTP delimitated much more mOTUs (41 and 42, respectively). Specifically, each of the following four pairs of morphospecies, namely Chrysoperla furcifera (Okamoto)/Mallada sp., Nineta vittata (Wesmael)/Nineta sp. 1, Cunctochrysa sp. 1/Cunctochrysa sp. 2, and Chrysoperla annae Brooks/Chrysoperla nipponensis (Okamoto), was delimited to be the same mOTU by mPTP. Similarly, ABGD and ASAP, respectively, recognized three pairs of morphospecies [i.e., Ch. annae/Ch. nipponensis, Anachrysa sp. 1/Anachrysa sp. 2, and Chrysopidia (Chrysopidia) remanei Hölzel/Chrysopidia (Chrysopidia) sp. 1] as the same mOTU, with prior maximal distance value as 0.022. Retipenna parvula Yang et al. and Apertochrysa jiuzhaigouana (Yang et al.), respectively, were separated into two mOTUs by bPTP. Cunctochrysa albolineatoides Tsukaguchi was divided into two mOTUs by GMYC.

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).

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 ranging from 2,087 to 3,579 m; Cunctochrysa albolineatoides, accounting for 13.61% of the total number of individuals, is distributed in all the sampling counties, ranging from 2,087 to 4,202 m in elevation; Chrysopa sp. 1 accounting for 10.36% of the total number of individuals, is dominated in Batang County and Derong County with elevation ranging from 2,093 to 2,409 m.

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² = 0.66, df = 11, P < 0.01; R² = 0.53, df = 12, P < 0.01; and R² = 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).

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 low-temperature 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.

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).

The previous studies on some organisms from the Hengduan Mountains (Gong et al., 2005; Rahbek, 2005; Fu et al., 2006, 2007; Li 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 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 mid-elevation 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, Bozdoğan (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 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 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).

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 high-temperature 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 high-temperature regions.