Dongting Lake (DL) (28°30′–29°38′N, 112°18′–113°15′E) is located in the middle reaches of the Yangtze River and in the north of Hunan Province (Li et al., 2013). The existing natural lake area of DL is about 2,625 km² (Ding and Li, 2011; Hu et al., 2020), which can be divided into three parts: East DL, South DL and West DL. The lake body is connected by numerous river networks and flood channels, and the bottom of the lake inclines slightly from northwest to Southeast (Ding and Li, 2011). DL belongs to the subtropical monsoon climate. The annual average temperature is 16.4–17°C, the frost-free period is 258–275 days, and the annual precipitation is 1,100–1,400 mm (Yang et al., 2009). The rainfall from april to June accounts for more than 50% of the total annual precipitation. Therefore, the water level of DL begins to rise in april, the flood season is from May to October, and the dry season is from November to april of the next year (Li et al., 2011). The huge drop of water level (up to 18 m) in wet and dry seasons produces unique natural conditions giving birth to diverse landscapes and rich natural resources, providing habitat for a large number of animals and plants (Guan et al., 2016).

In this study, in order to better explore the impact of PER on the habitat suitability of birds in DL, we integrated the traditional DL regional division method and the recent establishment and adjustment of nature reserve of DL, and divided DL into four regions, including East DL (E-DL), West DL (W-DL), South DL (S-DL) and Hengling Lake (HL). E-DL and W-DL are national nature reserves, S-DL and HL are provincial nature reserves (Figure 1).

Remote sensing is the most important data source in this work. In view of the need to splice two images covering the study area, this study selected 12 scenes (Supplementary Table S1) with the resolution being 30 m of 6 years from 2014-2019, which were collected from the USGS website (https://earthexplorer.usgs.gov/). We give priority to landsat8 (OLI) images, but our images also include four Landsat7 (ETM+) images because of the limitation of cloud cover, and the acquisition time span of these images is from the end of November to the beginning of February, but they are still in the winter dry season. All images are L1TP level data thus don’t need geometric correction (USGS, 2019).

In this study, eCognition Developer 9.01 (Definiens Imaging GmbH, Munich) based on object-oriented classification was used for land cover interpretation (Tian et al., 2008b). Compared with the traditional pixel-based method, the object-based method not only includes spectral information, but also gives the object shape and grain features, which can significantly reduce the “salt and pepper effect” (Vieira et al., 2012; Ceccarelli et al., 2013; Jia et al., 2013). In this study, the multi-scale segmentation provided by eCognition is used for image segmentation, and the ESP2 (Estimation of Scale Parameter 2) plug-in is used to determine the optimal segmentation scale to ensure the accuracy of classification and reduce the workload in the later stage (Drǎguţ et al., 2010). We select the threshold classification method to classify the land cover of DL, and then we use visual interpretation to correct the misclassified objects. Referring to the previous research, we divided the land cover of DL into six types: forest (poplar accounted for more than 95%), phragmites, grass, mudflat, shallow water and deep water (Yanxia et al., 2014; Wu et al., 2017). Finally, we use field survey data (30 points, located by GPS) during early December 2019 and Google Earth high-resolution historical images as the true values to establish the confusion matrix (Tian et al., 2019). The results show that the accuracy of classification is 86–90%, and the kappa value is between 0.83–0.86. It can be used for follow-up study (Tang et al., 2016; Hu et al., 2020).

Habitat suitability index (HSI) model is a commonly used model to study the change of habitat suitability of species (Li et al., 2019; Powers and Jetz, 2019). In this study, HSI model was used to quantify the habitat suitability changes of four major migratory bird guilds in DL, including Anatidae, Charadriidae, Ardeidae and Laridae. The first step of establishing HSI model is to select the habitat variables closely related to the survival of the evaluated species. According to previous studies, we selected four groups of habitat variables that affect the foraging and roost of four migratory bird guilds, including land cover type, altitude, slope and human disturbance (Wu et al., 2019b; Tang et al., 2016).

In this study, the preference and selection of migratory birds for a single habitat variable are determined based on people’s understanding of migratory birds’ habits (expert score) (Tian et al., 2008a; Seoane et al., 2006). For example, land cover, which directly determines the habitat type and food availability they can provide (Angradi et al., 2001). The suitability of different types of land cover for four migratory bird guilds are shown in Table 1. Similarly, we have developed a human disturbance, altitude and slope suitability grading system for migratory birds, as shown in Table 2, Human disturbance includes distance to resident district and distance to road, distance to residential area includes distance to urban resident districts and distance to township settlements, distance to road includes distance to first grade road (railway, highway, National Road, provincial road) and distance to secondary road (country road). The routes of these roads and the outline of residential areas were drawn in ArcGIS 10.7 according to Baidu map and Google map. The final score of distance to residential area is the minimum score of its two sub items, so as the distance to road. Elevation and slope data are from the digital elevation model (DEM). It is downloaded from the website: http://www.gscloud.cn/, and the resolution is also 30 m.

The score of single variable and combined score of suitability can be divided into four grades: poorly (0–25), generally (25–50), moderately (50–75) and highly (75–100). All data (including land cover type) are standardized as suitability index for HSI calculation (Zhang et al., 2017). The suitability of land cover type is standardized as the median of suitability classification interval (discontinuous value), and other indexes are standardized by piecewise linear model (continuous value) (Tian et al., 2008b).

The second step of establishing HSI model is to determine the weight of each habitat variable. In this study, the analytic hierarchy process (AHP) coupled with entropy weight method is used to determine the weight to enhance the robustness of the model (Aghajani et al., 2011; Dai et al., 2013). AHP is based on expert experience, which has certain subjectivity, but it can reduce the inconsistency of relative importance judgment through tests. The entropy weight rule provides an objective weighting method by calculating the degree of dispersion between indicators. Therefore, the final weight of each index is set as the average value of the weight determined by the two methods, as shown in Supplementary Table S2. Finally, the habitat suitability index of migratory birds in DL can be calculated by formula (1) H S I j = ∑ i = 1 n w i j f i j (1)

Among them, HSIj is the HSI of the jth bird (j = 1–4), n is the number of habitat variables, wij and fij (i = 1-n) represent the weight and suitability index of the ith habitat variables respectively.

The habitat suitability index calculated by the formula was reclassified to obtain the habitat suitability map of each migratory bird guild over the years. Then, the landscape index software Fragstats (version 4.2.1) was used to calculate the three indexes of classification scale (each grade), including total area (TA), largest patch index (LPI), and mean patch size (MPS).

The calculated landscape index was divided into two groups: before and after PER. The significance of the difference was analyzed by one way analysis of variance (ANOVA). Pearson correlation analysis was used to analyze the relationship between each index and the area of poplar planting. All the above analyses were analyzed by SPSS (version 25).

The results of land cover classification in 6 years before and after PER in DL are shown in Figure 2. It can be seen that poplar in DL is mainly distributed in W-DL and S-DL, while few in E-DL and HL. The areas of various land covers are shown in Supplementary Figure S1. It showed that Phragmites, grass and deep water are the main land cover types of DL. The area of each of them is more than 500 km², and the proportion of forest, mudflat and shallow water is small, with an area of less than 300 km². The poplar area from 2014 to 2016 is 235.76, 248.39 and 245.95 km² respectively, which is close to the official data of 260 km² (2016), showing the accuracy of our land cover interpretation work.

As can be seen from Supplementary Figure S1, the area of poplar remained stable before 2016, with a significant decrease (46.12 km²) in 2017, which shows the significant effect of PER in 2017. According to the report, three million poplar trees were actually cut in 2017 (Xu, 2018), with an area of about 50 km² based on 450–750 trees per hectare, which is consistent with our analysis. Then the poplar area remained stable in 2018, and decreased significantly in 2019 (the second stage of PER). In order to explore the impacts of PER more clearly on the land cover of DL, we extracted the PER area, as shown in Supplementary Figure S2. By calculating the area, we found that 3 years after the poplar forest felled, most of the area was converted into Phragmites (79.75%) and a small part transformed into grassland (15.16%), but very few parts converted into other types (5.09%).

In this paper, we selected the first 3 years and the last 3 years of PER as two groups of data for comparative analysis, and established the habitat suitability model to explore the impact of PER on four main migratory bird guilds in DL. We selected four groups of habitat variables, including land cover type, human disturbance (road distance and residential distance), altitude and slope, which are the four groups of habitat variables most used to evaluate bird habitat (Li et al., 2019; Tian et al., 2019; Wu et al., 2019b). Some studies also selected the normalized difference vegetation index (NDVI) as an indicator of bird food resources (Tian et al., 2008a); Tang et al., 2016; Ivajnšič et al., 2017). In this study, due to the band missing in some images and considering the strong correlation between NDVI and land cover types (Yanxia et al., 2014; Hu et al., 2020), land cover types can also be used as an indirect reflection to food resources (Gan et al., 2009). Therefore, we did not choose NDVI.

We chose the HSI model because it is simple to run and does not need the presence data and absence data of the research species (Hua et al., 2012; Frantz et al., 2018; Lebeau et al., 2019). Compared with maximum entropy (MaxEnt), general additive model (GAM), general liner model (GLM), boost regression tree (BRT) and other models that need these data (Zurell et al., 2012; Brambilla et al., 2020), HSI has no training process, so it needs to determine the weight in advance, and this weight directly affects the accuracy of the model, so we carefully determine the weight, combined with the analytic hierarchy process and entropy weight method. The importance of each habitat variable in AHP is determined according to the previous large number of studies in the same and adjacent areas (Tian et al., 2008a; Tang et al., 2016; Wu et al., 2019b), while the entropy weight is determined according to the characteristics of the habitat variables of DL. Hence, our model should be robust.

For the four bird guilds, the response of each grade habitat in E-DL and HL is not significant, because the poplar area in these two parts is very small, it has a certain relationship with the limited water tolerance of eastern cottonwoods (Cao and Conner, 1999; Glenz et al., 2006). The elevation of E-DL and HL is lower, and the annual flooding days are longer. The elevation of W-DL and S-DL is higher, so they are more suitable for poplar growth. Before PER, the area of poplar in HL accounted for only 3.39%, and that in E-DL was even lower, while that in W-DL and S-DL was 27.87 and 23.43% respectively. The W-DL and the S-DL together account for only 33.91% of the total area of DL, but the area of poplar trees distributed in them accounts for 93.35% of the total area of DL, and the area of poplar retreated from them accounts for 52.44 and 39.95% of the total retreated area of DL, respectively. Therefore, the habitats of the four bird species in E-DL and HL are mainly affected by the common interannual fluctuation caused by other factors, so the changes are not obvious, while the bird habitats in W-DL and S-DL are mainly affected by the land cover type conversion caused by PER, so there are great changes.

Although the land cover types of W-DL and S-DL have changed greatly after PER, only the habitat suitability of Anatidae and Ardeidae has been improved in the end. This is due to the different preferences of these birds for land cover types (Tian et al., 2008a). After poplar trees are cut down, the roots are easy to rot after being soaked in water, and they can return to the previous wetland state as soon as one or 2 years (Liu et al., 2020). The result of the image interpretation shows that most of the poplar felled areas are converted into Phragmites areas, and a small part into grassland, which can improve the habitat of Anatidae and Ardeidae, but for Charadriidae and Laridae, they prefer mudflat and shallow water habitats to Phragmites, and the mudflat and shallow water areas are more strongly affected by hydrological factors (Hu et al., 2015).

In this study, we focus on the effect of PER on the habitat suitability of four migratory bird guilds in DL. However, it is undeniable that our study is also affected by other factors, the most important of which is hydrology (Wang et al., 2019). As mentioned above, the habitat suitability of Charadriidae and Laridae has declined slightly. Through analysis, we find that this is caused by different water levels, and there are two main reasons: firstly, the actual hydrological situation of DL is different every year. Secondly, due to the limitation of the availability of the image, there is a certain interval in the time of our selection. According to the interpretation, the area of deep water area increased significantly, the shallow water area changed little, but the mudflat decreased significantly. This is obviously caused by different water levels corresponding to images in different years (Wu et al., 2020). For this cause, the habitat suitability of Charadriidae and Laridae decreased. The reason for the significant decrease of LPI of the high grade of Laridae habitat in the E-DL is that the habitat is separated by the deep water area. As for the increase of MPS of the high grade in HL, we speculate that the rising of water level may make the water level of HL more suitable for Laridae, because HL is located in the upper reaches of E-DL with higher elevation (Yao et al., 2019).