Endemism Patterns of Planthoppers (Fulgoroidea) in China

Studies on endemism are always of high interest in biogeography and contribute to better understanding of the evolution of species and making conservation plans. The present study aimed to investigate the endemism patterns of planthoppers in China by delimiting centers of endemism and areas of endemism. We collected 6,907 spatial distribution records for 860 endemic planthopper species from various resources. Centers of endemism were identified using weighted endemism values at 1° grid size. Parsimony analysis of endemicity and endemicity analysis were employed to detect areas of endemism at 1°, 1.5°, and 2° grid sizes. Six centers of endemism located in mountainous areas were identified: Taiwan Island, Hainan Island, eastern Yungui Plateau, Wuyi Mountains, western Qinling Mountains, and western Yunnan. We also delimited six areas of endemism, which were generally consistent with centers of endemism. Our findings demonstrated that mountainous areas have an essential role in facilitating the high level of endemism and formation of areas of endemism in planthoppers through the combined effects of complex topography, a long-term stable environment, and geological events. Dispersal ability and distribution of host plants also have important effects on the patterns of planthoppers’ endemism.


INTRODUCTION
The geographical distribution of endemic species (i.e., restriction of a species to a particular area) represents the highest degree of historical and ecological imprint of all biological entities (Casazza et al., 2008). Moreover, endemic areas have been generally recognized as priority areas for biodiversity conservation plans (Myers et al., 2000;Lamoreux et al., 2006;Huang et al., 2010Huang et al., , 2016Gomes-da-Silva et al., 2017;Zhao et al., 2020a). Therefore, studies related to patterns of endemism have always been a central topic of biogeography and biodiversity conservation (Laffan and Crisp, 2003;Orme et al., 2005;Posadas et al., 2006;Sandel et al., 2011;Huang et al., 2012;Feng et al., 2016;Li et al., 2017).
Studies on the spatial patterns of endemics often delimit the main distribution areas, and these areas are frequently described as centers of endemism (CoEs) and areas of endemism (AoEs). Although the drivers causing formation of CoEs and AoEs may show similar ecological or historical legacies (Hurdu et al., 2016), they differ because the two concepts are based on different identification methods and assumptions (Linder, 2001a). CoEs are areas including higher endemic richness or endemism than its surroundings (Crisp et al., 2001;Linder, 2001b). AoEs represent areas delimited by the congruent distribution of at least two endemic taxa (spatial homology) (Platnick, 1991;Morrone, 1994;Szumik et al., 2002). AoEs are extremely important because they can be used to design biogeographic regionalization schemes (Morrone, 2014a), to investigate organism-climate dynamics (Gámez et al., 2017), and as a priority target for conservation plans (Huang et al., 2010;Noroozi et al., 2018). Numerous biogeographers and evolutionary biologists have shown interest in evaluating the causes for the presence of CoEs and AoEs (Nelson and Platnick, 1981;Anderson, 1994;López-Pujol et al., 2011;Yuan et al., 2014;Noroozi et al., 2018).
Mountainous areas host a remarkable proportion of Earth's biodiversity, with several species completely restricted to these areas, although these areas only account for 16.5-27% of the land area (Fjeldså et al., 2012). Half of the high biodiversity areas identified so far are mountainous areas (Kohler and Maselli, 2009). The large biodiversity in mountainous areas is associated with their dual role as "species museums" (places of especially long-term persistence) and "species cradles" (places of especially rapid speciation). Therefore, high biodiversity in mountainous areas reflects two mechanisms: enhanced speciation rates and lineage persistence (Rahbek et al., 2019). Mountain biodiversity is characterized by deep-time evolution and ecological processes, which reflect a history worth protecting (Rahbek et al., 2019).
The planthoppers, members of superfamily Fulgoroidea (Insecta: Hemiptera), are a dominant group of herbivorous insects, consisting of 13,844 species in 33 families worldwide (Bourgoin, 2021). Planthoppers are obligatory phytophagous and mainly feed on the phloem tissue of woody or herbaceous plants (Wilson, 2005). Thus, host plants might have highly impacted their distribution. The majority of species in planthoppers lack the ability to disperse over long distances, resulting in a small distribution range. This property makes it a unique and excellent model for studying biogeography (Liang, 1998). China is one of the countries that is the richest in planthoppers, with over 1,300 described species. Unfortunately, although previous studies have investigated species richness patterns in planthoppers (Zhao et al., 2020b), the understanding of patterns of endemism in China remains in its infancy. Thus, to fully understand endemism patterns of Chinese planthoppers, there is a need to compile distributional data for all endemic planthoppers species. Research on spatial patterns of endemism not only contributes to understanding the evolution of planthoppers but also for identification of priority areas of planthopper biodiversity conservation.
The aim of the present study was to propose biogeographical patterns of planthopper endemism in China using the spatial distribution data for endemic species. We focused on the following research areas: (1) identifying the CoEs and delimiting AoEs; and (2) whether or not the CoEs and AoEs in this study are located in mountainous areas, as found in previous studies (e.g., Wang et al., 2017;Noroozi et al., 2018Noroozi et al., , 2019.

Species Distribution Data
The 6,907 spatial distribution records for 860 endemic planthoppers species (only recorded in China) were collected from the literature, books, MD/Ph.D. theses, zoological records, China Knowledge Resource Integrated Database, and museum specimens (Zoological Museum, Institute of Zoology, Chinese Academy of Sciences; Institute of Entomology, Guizhou University; the insect collection of China Agricultural University, Hebei University, Nankai University and Dali University; Entomological Museum, Northwestern A & F University; Shanghai Institute of Entomology, Chinese Academy of Sciences; Tianjin Museum of Natural History; Taiwan Agricultural Research Institute). Distribution sites at the city, county, or township levels were obtained from the original source. The distribution sites containing the latitude and longitude information were used directly, and distribution sites without the latitude and longitude information were expressed by the latitude and longitude of the corresponding administrative center.

Assessing Sampling Bias and Mapping Endemism Patterns
The inventory completeness for the study region was accessed using the species accumulation curve from the incidence-based bootstrap estimators (Zhao et al., 2020b). The presence (1) or absence (0) matrix for each endemic species in each 1 • grid (∼100 × 100 km) was constructed and was analyzed using EstimateS 9.1 with 100 randomizations (Colwell, 2013). To assess the completeness of the richness of each 1 • grid, we fitted a linear regression using the square-root transformed number of records and number of richness per grid. Spatial autocorrelation may inflate the rate of type I error (Diniz-Filho et al., 2003), so the P-value for linear regression was reported using geographically effective degrees of freedom (Dutilleul et al., 1993), evaluated using Spatial Analysis in Macroecology 4.0 (Rangel et al., 2010). The patterns of endemism were visualized by calculating the weighted endemism values of each 1 • grid. This process was conducted using Biodiverse 2.1 (Laffan et al., 2010). We defined the center of endemism, based on the criterion proposed by Crisp et al. (2001) and Linder (2001b).

Identifying Areas of Endemism
Parsimony analysis of endemicity (PAE) and endemicity analysis (EA) were used to identify AoEs. The performance of the two approaches was compared using the two criteria proposed by Escalante et al. (2009), namely the number of AoEs they delimit and the number of endemic species supporting them. Three grid sizes were used in two approaches: To perform PAE, we constructed presence (1) or absence (0) matrix based on the presence of each endemic species in each grid. A hypothetical outgroup "Root" with all zeros was added to the matrix to root the resulting tree (Morrone, 1994(Morrone, , 2014b. All matrices were performed in TNT v1.1 (Goloboff et al., 2008). New Technology algorithms were based on maximum trees with 1,000 and used sectorial search and tree fusing, with 10 initial addseqs. We obtained a strict consensus tree for each analysis. The relative support for each branch was estimated using bootstrap analysis with 1,000 replicates (Felsenstein, 1985). Branches with bootstrap values ≥ 50% were selected as candidates for AoEs. Finally, AoEs (clades of grids), containing at least two species restricted by these areas, were delimited and mapped.
EA was conducted using NDM/VNDM v 3.1 (Goloboff, 2016) under three grid sizes. This method is used to identify AoEs using an optimality criterion, which explicitly considers the spatial location of the species in a given area (Szumik et al., 2002). The parameters used here were performed by saving temporary sets within 0.99 of the current score; sets containing at least two endemic species and scores above 3.0 were saved. The search was repeated 10 times, keeping overlapping subsets when 60% of their defining species were unique. From the sets obtained, we chose species with a minimum score of 0.5. The strict rule was used to calculate the consensus areas at a cut-off of 40% similarity in species. To obtain the final AoEs, the consensus areas among the different grid sizes were overlapped (do Prado et al., 2015;Gao et al., 2018;Du et al., 2020).

RESULTS
The species accumulation curve using bootstrap estimators obtained the "true" species number, namely 1,006 (Figure 1), and the data completeness degree is 85.4%. Furthermore, the ratio of observed species richness to those predicted by linear regression models for each grid cell was >78.8% (Figure 2). These results show that planthoppers were adequately sampled.
Six CoEs were identified (Figure 3), located in Taiwan Island, Hainan Island, eastern Yungui Plateau, Wuyi Mountains, western Qinling Mountains, and western Yunnan. A total of three AoEs were identified using three grid sizes through PAE analysis (Figure 4). These areas are supported by single or multiple grid sizes. Taiwan Island and Hainan Island was  detected under all three grid sizes. Western Qinling Mountain was detected to be AoEs only at 2 • grid size. The AoEs under each grid size and the species supporting them were identified (Supplementary Table 1).
The consensus areas and supported endemic species obtained by EA in each grid size are summarized in Supplementary Material. Analysis using 1 • , 1.5 • , and 2 • grid sizes yielded 4, 10, and 9 consensus areas, respectively. Through an overlapping pattern of the consensus areas in different grid sizes, five AoEs were identified (Figure 5).

Causes of High Endemism and AoEs
In this study, we found that all CoEs and AoEs were located in mountainous areas, consistent with the findings in previous studies for insects, such as leafhoppers, scale insects, and aphids (Yuan et al., 2014;Wang et al., 2017;Gao et al., 2018), birds (Lei et al., 2003;Huang et al., 2010), mammals (Tang et al., 2006, and plants (Huang et al., 2012;Noroozi et al., 2019). The topographic complexity, long-term stable environment, and geological events experienced in mountainous areas are generally considered to be the main driving forces for the high level of endemism and formation of AoEs in these areas (Tribsch, 2004;López-Pujol et al., 2011;Wang et al., 2017;Gao et al., 2018;Noroozi et al., 2018Noroozi et al., , 2019. With increasing topographic complexity, the amount of habitat diversity is expected to increase, promoting available niche space, and thereby increasing ecological speciation via adaptation to different niches (Rundle and Nosil, 2005;Hendry et al., 2007) and the coexistence of species (Tews et al., 2004;Hortal et al., 2009). Furthermore, topographical complexity can act as a barrier to gene-flow among diverging populations, resulting in supporting reproductive isolation, increasing differentiation, and speciation (Gillespie and Roderick, 2014), and ultimately high endemism. Mountainous areas have longterm climate stability, which is conducive to the persistence of relict lineages, specialization, speciation of small-ranged species, and reduction of extinction probability (Fjeldså and Lovett, 1997;Dynesius and Jansson, 2000;Jansson, 2003;López-Pujol et al., 2011), thus further increasing diversification.
High level of endemism in mountainous areas is also significantly related to the geological events they have experienced. Orogenic processes promote the emergence of new lineages and radiation of species due to the repeated formation, connectivity, and disappearance of habitats within mountain ranges (Rahbek et al., 2019). We found that Taiwan Island and Hainan Island had the first and second highest weighted endemism values and endemicity scores, respectively. Taiwan Island and Hainan Island are China's largest and second largest islands, respectively, with typical mountain habitats affected by tropical/subtropical climates. The highest mountain in Taiwan is Mt. Yunshan (3,997 m), which is also the highest mountain in southeastern China. The highest mountain of Hainan Island is Mt. Wuzhishan. Taiwan Island and Hainan Island are isolated from the mainland, generating limited interchange between the islands and the mainland (Yuan et al., 2014). Furthermore, volcanoes often erupted on Hainan Island during the Pleistocene (Wang, 1985), causing populations to be isolated on the local and over decadal timescales, thereby increasing the development of different gene pools and the appearance of new genetic combinations (Gillespie and Roderick, 2014). These events were beneficial for promoting the emergence of novel lineages and a high degree of endemism in the two islands. The levels of endemicity in Taiwan Island are significantly higher than Hainan Island, probably because it is farther away from the mainland and had formed and separated earlier.
Planthoppers are very weak flyers and move by shortdistance jumping (Liang, 1998). Therefore, planthoppers lack a strong ability to disperse. This characteristic may allow them to remain in places of speciation and glacial refuges, particularly when these places have high altitudes and complex terrain, and thus build high endemic species numbers and AoEs over time. The vast majority of various planthoppers feed on woody or herbaceous plants (Wilson, 2005), such as Poaceae, Cyperaceae, Rosaceae, Oleaceae, and Polygonaceae (Ding, 2006;Chen et al., 2014). Host plants, apart from being a food source, are also used by some planthoppers for mating, oviposition, protection in winter, and avoiding their natural enemies (Wilson et al., 1994). Many planthoppers have high degrees of specialization with respect to host plants (Wilson et al., 1994;Liang, 1998). The family Delphacidae, which has the largest number of species among planthoppers, mainly feeds on a small number of host plants species (concentrated to Poaceae and Cyperaceae), and some even feed only on a single host plant species (Wilson et al., 1994;Ding, 2006). Most of the species in the family Caliscelidae feed on bamboo and a few on other plants (Chen et al., 2014). In China, we found that the diversity centers of many host plants are consistent with the CoEs/AoEs identified in the present study. For example, high bamboo diversity is mainly concentrated in southern Yunnan and Wuyi Mountains (Xu et al., 2019). There is high species richness of Rosaceae in northern Yunnan and the Qinling Mountains (Zou et al., 2019). This consistency shows that the distribution pattern of planthoppers is highly associated with the diversity patterns of host plants.

The Effects of Grid Size and Approach on Results
PAE and EA use grids as basic operating units in delimiting AoEs; therefore, the grid size affects the results in a deterministic manner. Grid size can affect the number, dispersion, and range of the localities. Smaller grid sizes will produce narrower AoEs, allowing a small number of endemic species to be included, whereas larger grid sizes identify wider AoEs and accommodate more endemic species (Casagranda et al., 2009;do Prado et al., 2015;DaSilva et al., 2015).
We found the greater number of AoEs and supported endemic species were identified by EA at larger grid cell (1.5 • and 2 • ), as was PAE (2 • ). EA produced five and four AoEs in 1.5 • and 2 • grids and PAE delimited three AoEs in 2 • grids. Combining different grid sizes is necessary to define the AoEs. Using several grid sizes facilitates the discovery of different AoEs in regions with topographical complexity (Elías and Aagesen, 2016), such as western Yunnan, and eastern Yungui Plateau is only delimited by at 1.5 • and 2 • grids. AoEs that recovered under different grid sizes were more robustly and clearly supported by data (Aagesen et al., 2009Szumik et al., 2012). Additionally, the use of several grids sizes was slightly different in location, area, and boundary of the resulting AoEs, and minimized the effects of sampling biases (DaSilva et al., 2015). The use of grids brings has some drawbacks, such as the inability to identify AoEs smaller than grid sizes and difficulty in representing the fuzzy edges of the AoEs. Geographical Interpolation of Endemism (GIE) (Oliveira et al., 2015), an approach independent of grids, can be valuably used in future work to obtain more detailed information about the AoEs of planthoppers. The performance of PAE and EA has been widely discussed using real or hypothetical distributions (e.g., Moline and Linder, 2006;Carine et al., 2009;Casagranda et al., 2012;Escalante, 2015). However, there is no consensus on which approach is most suitable and their exploration should be continued (Escalante, 2015). Using a total of three grids, we found that EA identified more AoEs than PAE. The observed difference is because species that support the AoEs in PAE require higher and stricter sympatry, resulting in a small number of AoEs. Furthermore, western Qinling Mountains delimited by the PAE are not identified in the EA. In this way, combining the two methods may provide a more comprehensive knowledge of AoEs.

CONCLUSION
The spatial structure of planthoppers' endemism in China was analyzed by identifying CoEs and AoEs. We found that the CoEs and AoEs identified fell in the mountainous areas, mainly due to their complex topography, stable eco-climatic conditions, and related geological events that promote long-term persistence, speciation, and species accumulation. Furthermore, we found that the dispersal ability of planthoppers and diversity pattern of host plants also are responsible for high planthoppers endemism in the mountainous areas. Future studies can add temporal information based on phylogenetic information from endemic species to better interpret AoEs.

DATA AVAILABILITY STATEMENT
The original data of this study are available from the corresponding author upon reasonable request.

AUTHOR CONTRIBUTIONS
ZZha analyzed the data and wrote the article. LNY, JL, ZC, ZZho, YZ, LAY, HL, YS, NG, and XW collected the data. XC revised the article. All authors contributed significantly to the drafts and gave final approval for publication.