Species Monitoring Using Unmanned Aerial Vehicle to Reveal the Ecological Role of Plateau Pika in Maintaining Vegetation Diversity on the Northeastern Qinghai-Tibetan Plateau

Plateau pika (Ochotona curzoniae, hereafter pika) is considered to exert a profound impact on vegetation species diversity of alpine grasslands. Great efforts have been made at mound or quadrat scales; nevertheless, there is still controversy about the effect of pika. It is vital to monitor vegetation species composition in natural heterogeneous ecosystems at a large scale to accurately evaluate the real role of pika. In this study, we performed field survey at 55 alpine grassland sites across the Shule River Basin using combined methods of aerial photographing using an unmanned aerial vehicle (UAV) and traditional ground measurement. Based on our UAV operation system, Fragmentation Monitoring and Analysis with aerial Photography (FragMAP), aerial images were acquired. Plot-scale vegetation species were visually identified, and total pika burrow exits were automatically retrieved using the self-developed image processing software. We found that there were significant linear relationships between the vegetation species diversity indexes obtained by these two methods. Additionally, the total number of identified species by the UAV method was 71, which was higher than the Quadrat method recognition, with the quantity of 63. Our results indicate that the UAV was suitable for long-term repeated monitoring vegetation species composition of multiple alpine grasslands at plot scale. With the merits of UAV, it confirmed that pika’s disturbance belonged to the medium level, with the density ranging from 30.17 to 65.53 ha−1. Under this density level, pika had a positive effect on vegetation species diversity, particularly for the species richness of sedge and forb. These findings conclude that the UAV was an efficient and economic tool for species monitoring to reveal the role of pika in the alpine grasslands.


Study Area
During the period from the end of July to the middle of August in 2019, field sampling and aerial photographing were performed in a transect that traversed the typical vegetation within the source region of Shule River Basin, which is situated at the northeastern edge of the QTP, China (Figure 1). The study area belongs to a typical semi-arid region with the mean annual temperature of −4°C and mean annual precipitation ranging from 200 to 400 mm. Notably, more than 90% of precipitation falls during the summer monsoon season. There are six permafrost types, including extreme stable permafrost, stable permafrost, sub-stable permafrost, transition permafrost, unstable permafrost and seasonal frost, accounting for approximately 18%, 16%, 18%, 13%, 20% and 25% of the total area of the source region of Shule River Basin ( Table 1). The average annual ground temperature of these six types of permafrost ranges from > 0.5 to < −5.0°C (Table 1). Four typical alpine grasslands, i.e., alpine swamp meadow (ASwM), alpine meadow (AM), alpine steppe meadow (AStM) and alpine steppe (AS) are widely distributed in our study area. The vegetation species composition of each type of alpine grassland can be found in Table A1.  We set up 55 sites among the gradients at an interval of 1 to 5 km dependent on the road availability for walking and vehicle ( Figure 1). There were 7 sites in ASwM, 14 sites in AM, 10 sites in AStM and 24 sites in AS, respectively. Three for ASwM, 5 for AM, 4 for AStM and 9 for AS were selected randomly to perform aerial photographing (UAV-based method) and field sampling (Quadrat method) simultaneously, and the other 34 sites were performed using the UAV-based method.

Study Scheme
The flowchart of the study scheme is shown in Figure 2. It contains five steps: (1) take aerial images by the UAV at 2 and 20 m height, respectively; (2) retrieve plot-scale total pika burrow exits based on the images that were taken at 20 m height and identify vegetation species from each image that was taken at 2 m height; (3) calculate the coefficients of active pika burrows and pika density, and test the feasibility of the UAV for vegetation species monitoring; (4) acquire plot-scale pika density and vegetation species diversity indexes and (5) explore the relationship between plot-scale pika density and vegetation species diversity indexes to reveal the role of pika. The detailed information about aerial photographing, field sampling and aerial image analysis is listed in the following sections.
Remote Sens. 2020, 12 We set up 55 sites among the gradients at an interval of 1 to 5 km dependent on the road availability for walking and vehicle ( Figure 1). There were 7 sites in ASwM, 14 sites in AM, 10 sites in AStM and 24 sites in AS, respectively. Three for ASwM, 5 for AM, 4 for AStM and 9 for AS were selected randomly to perform aerial photographing (UAV-based method) and field sampling (Quadrat method) simultaneously, and the other 34 sites were performed using the UAV-based method.

Study Scheme
The flowchart of the study scheme is shown in Figure 2. It contains five steps: (1) take aerial images by the UAV at 2 and 20 m height, respectively; (2) retrieve plot-scale total pika burrow exits based on the images that were taken at 20 m height and identify vegetation species from each image that was taken at 2 m height; (3) calculate the coefficients of active pika burrows and pika density, and test the feasibility of the UAV for vegetation species monitoring; (4) acquire plot-scale pika density and vegetation species diversity indexes and (5) explore the relationship between plot-scale pika density and vegetation species diversity indexes to reveal the role of pika. The detailed information about aerial photographing, field sampling and aerial image analysis is listed in the following sections.

UAV and Aerial Photographing
A DJI drone (Mavic 2 Zoom, DJI Innovation Company, China) is used in this study. It is a lightweight four-wheel drone equipped with an omnidirectional sensing system and autopilot system. Vertical and horizontal accuracies of this drone are ± 0.5 and ± 1.5 m, respectively. The maximum flying height could reach 6000 m above sea level. The duration of flight is 31 minutes (at a constant speed of 25 km h -1 ) and the hovering endurance is 29 minutes. A lightweight sensor is integrated in the front body of Mavic 2 Zoom. It is a 24-48 mm digital optical zoom camera with a resolution of 12 Mega pixels (maximum image size is 3000 × 4000) and a three-axis gimbal. The field vision is 83° at 24 mm focus and 48° at 48 mm focus. The sensor can receive the irradiance of three  Figure 2. The flowchart of study scheme. N APB , N TPB and N PD mean the number of active pika burrows, total pika burrows and pika density in a hectare, respectively.

UAV and Aerial Photographing
A DJI drone (Mavic 2 Zoom, DJI Innovation Company, China) is used in this study. It is a lightweight four-wheel drone equipped with an omnidirectional sensing system and autopilot system. Vertical and horizontal accuracies of this drone are ±0.5 and ±1.5 m, respectively. The maximum flying height could reach 6000 m above sea level. The duration of flight is 31 min (at a constant speed of 25 km h −1 ) and the hovering endurance is 29 min. A lightweight sensor is integrated in the front body of Mavic 2 Zoom. It is a 24-48 mm digital optical zoom camera with a resolution of 12 Mega pixels (maximum image size is 3000 × 4000) and a three-axis gimbal. The field vision is 83 • at 24 mm focus and 48 • at 48 mm focus. The sensor can receive the irradiance of three spectral bands (red, green and Remote Sens. 2020, 12, 2480 5 of 21 blue) and they are stored in JPEG format with digital number from 0 to 255. Information about the longitude and latitude of the central position was recorded in each aerial image.
Vegetation species were identified from traditional quadrat and aerial images that were taken at a height of 2 m using UAV to evaluate plot-scale vegetation diversity. Total pika burrow exits were retrieved from each aerial image that was taken at a height of 20 m by UAV to estimate plot-scale pika density. Briefly, on each site, we set a GRID flight way (200 m × 200 m) (Figure 3a). In each GRID flight way, we set up 3 BELT flight ways (40 m × 40 m) (Figure 3b). The UAV was auto-piloted to 16 preset way points to take photos at a height of 20 (GRID flight way) and 2 m (BELT flight way) with the camera looking vertically down using Fragmentation Monitoring and Analysis with aerial Photography (FragMAP) [41]. The flight speeds were preset to 8 and 4 m s −1 in the GRID flight way and BELT flight way, respectively. The UAV hovered for 3 s in each way point and took photographs within 1 s in both the GRID and BELT flight way. Overall, 880 aerial images were taken to extract total pika burrow exits, and another 2640 images were taken to identify vegetation species visually.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 22 spectral bands (red, green and blue) and they are stored in JPEG format with digital number from 0 to 255. Information about the longitude and latitude of the central position was recorded in each aerial image. Vegetation species were identified from traditional quadrat and aerial images that were taken at a height of 2 m using UAV to evaluate plot-scale vegetation diversity. Total pika burrow exits were retrieved from each aerial image that was taken at a height of 20 m by UAV to estimate plot-scale pika density. Briefly, on each site, we set a GRID flight way (200 m × 200 m) (Figure 3a). In each GRID flight way, we set up 3 BELT flight ways (40 m × 40 m) (Figure 3b). The UAV was auto-piloted to 16 preset way points to take photos at a height of 20 (GRID flight way) and 2 m (BELT flight way) with the camera looking vertically down using Fragmentation Monitoring and Analysis with aerial Photography (FragMAP) [41]. The flight speeds were preset to 8 and 4 m s -1 in the GRID flight way and BELT flight way, respectively. The UAV hovered for 3 seconds in each way point and took photographs within 1 second in both the GRID and BELT flight way. Overall, 880 aerial images were taken to extract total pika burrow exits, and another 2640 images were taken to identify vegetation species visually.

Synchronous Survey of UAV and Ground Measurement
In each BELT flight way of the abovementioned 21 sites, four 0.5 m × 0.5 m quadrats were placed in the way points 6, 7, 10 and 11 before the UAV arrives ( Figure 3c). After aerial photographing, all of the vegetation species occurring within each quadrat were identified visually and classified into four functional groups, i.e., grass, legume, sedge and forb. Then, each species was clipped at the ground level, dried (oven-dried at 65 °C for 48 h) and weighed in laboratory. The dry matter weight of each species was further used to calculate the biomass of the four functional groups.
Active pika burrows were determined by covering all pika burrows with soil and then counted the open burrows after 72 h. Pikas were live-captured by string nooses anchored in all active pika burrows. The detailed information for investigating active pika burrows and pika density can be found in our previous study [35].

Aerial Image Preprocessing
Prior to aerial image analysis, Photoshop Camera Raw (Adobe Systems Incorporated, USA) was used to correct geometric distortion. For images that were taken from the air at 20 m, total pika

Synchronous Survey of UAV and Ground Measurement
In each BELT flight way of the abovementioned 21 sites, four 0.5 m × 0.5 m quadrats were placed in the way points 6, 7, 10 and 11 before the UAV arrives ( Figure 3c). After aerial photographing, all of the vegetation species occurring within each quadrat were identified visually and classified into four functional groups, i.e., grass, legume, sedge and forb. Then, each species was clipped at the ground level, dried (oven-dried at 65 • C for 48 h) and weighed in laboratory. The dry matter weight of each species was further used to calculate the biomass of the four functional groups.
Active pika burrows were determined by covering all pika burrows with soil and then counted the open burrows after 72 h. Pikas were live-captured by string nooses anchored in all active pika burrows. The detailed information for investigating active pika burrows and pika density can be found in our previous study [35].

Aerial Image Preprocessing
Prior to aerial image analysis, Photoshop Camera Raw (Adobe Systems Incorporated, San Jose, CA, USA) was used to correct geometric distortion. For images that were taken from the air at 20 m, total pika burrow exits, fractional vegetation cover (FVC) and bare patch area fraction were retrieved using the image processing software Proposal-Classifier and Pixel-Classifier [41]. Total pika burrow exits were first identified by automatic recognition and then were verified and corrected by manpower. The threshold method based on excess green index (EGI = 2G-R-B; with R, G, B being red, green, and blue bands, respectively) of each pixel was applied to calculate FVC and bare patch area fraction from original digital picture without georeferencing [35]. The specific procedures were to: (1) calculate an initial value of EGI threshold and compare it with each pixel; (2) classify each pixel in a photograph as either a green vegetation pixel or a bare pixel based on this threshold; (3) draw the contours based on the locations of green vegetation pixels using the OpenCV package; (4) overlay on and compare the contours with the original photograph; and (5) if the contours did not fit with the shape of the vegetation patches, then the threshold was adjusted and the process repeated from procedure 2 until the contours fit well with the vegetation patches in the original photograph. All the species occurring within each aerial image that was taken at 2 m were also identified visually and classified into four functional groups as the quadrat.
We were aware that there is a little bit of distortion around the margin of aerial images, even though we took them vertically down. However, it is reasonable to extract data directly (no referencing and rectification) from the aerial images as they were collected along the same method using the same type UAV, and it did not affect the processes of identifying species and pika burrows based on former studies [35,37]. In addition, the GRID and BELT modes were set to match traditional ecological study areas, without overlap among aerial images [41]. Furthermore, the cover area of each aerial image is small (26 m × 35 m and 2.5 m × 3.5 m), meaning that it is hard and inefficient to rectify each one by setting ground reference points, but it meets the study accuracy requirements [28,35,37].

Pika Density at Plot Scale
Two ratios were used to estimate the active pika burrows and pika density at plot scale [35]. The first one (R 1 ) was the ratio between the number of active pika burrows and total pika burrows, and the second (R 2 ) was the ratio between the number of pika caught and the number of active pika burrow. Both of these two ratios were calculated based on the field data for each grassland type. Plot-scale active pika burrows and pika density were then estimated by the Equations (1) and (2).
where N APB , N TPB and N PD are the number of active pika burrows, total pika burrows and pika density in a hectare, respectively.

Calculation of Vegetation Species Diversity Indexes
To test the feasibility of the UAV method for monitoring vegetation species composition, four common species diversity indexes-species richness (N), Shannon index (H), Simpson index (D) and Pielou's J index (E)-were selected and used to compare with the traditional Quadrat method. All of these indexes were calculated based on the data of vegetation species identified from aerial images and ground quadrats. Species richness (N) was the total number of species occurring within each sampling unit, i.e., 4 ground measured quadrats or the coverage area of 16 aerial images.
where pi is the frequency proportion of species i in 4 quadrats or 16 photographs of each sampling unit. Pielou' s J index (E) was calculated by Equation (5) [44].
where H and N are the Shannon index and species richness, respectively.

Data Statistics Analysis
One-way analysis of variance (ANOVA) was used to determine statistical differences of pika density, vegetation diversity indexes and biomass of functional groups at p = 0.05 level. Firstly, we tested the normal distribution and homogeneity of variances using the Shapiro-Wilk test and Bartlett test. Based on the P-value of these two tests, the Tukey HSD model or Kruskal-Wallis H model were then selected for multiple comparisons [45]. Regression models were selected dependent on the comparative results by curve estimation. Both the linear regression model and logarithmic regression model passed the significant test at 0.01 level. However, the linear regression model performed well for the correlations of vegetation diversity indexes estimated by two methods, and the logarithmic regression model expressed the relations between vegetation diversity indexes and pika density better. Therefore, the relationships between vegetation diversity indexes estimated by the Quadrat method and UAV method were analyzed by the linear regression model (y = ax + b), while the relationships between vegetation diversity indexes and pika density were analyzed by the logarithmic regression model (y = aln(x) + b). The correlation of pika density and species richness of functional groups was further studied using CANOCO software version 5 (Biometrics). A detrended correspondence analysis (DCA) showed that the values of the extracted gradient lengths in the data were less than 3 SD and thus redundancy analysis (RDA) was used [46]. All statistical analyses were performed using R version 3.5.3 (R core team, 2018) [47] and CANOCO software version 5 (Biometrics) [46].

Monitoring Vegetation Species Composition and Pika Density at Plot Scale
The aerial image that was taken by the UAV at 2 m flight height covered ground area of 2.6 m × 3.5 m and the resolution of each single image was~1 mm. Plot-scale vegetation species composition could be acquired based on the aerial images of one route. The single aerial image that was taken at 20 m flight height, with the coverage of 26 m × 35 m, had the resolution of approximate 1 cm. The total number of pika burrow exits was easily identified and retrieved to estimate plot-scale pika density.

Accuracy of Vegetation Species Composition Monitoring based on UAV
Four common species diversity indexes based on the data of vegetation species identified from aerial images and ground quadrats were used to evaluate the monitoring accuracy of the UAV-based method. The linear regression model showed that there were significant relationships between vegetation diversity indexes (species richness, Shannon index, Simpson index and Pielou's J index) estimated by the Quadrat method and UAV-based method (p < 0.01; Figure 4). At each sampling unit, the average numbers of identified species were 22 and 21 in the alpine swamp meadow, 22 and 25 in the alpine meadow, 22 and 28 in the alpine steppe meadow, 16 and 19 in the alpine steppe, as measured by the Quadrat method and UAV method, respectively. Generally, the UAV method monitored more species than the Quadrat method in three out of the four grassland types ( Figure 5), although the difference was significant only in the alpine steppe (p < 0.05). Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 22

Vegetation Species Composition at Plot Scale
At the plot scale, differences between vegetation diversity indexes calculated by the UAV method were not significant among different types of grasslands (Table 2). Species richness, Shannon index, Simpson index and Pielou's J index ranged from 13.07 ± 5.20 to 16.75 ± 5.39, 2.15 ± 0.43 to 2.45 ± 0.35, 0.75 ± 0.05 to 0.78 ± 0.04 and 0.29 ± 0.03 to 0.32 ± 0.05 in the four types of alpine grasslands, respectively. The highest species richness, Shannon index and Simpson index were all found in the alpine steppe meadow, and the lowest were in the alpine steppe (Figure 6a-c). Nevertheless, the opposite trend was found for Pielou's J index (Figure 6d). Generally, the UAV method could monitor more species than the Quadrat method. A total of 63 and 71 species were identified by the Quadrat and UAV method (Table A1), respectively.

Vegetation Species Composition at Plot Scale
At the plot scale, differences between vegetation diversity indexes calculated by the UAV method were not significant among different types of grasslands (Table 2). Species richness, Shannon index, Simpson index and Pielou's J index ranged from 13.07 ± 5.20 to 16.75 ± 5.39, 2.15 ± 0.43 to 2.45 ± 0.35, 0.75 ± 0.05 to 0.78 ± 0.04 and 0.29 ± 0.03 to 0.32 ± 0.05 in the four types of alpine grasslands, respectively. The highest species richness, Shannon index and Simpson index were all found in the alpine steppe meadow, and the lowest were in the alpine steppe (Figure 6a-c). Nevertheless, the opposite trend was found for Pielou's J index (Figure 6d). Generally, the UAV method could monitor more species than the Quadrat method. A total of 63 and 71 species were identified by the Quadrat and UAV method (Table A1), respectively.

Pika Burrows and Density at Plot Scale
There were large variations of pika burrows and density even within one type of grassland (Figure 7), despite the fact that the significant difference was only found in total pika burrows among different grasslands (Table 2). Total pika burrows, active pika burrows and pika density were in a range from 12.36 to 1260 ha −1 , 4.50 to 415.36 ha −1 , and 0.81 to 160.78 ha −1 in all grassland types. However, the mean total pika burrows, active pika burrows and pika density were all the largest in the alpine meadow, with values of 462.02 ha −1 (Figure 7a), 152.39 ha −1 (Figure 7b) and 65.53 ha −1 (Figure 7c), respectively, and they were 200 to 300 ha −1 for the mean total pika burrows (Figure 7a), 96 to 128 ha −1 for the mean active pika burrows (Figure 7b) and 30 to 44 ha −1 for the mean pika density (Figure 7c) in the other types of grasslands. There were large variations of pika burrows and density even within one type of grassland (Figure 7), despite the fact that the significant difference was only found in total pika burrows among different grasslands (Table 2). Total pika burrows, active pika burrows and pika density were in a range from 12.36 to 1260 ha -1 , 4.50 to 415.36 ha -1 , and 0.81 to 160.78 ha -1 in all grassland types. However, the mean total pika burrows, active pika burrows and pika density were all the largest in the alpine meadow, with values of 462.02 ha -1 (Figure 7a), 152.39 ha -1 (Figure 7b) and 65.53 ha -1 ( Figure  7c), respectively , and they were 200 to 300 ha -1 for the mean total pika burrows (Figure 7a), 96 to 128 ha -1 for the mean active pika burrows (Figure 7b) and 30 to 44 ha -1 for the mean pika density ( Figure  7c) in the other types of grasslands. The blue bold solid horizontal lines are the median value. ASwM, AM, AStM and AS indicate alpine swamp meadow, alpine meadow, alpine steppe meadow and alpine steppe, respectively.

Fractional Vegetation Cover, Bare Patch Area Fraction and Biomass of Functional Groups
Fractional vegetation cover was 49.39 ± 17.17% in the alpine swamp meadow, 46.70 ± 19.79% in the alpine meadow, 47.22 ± 16.04% in the alpine steppe meadow and 35.06 ± 15.14% in the alpine steppe, respectively (Figure 8a). Bare patch area fraction was nearly 65% in the alpine steppe, while it was approximately 50% for the other types of grasslands (Figure 8b). Although the difference in fractional vegetation cover and bare patch area fraction among different types of grasslands was not significant (P > 0.05; Table 2), the alpine steppe meadow had the higher fractional vegetation cover and fewer bare patches than the other grasslands ( Figure 8). Vegetation biomass of functional groups showed sharp variation among the different types of grasslands (Table 3 and Figure 9). Both the biomass of sedge (59.55-61.16 kg m −2 ) and forb (39.23-74.84 kg m −2 ) was higher than legume (1.72-9.91 kg m −2 ) and grass (1.36-15.90 kg m −2 ) in the alpine swamp meadow and alpine meadow (P > 0.05; Figure 9a

Fractional Vegetation Cover, Bare Patch Area Fraction and Biomass of Functional Groups
Fractional vegetation cover was 49.39 ± 17.17% in the alpine swamp meadow, 46.70 ± 19.79% in the alpine meadow, 47.22 ± 16.04% in the alpine steppe meadow and 35.06 ± 15.14% in the alpine steppe, respectively (Figure 8a). Bare patch area fraction was nearly 65% in the alpine steppe, while it was approximately 50% for the other types of grasslands (Figure 8b). Although the difference in fractional vegetation cover and bare patch area fraction among different types of grasslands was not significant (p > 0.05; Table 2), the alpine steppe meadow had the higher fractional vegetation cover and fewer bare patches than the other grasslands ( Figure 8). Vegetation biomass of functional groups showed sharp variation among the different types of grasslands (Table 3 and Figure 9).         The bold solid horizontal line means the median value. Grass is monocotyledon, which is characterized by an alternate leaf, stem section and caryopsis. Legume, belonging to the family of Rosalesis, is dicotyledon, which can fix nitrogen by symbiotic rhizobia. Sedge is also a monocotyledon with a solid stem, leaf sheath and small flower clustered into a spikelet. Forb is the general name of all herbage plants other than grass, legume and sedge.

Relationships between Pika Density and Vegetation Diversity Indexes
Vegetation species richness, Shannon index and Simpson index had a positive correlation with pika density (p < 0.01; Figure 10a-c), whereas Pielou's J index had a negative correlation with pika density (p < 0.01; Figure 10d). For the richness of functional groups, the most grass and legume species were in the alpine steppe, while the alpine meadow had the maximum species of sedge and forb ( Figure 11). Noticeably, sedge and forb had a higher correlation with pika density than grass and legume (Figure 11). Monte Carlo permutation tests showed that total pika burrows, active pika burrows and pika density separately explained 13%, 12.9% and 11.3% of the variation of the diversity of functional groups; however, they accounted for 20.5% of the variation of the diversity of functional groups in total (Table 4). Although pika density explained only 13% of the variation, both the simple and conditional effects showed that pika density had significant influence on the vegetation species richness of functional groups (p < 0.01; Table 4).

Relationships between Pika Density and Vegetation Diversity Indexes
Vegetation species richness, Shannon index and Simpson index had a positive correlation with pika density (P<0.01; Figure 10a-c), whereas Pielou's J index had a negative correlation with pika density (P<0.01; Figure 10d). For the richness of functional groups, the most grass and legume species were in the alpine steppe, while the alpine meadow had the maximum species of sedge and forb ( Figure 11). Noticeably, sedge and forb had a higher correlation with pika density than grass and legume ( Figure 11). Monte Carlo permutation tests showed that total pika burrows, active pika burrows and pika density separately explained 13%, 12.9% and 11.3% of the variation of the diversity of functional groups; however, they accounted for 20.5% of the variation of the diversity of functional groups in total (Table 4). Although pika density explained only 13% of the variation, both the simple and conditional effects showed that pika density had significant influence on the vegetation species richness of functional groups (P<0.01; Table 4).

Monitoring Vegetation Species Composition at Plot Scale
Long-term repeated vegetation species composition monitoring at a large scale is a foundation for improving our knowledge of the role of plant diversity in maintaining the stability and multifunctionality of ecosystems. Nevertheless, study of monitoring vegetation species composition regarding alpine grasslands is scarce, especially at the plot scale with a mass of sampling sites because of the harsh conditions and high altitude [48]. In this study, we test a practical method for monitoring vegetation species composition by the UAV at plot scale. Our results clearly show that there are significant linear relationships of vegetation diversity indexes between the Quadrat method and the UAV method (Figure 4), which indicates that the UAV method performs well in monitoring vegetation species composition compared with the traditional field observation (Figure 12a-e). Due to the high canopy density (Figure 8a), the understory is difficult to find from aerial images and thus the UAV method monitors fewer species than Quadrat method at each sampling unit in the alpine swamp meadow. However, each plot has 16 photographs, which is much more than the traditional 3-5 quadrats. The total number of identified vegetation species by the UAV method is higher than the Quadrat method in all types of grasslands at the plot scale (Table A1). We therefore conclude that the UAV method can be used to monitor the vegetation species composition for multiple types of alpine grasslands at plot scale.
The UAV method overcomes the limitation of the excessive disturbance by traditional field survey (Figure 12e, f) and the difficulty in identifying individual species by satellite remote sensing Figure 11. Redundancy analysis (RDA) triplot presenting the relationship of vegetation species richness of functional groups with pika burrows and pika density. The blue solid arrows are vegetation species richness of functional groups. The hollow red arrows are pika burrows and pika density. Green, orange, yellow and blue solid cicles are all plots for vegetation species richness of functional groups investigating in the alpine swamp meadow, alpine meadow, alpine steppe meadow and alpine steppe, respectively. TPB, APB and PD are total pika burrows, active pika burrows and pika density, respectively.

Monitoring Vegetation Species Composition at Plot Scale
Long-term repeated vegetation species composition monitoring at a large scale is a foundation for improving our knowledge of the role of plant diversity in maintaining the stability and multi-functionality of ecosystems. Nevertheless, study of monitoring vegetation species composition regarding alpine grasslands is scarce, especially at the plot scale with a mass of sampling sites because of the harsh conditions and high altitude [48]. In this study, we test a practical method for monitoring vegetation species composition by the UAV at plot scale. Our results clearly show that there are significant linear relationships of vegetation diversity indexes between the Quadrat method and the UAV method (Figure 4), which indicates that the UAV method performs well in monitoring vegetation species composition compared with the traditional field observation (Figure 12a-e). Due to the high canopy density (Figure 8a), the understory is difficult to find from aerial images and thus the UAV method monitors fewer species than Quadrat method at each sampling unit in the alpine swamp meadow. However, each plot has 16 photographs, which is much more than the traditional 3-5 quadrats. The total number of identified vegetation species by the UAV method is higher than the Quadrat method in all types of grasslands at the plot scale (Table A1). We therefore conclude that the UAV method can be used to monitor the vegetation species composition for multiple types of alpine grasslands at plot scale. and thus the UAV method can survey a larger area. Secondly, the UAV method is efficient [41]. The evidence from field work demonstrated that the UAV method needs merely one tenth of the time to monitor vegetation composition in the same sampling unit than the traditional in situ observation [37]. Thirdly, the UAV method is non-destructive and easy to perform. In our field investigation, less than 30 minutes is expended to fly over triple 16 preset way points in a 200 m × 200 m plot and take photos automatically. Additionally, the UAV flying depending on the FragMAP is suitable for the long-term fixed monitoring of vegetation species composition [41].

Effects of Pika's Disturbance on Vegetation Species Diversity
Our results show that pika's disturbance increases vegetation species diversity, which may be attributed to its clipping and foraging activities improving the species of sedge and forb. Substantial studies in other similar burrowing animals also found that their disturbances enhance vegetation The UAV method overcomes the limitation of the excessive disturbance by traditional field survey (Figure 12e,f) and the difficulty in identifying individual species by satellite remote sensing images. Firstly, the high resolution of aerial images that were taken by the UAV is capable of identifying individual species versus satellite remote sensing images. To our knowledge, the high-resolution of satellite remote sensing images can reach up to tens of centimeters [49]. The dominant vegetation species of alpine grasslands are typically small in size and highly mixed. They are difficult to distinguish, even by the high-resolution of satellite remote sensing images. Although field survey gives a guarantee to investigate each species occurred in the sampling unit, it is hard to perform the long-term fixed monitoring at a large scale. Compared with the traditional quadrat sampling range (1 m × 1 m), the UAV can gather large amounts of information on the distribution of vegetation. For example, each aerial image that was taken at 2 m height covers approximately 2.6 m × 3.5 m of ground and thus the UAV method can survey a larger area. Secondly, the UAV method is efficient [41]. The evidence from field work demonstrated that the UAV method needs merely one tenth of the time to monitor vegetation composition in the same sampling unit than the traditional in situ observation [37]. Thirdly, the UAV method is non-destructive and easy to perform. In our field investigation, less than 30 min is expended to fly over triple 16 preset way points in a 200 m × 200 m plot and take photos automatically. Additionally, the UAV flying depending on the FragMAP is suitable for the long-term fixed monitoring of vegetation species composition [41].

Effects of Pika's Disturbance on Vegetation Species Diversity
Our results show that pika's disturbance increases vegetation species diversity, which may be attributed to its clipping and foraging activities improving the species of sedge and forb. Substantial studies in other similar burrowing animals also found that their disturbances enhance vegetation species richness [17,[50][51][52], such as prairie dog [53] and Daurian pika [54]. The most grazed forage by pika are grass and legume, follows by forb and sedge [55], and our results also fit this pattern (Figure 9). Low and creeping vegetation under the canopy is generally dependent on the resource availability [56]. Canopy interception decreases the water and light availability of understory [57,58], and thus slows down the decomposition rates of litter and soil organic matter [59]. Both of these factors inhibit the growth and establishment of vegetation and hence reduce species richness [60]. Pika frequently clips tall plants to avoid the risk of predation [61]. Meanwhile, its foraging can decrease the stand withered leaves and remove apical dominance of vegetation [17,62]. Both clipping and foraging activities open the canopy and thus enhance the infiltration of water and light [63,64], and meanwhile pika's burrowing promotes soil nutrition circulation [65,66]. These processes potentially alleviate understory resource limitation [67] and promote the species richness of forb [68]. The high abundance of forb and the great floral density on habitats can in turn enhance vegetation species diversity by increasing the richness of pollinating insects [52].
Ecosystems with greater species diversity can maintain higher community productivity, stability and ability to resist invasive species [69][70][71]. Our results support the idea that medium disturbance by pika leads to the highest levels of vegetation species diversity [72,73], indicating that medium pika disturbance intensity may be beneficial to maintain the stability and the multi-functionality of alpine grasslands. The average pika density of four types of alpine grasslands is in a range of 30.17 to 65.53 ha −1 . According to the pika density classification, pika density belongs to the medium level in this study [74]. Under this density level, vegetation species richness and diversity indexes have a significant positive relationship with pika density (Figure 10). This finding is different from those of previous studies [24,63,75], which suggest that unimodal curve function is found between species diversity indexes and pika density. Unlike these studies used few plots at small scale (10 m × 10 m~25 m × 25 m), our field survey is carried out at large scale (200 m × 200 m) with 55 plots. The larger scale and more plots reduce sampling bias and uncertainty in heterogeneous distribution of pika and vegetation species composition [35,37], and thus we believe that our results are robust to accurately estimate the real relationship between pika density and vegetation species diversity.

Limitation of the Current Monitoring Method by the UAV and Future Work
To our knowledge, this is the first attempt to monitor vegetation species composition at large scale with a mass of plots and analyzes the relationship between pika density and vegetation diversity indexes in real-world alpine grasslands. However, we do acknowledge there is still room for improvement in this study. First of all, the accuracy of species identification needs to be improved. Three species, i.e., Gentiana squarrosa, Potentilla acaulis and Sibbaldia procumbens failed to be identified from UAV photographs (Table A1). Gentiana squarrosa and Potentilla acaulis are rare species (only found in 3 plots) with the biomass being less than 1%, and Sibbaldia procumbens is auxiliary species (found in 12 plots) with the biomass being approximately 3%. Given that both the Shannon index and Simpson index have significant linear relationships between the Quadrat method and the UAV method (p < 0.01, Figure 4b,c), we therefore conclude that the error derived from unidentified species is limited [37]. Although the error for species identification by the UAV is weak, observation in situ still has to be considered to investigate the species composition in grasslands with dense canopy such as the alpine swamp meadow (Figure 8a) for the fine study at plot scale. To date, vegetation species recognition has depended on visual identifying by manpower. It requires substantial time for aerial images analysis. The application of automate species recognition has been in view [37], which will further improve the efficiency of the UAV.
Secondly, the influence of pika's disturbance on seed dispersal and vegetation succession remains to be developed. Other than foraging, pika's burrowing not only intensifies the spatial heterogeneity of the land surface [76], but also creates open gaps for vegetation invasion [77,78]. Additionally, it is considered that pika directly disperses seeds via carrying behavior and feces [63], or indirectly disperses them through the birds living in discarded burrows [79]. However, there is no direct evidence regarding the influence of pika's disturbance on seed dispersal and the vegetation succession process, especially on the pika burrow and the associated bare patch. The long-term repeated monitoring, therefore, is still needed to (1) explore the change of pika density and the corresponding response of plant community structure, species and functional diversity; (2) investigate the seed storage in pika pile and bare patch originated from the pika pile; and (3) reveal the potential vegetation succession process on the open gaps.

Conclusions
The aerial image that was taken by the UAV at 2 m flight height covers a ground area of approximately 2.6 m × 3.5 m with the resolution of millimeter-level, which overcomes the limitation of the difficulty in identifying individual species by satellite remote sensing images. Additionally, no more than 30 min is used to fly over triple 16 preset way points in a 40,000 m 2 plot and take photos automatically, which can effectually avoid the excessive disturbance caused by traditional field survey. At the plot scale, we find that UAV method monitors more species than the Quadrat method, indicating that UAV can be used to monitor the vegetation species composition for multiple types of alpine grasslands at a large scale. The coverage range of a single aerial image at 20 m flight height approximately matches the satellite image pixel (30 m × 30 m resolution). The resolution of a single image at this height is~1 cm. Pika burrow exits, with a diameter of~6 cm, can be easily identified from the aerial image and used for estimating plot-scale pika density. Based on the pika density at plot scale, our result confirms that pika's disturbance belongs to the medium level. Under this density level, pika's disturbance promotes vegetation species biodiversity by increasing the species of sedge and forb. Given that the influence mechanisms of pika's disturbance on seed dispersal, vegetation succession and functional diversity remain unclear, long-term repeated monitoring is still required, and UAV is an efficient and economic tool to fulfill the requirement.     Total species did not indentify 11 3