Response of Anatidae Abundance to Environmental Factors in the Middle and Lower Yangtze River Floodplain, China

Understanding and predicting animal distribution is one of the most elementary objectives in ecology and conservation biology. Various environmental factors, such as habitat area, habitat quality, and climatic factors, play important roles in shaping animal distribution. However, the mechanism underlying animal distribution remains unclear. Using generalized additive mixed models, we analyzed the effects of environmental factors and years on the population of five Anatidae species: Tundra swan, swan goose, bean goose, greater and lesser white-fronted goose, across their wintering grounds along the Middle and Lower Yangtze River floodplain (MLYRF) during 2001–2016. We found that: (1) All populations decreased except for that of the bean goose. (2) The patch area was not included in any of the best models. (3) NDVI was the most important factor in determining the abundance of grazing geese. (4) Climatic factors had no significant effect on the species in question. Our results suggest that, when compared to habitat area, habitat quality is better in predicting Anatidae distribution on the basin scale. Thus, to better conserve wintering Anatidae, we should keep a sufficiently large area at the single lake, as well as high quality habitat over the whole basin. This might be achieved by developing a more strategic water plan for the MLYRF.


Introduction
Understanding and predicting animal distribution is one of the most elementary objectives in ecology and conservation biology. The animal distribution could be affected by various environmental factors, such as habitat area, habitat quality and climatic factors [1]. Despite endeavors to unravel such a fundamental problem during past decades, it still remains not fully understood. Changes in species abundance and distribution might be driven rapidly by local or global environmental changes [2,3], so addressing which is the main driver of species changes at multiple spatial scales is critical for taking conservation measures.
Eastern China is one of the most important wintering regions for migratory waterbirds, supporting over two million individuals, including more than one million Anatidae [4,5]. About 80% of these Anatidae winter within the wetlands in the Middle and Lower Yangtze River floodplain (MLYRF) [5]. However, the Anatidae population size has been declining continuously in this area [5,6]. Understanding how environmental factors affect animal population changes at large tempo-spatial scales, as shown for Anatidae in the MLYRF, is crucial for their conservation [7]. Some case studies have found that habitat area and quality play an important role in determining waterbird distribution in the MLYRF. Zhang et al. demonstrated that a larger patch area concentrated the larger number of grazing Anatidae at Shengjin Lake [8]. Guan et al. found a positive effect of habitat quality (including winter season vegetation index, habitat area and topographic wetness index) on geese density at Dongting Lake [9]. However, all these studies were conducted for one individual lake, and few studies were conducted at larger tempo-spatial scales within the MLYRF [10]. Such a comprehensive study is urgently needed.
The individual-area relationship, which describes the relationship between population density and patch area, has been successfully used to explain the distribution of Anatidae populations in wetland [8]. In MLYRF, grazing Anatidae mainly feed on recessional meadows after the drawdown of the water level, and tuber-feeding birds mainly forage on tubers of submerged macrophytes during wintering season [11,12]. Thus, we expect positive correlations between land/water area and Anatidae populations. As suggested by optimal foraging theory [13][14][15], food quantity and habitat heterogeneity could determine Anatidae distribution through forage intake limitations. As difficulties in handling long leaves and in locating bites with increasing biomass, the intake rate of Anatidae follows a bell-shaped curve with increasing food quantity, which maximizes at the intermediate biomass [15,16]. Habitat heterogeneity, such as the spatial variance of vegetation biomass or structure, generally increases the time needed to search for and handle food, thus, negatively affecting individual intake rates [17,18]. In addition, according to the allometry scaling law, the effects of food quantity and habitat heterogeneity would also differ among species with different body sizes [15,16,19].
Besides habitat, the distribution of birds can be affected by climatic factors, such as temperature and precipitation [20]. For example, a warmer climate is preferred by wintering birds, as more energy is consumed to compensate for colder ambient temperatures [21]. Precipitation is positively correlated with primary productivity within wetlands, which favors herbivores [22]. A high water level disfavors birds, such as tuber-feeders because a higher water level denies them access to the food hidden under the water [23].
In recent years, more and more ornithologists prefer to use contemporary statistical methods, such as generalized additive model. Zenzal et al. used generalized additive mixed model (GAMM) to investigate hummingbird migration across years [24]. La Sorte et al. used GAMM to investigate how land-use change affect avian species [25]. Wen et al. used GAMM to analyze the effects of climate and hydrology on waterbird population dynamics [2].
In this research, we aimed to investigate relationships between the environmental variables and population sizes of five Anatidae species in four key wintering areas in the MLYRF during 2001-2016: Tundra swan Cygnus columbianus, swan goose Anser cygnoides, bean goose A. fabalis, greater white-fronted goose A. albifrons and lesser white-fronted goose A. erythropus. Tundra swan and swan goose mainly feed on tubers of Vallisneria spiralis [12], whereas, the latter three species graze on recessional grasslands in wintering areas.

Study Area
This study was conducted in 43 wetlands in four key wintering areas (Poyang Lake, Dongting Lake, Shengjin Lake and Anqing Lakes) in the MLYRG (see Figure 1 and supplementary materials Table S1). The total area of the 43 lakes is ca. 2800 km 2 . This region has a subtropical monsoon climate that produces a large variation in precipitation among seasons. Wetlands here are characterized by annual recharge from summer rainfall, followed by water level recession in winter, which provides habitats for hundreds of thousands of wintering waterbirds.  Table S1 according to numbers.
Poyang Lakes, which are located in Jiangxi Province 28°24′-29°46′ N, 115°49′-116°46′ E), maintain a direct connection with the Yangtze River, which has a naturally fluctuating water level. The submerged area in the wet season (July to September, ca. 3600 km 2 ) can be more than two times that in the dry season (November to February, ca. 1500 km 2 ) [26]. This variation caused by water level fluctuation produces many types of habitats, such as ephemeral grasslands and mudflats.
Dongting Lakes lie in Hunan Province (28°44′-29°35′ N, 111°53′-113°05′ E), with an annual seasonal water level fluctuation up to 18 m (the variation between maximal and minimal water level). Similar to the Poyang Lakes, the water level fluctuation has not been disturbed in Dongting Lakes, so various habitat types are formed, such as marshes and mudflats. However, the construction of Three Gorges Dam may have prolonged the drought period in Dongting Lakes, since its operation in 2003 [27].
Shengjin Lake and Anqing Lakes are in Anhui Province (29°50′-30°58′ N, 116°07′-117°44′ E). These lakes are much smaller (the biggest is Huangda Lake with an area of ca. 270 km 2 ) and lost their direct connection to the Yangtze River during the 1960s. Due to the different types of management at different lakes, various types of habitats exist as different lakes, such as recessional meadows in Shengjin Lake, and mudflats in Baidang Lake and Caizi Lake.

Census Data
Census data for the five studied Anatidae species were obtained from systematic surveys of these four areas. Wintering waterbirds were counted systematically from 2001 to 2016 (supplementary materials Table S1). In all the surveys, the "look-see" [28] method was used to identify and count all present waterbirds. The observers mostly reached the counting spots by car or  Table S1 according to numbers.
Poyang Lakes, which are located in Jiangxi Province 28 • 24 -29 • 46 N, 115 • 49 -116 • 46 E), maintain a direct connection with the Yangtze River, which has a naturally fluctuating water level. The submerged area in the wet season (July to September, ca. 3600 km 2 ) can be more than two times that in the dry season (November to February, ca. 1500 km 2 ) [26]. This variation caused by water level fluctuation produces many types of habitats, such as ephemeral grasslands and mudflats.
Dongting Lakes lie in Hunan Province (28 • 44 -29 • 35 N, 111 • 53 -113 • 05 E), with an annual seasonal water level fluctuation up to 18 m (the variation between maximal and minimal water level). Similar to the Poyang Lakes, the water level fluctuation has not been disturbed in Dongting Lakes, so various habitat types are formed, such as marshes and mudflats. However, the construction of Three Gorges Dam may have prolonged the drought period in Dongting Lakes, since its operation in 2003 [27].
Shengjin Lake and Anqing Lakes are in Anhui Province (29 • 50 -30 • 58 N, 116 • 07 -117 • 44 E). These lakes are much smaller (the biggest is Huangda Lake with an area of ca. 270 km 2 ) and lost their direct connection to the Yangtze River during the 1960s. Due to the different types of management at different lakes, various types of habitats exist as different lakes, such as recessional meadows in Shengjin Lake, and mudflats in Baidang Lake and Caizi Lake.

Census Data
Census data for the five studied Anatidae species were obtained from systematic surveys of these four areas. Wintering waterbirds were counted systematically from 2001 to 2016 (supplementary  materials Table S1). In all the surveys, the "look-see" [28] method was used to identify and count all Sustainability 2019, 11, 6814 4 of 10 present waterbirds. The observers mostly reached the counting spots by car or on foot. Boats were sometimes used. Large Anatidae (like the swans and geese considered in this study) usually form large flocks during the non-breeding season, allowing their easy identification and counting [29]. Surveys were conducted by staffs from the nature reserve and by the authors using the same survey methods. Multiple teams were made up, and there were at least one experienced observer to guide the survey. It took 5-20 days to complete surveys. See Barter et al. [30] for more details.

Environmental Data
We used Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) surface reflectance (SR) products to extract the habitat features that might potentially affect the waterbird abundance in these areas. The current SR products are provided by the Landsat ecosystem disturbance adaptive processing system (LEDAPS) with satisfying quality and consistency [31].
We selected images that were acquired around our survey dates with less than 10% cloud cover (supplementary materials Table S2). Before processing, we employed a gap-filling method based on local linear histogram [32] matching to remove duplicated data and correct the loss in the ETM+ images, due to the failure of the scan line corrector in 2003. Afterwards, we geometrically corrected each image with an accuracy of fewer than 0.5 pixels.
Before image processing, we delineated the boundaries for all concerned wetlands according to their natural and artificial landscapes ( Figure 1). We then applied supported vector machines (SVMs) to discriminate between water and land within each boundary. SVMs are supervised non-parametric statistical learning techniques that can produce higher accuracy classifications, even without larger training data sets [33]. For each image, we visually selected training sets based on the contrasting differences between pixels containing land or water. All pixels within the study sites were classified as land or water. Patch areas of land or water were then calculated by counting the related pixels. The normalized difference vegetation index (NDVI) derived from the SR products was, thus, extracted from the land as a surrogate of food abundance. We calculated the coefficient of variation (CV) of the NDVI as an index of habitat heterogeneity. Images were processed in ENVI 5.0 (Harris Geospatial Solutions, Broomfield, CO, USA) and ArcGIS 10.2 (Environmental Systems Research Institute, Redlands, CA, USA). We obtained averaged temperature and accumulated precipitation in January for 2001-2016 from Chinese Meteorological Administration.

Statistical Analysis
We used generalized additive mixed models (GAMMs) to model the effects of the patch area, food quantity, habitat heterogeneity, and climatic factors on swans and geese abundance with the survey site as the random factor. The patch area used for grazers (bean goose and greater and lesser white-fronted goose) was the land area, whereas, that used for tuber-feeders (tundra swan and swan goose) was the water area. "Year" was also included to account for the possible temporal correlation of surveys among the years. Usually, the count data contains excess zero records. Poisson distribution is often used to model such count data, but provides poorer performance compared with zero-inflated Poisson or negative binomial models [34]. Thus, based on the diagnosis by Wen et al. [2], we determined the probability distribution function of our count data to be negative binomial type I (NBI) or II (NBII). The general GAMM formulation is as follows: where g(·) is a link function (e.g., NBI or NBII). E(Y) is the expected value of response variable Y. X i is the matrix of explanatory variables with f i (·) as the cubic spine (cs) smoothing functions, and α is a constant. We then followed a stepwise model selection procedure (for the limitations of this approach see Whittingham et al. 2006 [35]), initiated from a full model, including all predictor variables. Since some research suggested that variables with a p-values slight larger than 0.05 also included some useful information when modeling, hence, here variables with p-values higher than 0.1, estimated by a Chi-Squared test, were, thus, omitted. Afterwards, the model was compared with the full model based on generalized Akaike information criterion (GAIC). If the new model has a lower GAIC than the full model, the variable is neglected in subsequent procedures [2]. The goodness of fit of the selected models was evaluated with the pseudo-coefficient of determination (R 2 ), which was used as a surrogate for R 2 . All statistical analyses were conducted in R 3.2.4 [36] with the GAMLSS package [37].

Population Trends
"Year" was included in all the best models except for bean goose. During 2001-2016 lesser white-fronted goose populations increased slightly (p < 0.05, Table 1) within the study area, while swan goose and greater white-fronted goose decreased (p < 0.05, see Table 1). Tundra swan tended to decrease as well, though the year effect did not reach statistical significance (p < 0.1, see Table 1). For swan goose and greater white-fronted goose, the models showed the same pattern: The population size stopped decreasing around 2005, followed by a rise in 2006-2007, and then a sharper decline after 2007-2008 (Figure 2b,c). Tundra swan showed a similar trend, but compared with swan goose and greater white-fronted goose, tundra swan declined later, starting in 2010 (Figure 2a). The rate of decline was faster for swan goose ( Table 1). The population of lesser white-fronted goose increased until 2010, then subsequently decreased (Figure 2d).
For swan goose and greater white-fronted goose, the models showed the same pattern: The population size stopped decreasing around 2005, followed by a rise in 2006-2007, and then a sharper decline after 2007-2008 (Figure 2b,c). Tundra swan showed a similar trend, but compared with swan goose and greater white-fronted goose, tundra swan declined later, starting in 2010 (Figure 2a). The rate of decline was faster for swan goose ( Table 1). The population of lesser white-fronted goose increased until 2010, then subsequently decreased (Figure 2d).

Impacts of Habitat-Related Variables
For grazing Anatidae (bean goose and greater and lesser white-fronted goose), NDVI was the most significant factor positively related to population abundance ( Table 1). The three grazing species responded differently to the increasing NDVI: The population density of the larger species (bean goose and greater white-fronted goose) stopped increasing when NDVI reached 0.4 (Figure 3a,b), whereas, the smaller species, the lesser white-fronted goose, stopped at a lower NDVI value (0.28, Figure 3c).
The NDVI CV, which represented habitat heterogeneity, was negatively correlated to bean goose and lesser white-fronted goose (Figure 3e). The population of lesser white-fronted goose decreased until the NDVI CV was three (Figure 3e). The patch area (land/water) was not selected in the best models (Table 1). Besides, using a linear mixed model with the lake as a random factor, we found no significant relationship between the log-transformed land/water area and the year. (p = 0.23 for log(Land)~year and p = 0.77 for log (Water)~year).

Impacts of Habitat-Related Variables
For grazing Anatidae (bean goose and greater and lesser white-fronted goose), NDVI was the most significant factor positively related to population abundance ( Table 1). The three grazing species responded differently to the increasing NDVI: The population density of the larger species (bean goose and greater white-fronted goose) stopped increasing when NDVI reached 0.4 (Figure 3a,b), whereas, the smaller species, the lesser white-fronted goose, stopped at a lower NDVI value (0.28, Figure 3c).
The NDVI CV, which represented habitat heterogeneity, was negatively correlated to bean goose and lesser white-fronted goose (Figure 3e). The population of lesser white-fronted goose decreased until the NDVI CV was three (Figure 3e). The patch area (land/water) was not selected in the best models (Table 1). Besides, using a linear mixed model with the lake as a random factor, we found no significant relationship between the log-transformed land/water area and the year. (p = 0.23 for log(Land)~year and p = 0.77 for log (Water)~year).

Effects of Climatic Variables
Precipitation was positively related to the numbers of tundra swan and bean goose, whereas, temperature was negatively related to numbers of swan goose. Though included in the best models, the effects of precipitation or temperature were not significant (p > 0.1; Table 1).

Effects of Climatic Variables
Precipitation was positively related to the numbers of tundra swan and bean goose, whereas, temperature was negatively related to numbers of swan goose. Though included in the best models, the effects of precipitation or temperature were not significant (p > 0.1; Table 1).

Discussion
We analyzed the effects of environmental factors and year on five Anatidae species during 2001-2016. The results suggest that: (1) During the study period, most populations have been decreasing, except for bean goose. Swan goose declined the fastest over the study period, whereas, the number of lesser white-fronted goose declined quickly since around 2010. (2) The patch area was not included in any of the best models. (3) NDVI was the most important factor in determining the abundance of grazing geese. (4) Climatic factors had no significant effect on any species.
Populations of tundra swan, swan goose and greater white-fronted goose showed a clearly decreasing trend, and lesser-fronted goose population stopped increasing until 2010 ( Figure 2). In eastern China, more than 80% of the population of these five species spend their non-breeding period in the four key areas studied here [6]. The generally decreasing trends suggest possible degradation within these four areas. Zhang et al. found that the population in the National Natural Reserves had a lower decreasing rate [10]. However, only 6% of wetlands are protected by National Natural Reserves in the MLYRF [38,39], which seems insufficient to improve the generally decreasing population trends, especially considering their current status. Notably, the population of tuber feeders (i.e., tundra swan and swan goose) decreased faster than that of grazers. Tundra swan and swan goose mainly forage on tubers of submerged macrophytes, which have been dying out in the wetlands of the MLYRF, due to intensive human activities in recent years [12]. The diminishing food resources likely explain the decline in tuber feeders. Grazing geese are affected by habitat changes, which might be a consequence of the hydrological changes caused by dam operation [40,41]. Our results suggest that tuber feeders are more sensitive to habitat degradation. The bean goose diet is

Discussion
We analyzed the effects of environmental factors and year on five Anatidae species during 2001-2016. The results suggest that: (1) During the study period, most populations have been decreasing, except for bean goose. Swan goose declined the fastest over the study period, whereas, the number of lesser white-fronted goose declined quickly since around 2010. (2) The patch area was not included in any of the best models. (3) NDVI was the most important factor in determining the abundance of grazing geese. (4) Climatic factors had no significant effect on any species.
Populations of tundra swan, swan goose and greater white-fronted goose showed a clearly decreasing trend, and lesser-fronted goose population stopped increasing until 2010 ( Figure 2). In eastern China, more than 80% of the population of these five species spend their non-breeding period in the four key areas studied here [6]. The generally decreasing trends suggest possible degradation within these four areas. Zhang et al. found that the population in the National Natural Reserves had a lower decreasing rate [10]. However, only 6% of wetlands are protected by National Natural Reserves in the MLYRF [38,39], which seems insufficient to improve the generally decreasing population trends, especially considering their current status. Notably, the population of tuber feeders (i.e., tundra swan and swan goose) decreased faster than that of grazers. Tundra swan and swan goose mainly forage on tubers of submerged macrophytes, which have been dying out in the wetlands of the MLYRF, due to intensive human activities in recent years [12]. The diminishing food resources likely explain the decline in tuber feeders. Grazing geese are affected by habitat changes, which might be a consequence of the hydrological changes caused by dam operation [40,41]. Our results suggest that tuber feeders are more sensitive to habitat degradation. The bean goose diet is relatively more diverse [42], and less sensitive to habitat changes. Thus, the population of bean goose showed no clear temporal trends.
The patch areas were not featured in the best models for any of the five species. This result is inconsistent with those reported by Connor et al. [43] and Zhang et al. [8], who found a positive effect for the patch area. We found no significant relationship between the log-transformed land/water area and the year, which might explain the lesser contribution from patch areas. We found NDVI to be the most important factor influencing grazing geese, which is similar to the findings of Zhang et al. [10]. The population size firstly increased, and then decreased with an increase in NDVI. Lesser white-fronted goose stopped increasing when the NDVI reached 0.28, whereas, bean goose and greater white-fronted goose stopped increasing at 0.4. Compared with the bigger bean goose and greater white-fronted goose, the NDVI value at which the species stopped increasing was lower for the lesser white-fronted goose, implying a possible effect of the allometry law [44,45]. The NDVI CV was generally negatively correlated with the bean goose population, but this correlation was not significant; this result is also in line with that of Zhang et al. [8]. As expected, the response of the lesser white-fronted goose population to the increasing NDVI CV was an inverted bell curve. Unlike the widely distributed sedge Carex meadow, the grass Alopecurus aequalis and the spikerush Eleocharis migoana, which are the main food sources of the wintering lesser white-fronted goose, are seldom found beyond the heterogeneous recessional grassland in East Dongting Lake [41]. Combined with heterogeneity negatively affecting the foraging efficiency of geese by increasing searching and handling time [17,18], we conclude that lesser white-fronted goose would favor a habitat with mild heterogeneity. In general, such results suggest that, mediated by the allometry law, geese might prefer habitats of higher quality to larger areas at basin scales [46,47].
Our study showed that wintering Anatidae in the MLYRF have generally declined. We also showed the important effects of NDVI (i.e., habitat quality) on grazing Anatidae at the basin scale. At this scale, patch area or climatic factors seem to play insignificant roles in determining Anatidae populations. According to Zhang et al., increasing the patch area by managing the water level of a lake would effectively protect wintering geese [8]. However, our results suggest that such management actions might not benefit wintering geese at the basin scale. To better conserve the wintering Anatidae, we need to integrate priorities at both the lake and basin scales, which means that enough patch areas should be retained at a single lake, as well as high habitat quality over the whole basin. As suggested by Zhao et al., reductions in grazing goose numbers at East Dongting Lake were correlated with declines in the availability of suitable sedge swards, caused by earlier water table recession, due to the commissioning of the Three Gorges Dam since mid-2003 [40]. Thus, habitat quality could be improved by hydrological regimes, especially the timing of water draw down [9]. In conclusion, we might benefit wintering Anatidae by developing a more strategic water plan for the MLYRF.