Declines in Common and Migratory Breeding Landbird Species in South Korea Over the Past Two Decades

Population declines in terrestrial bird species have been reported across temperate regions in the world and are attributed to habitat loss, climate change, or other direct mortality sources. North American and European studies indicate that long-distance migrants, common species, and species associated with grasslands and agricultural lands are declining at the greatest rates. However, data from East Asia on avian population trends and associated drivers are extremely sparse. We modeled changes in occupancy of 52 common breeding landbird species in South Korea between 1997–2005 and 2013–2019. Thirty-eight percent of the species showed evidence of declines, and seven of these were declining severely (46–95%). Occupancy of Black-capped Kingfisher (Halcyon pileata) populations have dropped the most precipitously over the study period. Among declining species, long-distance migrants (9/20) and common species (14/20) showed more rapid declines than other groups. Declines of five species were associated with climate change, and two species appeared to be affected by land-cover change. However, causes of change in occupancy of other species (46/52) remains cryptic. Based on our results, we suggest an immediate re-evaluation of species’ conservation status and legal protection levels for seven severely declining species in South Korea, and a dedicated survey design and analysis effort for the continued monitoring landbird populations. Because many species exhibiting declines migrate from beyond national boundaries, international collaborations will be required to better quantify population trends across the full annual cycle, and to understand mechanisms for these declines.


INTRODUCTION
In recent decades, a rapid decline in avian biodiversity has been reported around the world. Studies from Europe (Inger et al., 2014;Gregory et al., 2019), Canada, and the United States (Rosenberg et al., 2019) show rapid and substantial declines in common and widespread species across wide geographic ranges. Habitat loss, degradation and climate change are expected to be the most substantial anthropogenic drivers of avian biodiversity loss and population declines (Thomas et al., 2004;Sekercioglu et al., 2008;Pearce-Higgins et al., 2015).
Increasing pressure on land-use from human population and economic growth constrains biodiversity conservation efforts, especially where human population density is high and economic growth has increased (Lambin and Meyfroidt, 2011;Di Marco et al., 2018). The loss of avian biodiversity in common and widespread species from land-cover and climate change can degrade ecosystem services and ecological functions provided by avifauna (Şekercioǧlu et al., 2004;Cardinale et al., 2012).
Agricultural land covers at least a third of the earth's icefree land area in the Anthropocene (Ramankutty et al., 2018). While the conversion of native vegetation for agricultural use (both cropland and pasture) poses a dominant threat to avian biodiversity (Gaston et al., 2003), species that are adapted to open vegetation and agricultural landscape are also declining in more populated regions of the world (Stanton et al., 2018;Reif and Hanzelka, 2020). Declines in bird populations that breed in agricultural areas are known to be related to intensive management practices such as the increase in pesticide use or loss of functional diversity from the homogenization of vegetation (Hallmann et al., 2014;Šálek et al., 2018). However, abandonment and conversion of farmlands, especially the loss of low-intensity, small scale traditional farming, has resulted in local biodiversity loss in different studies (MacDonald et al., 2000;Queiroz et al., 2014;Katayama et al., 2015).
Migratory and insectivorous species are declining in temperate regions (Sanderson et al., 2006;Bowler et al., 2019;Rushing et al., 2020). Climate and land-use change often pose greater threats to migratory species than resident species due to the potential effects of these factors on the availability of arthropod prey and thermoregulation during periods of high energy demands (Both et al., 2010;McKechnie and Wolf, 2010;Pearce-Higgins et al., 2015). Declines in the availability of invertebrate foods have been hypothesized to be related to changing thermal conditions (Lister and Garcia, 2018) and intensified land use (Sánchez-Bayo and Wyckhuys, 2019), potentially leading to the functional collapse of food webs (Hallmann et al., 2017;Seibold et al., 2019). However, these patterns have been described mainly for temperate ecosystems in North America and Europe, where there are long-term records of avian biodiversity data (Van Strien et al., 2001;Sauer et al., 2013). Many parts of the world still need systematic collection and analysis of long-term data to quantify and understand these changes (Proença et al., 2017). In East Asia, there is a long history of anthropogenic modification of landscapes through rice-paddy dominated agriculture, forest exploitation and recent recovery, and conversion of agricultural lands to urban structures (Aikens and Lee, 2013;Ramankutty et al., 2018).
Up-to-date systematic assessment of breeding land birds has been rare in this region (but see: Yamaura et al., 2009;Lin and Pursner, 2020) except for regional long-term waterbird monitoring (Mundkur et al., 2017;Amano et al., 2018). While declines in the global population of once-common species of buntings in Asia have recently been documented (Kamp et al., 2015;Tamada et al., 2017), the population status of most species is uncertain (Yong et al., 2015). Given past and current changes in land use (Song et al., 2018) and global climate change effects in this region (Gu et al., 2018), systematic monitoring and analysis are necessary for detecting changes which could inform conservation and management actions.
In South Korea, rice paddies have decreased by more than 33% between 1980 and 2014 (Choi et al., 2016;National Geography Information Institute, 2019). Open agricultural fields are often converted to urban land cover (125% increase between 1990 and 2010) and greenhouse facilities (85% increase between 1991 and 2012; National Geography Information Institute, 2019). Meanwhile, forest cover in Korea has remained relatively stable following a rapid increase between the 1960s and 1990s, due to intensive reforestation and forest protection efforts after the Korean War (1950)(1951)(1952)(1953) (Tak et al., 2007;Bae et al., 2012). If occupancy trends are a function of habitat loss and degradation, we predicted that these changes (Figures 1B-D) in land use and climate should be associated with large-scale changes in the occupancy of breeding landbird populations in South Korea. More specifically, we expected that occupancy of species that use forest and urban areas as breeding habitat should be stable or increasing. However, species that rely on agricultural ecosystems should have declined, especially those associated with rice paddies. Migratory species that have been affected by climate change and land-use change across their migratory cycle may have been declined more severely. Based on previous studies and recent environmental changes in Korea, we suspected that there could be substantial declines in breeding populations of many species, that would be reflected in the estimated occupancy dynamics.
We used the South Korean national bird survey data to model the broad-scale occupancy trends of 52 breeding landbird species between two survey periods (1997-2005 and 2013-2019). We also tested whether land-cover and climate change, as well as species' ecological traits, could explain inter-species variation in trends. To our knowledge, this is the first comprehensive multi-taxa analysis of long-term avian occupancy dynamics for the breeding landbird populations in East Asia.

Korean Bird Atlas Data, Species Criteria, and Spatial Units
The National Natural Environment Survey (NES) in South Korea is one of the first systematic and legally binding nationwide surveys that included terrestrial avifauna (National Institute of Ecology, 2017; . The first survey started in the late 1980s, reporting checklists for species occurrence on a broad provincial scale, without recording the search effort (time and space) or consistency in survey methods. Such a lack of baseline information poses substantial challenges for data analysis. The second survey phase (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) covered all of South Korea. In this phase, the field surveys were conducted in watershedbased survey units with streams, coastlines, and drainage divide lines separating the sampling units in different sizes and shapes (National Institute of Ecology, 2017). For the 1997-2005 period, experts in field ornithology were selected from universities and other institutions; for more recent surveys, members of public who are certified through the Korean Ministry of Environment (K-MOE)'s training program were included in the surveyor pool . Surveyors were asked to visit and record species, breeding status (in categories), environment and their relative abundance based on repeated visits in different seasons to reflect breeding and migration periods. Each surveyor group was assigned to multiple adjacent units and reported a list of birds observed from each of these watershed-based units. However, some survey reports did not provide explicit lists for each unit and instead, reported a combined list for all units they surveyed. We excluded these lists from the analysis. Even with multiple visits, the published data were combined in a single table for each unit, making NES a single-visit, atlas-type survey for the 1997-2005 period. For the third (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) and fourth survey (2014-2019) phases, sampling units were set to rectangular plots with nine grid cells of 994 national reference maps (at a scale of 1:25,000; 2 3 × 2 3 grids in latitude and longitude) across South Korea (National Institute of Ecology, 2017). Survey reports from this period provide date of each observation, but do not provide exact survey duration or time of observation.
However, the number of "visits" were reported for each grid (post-2006) or survey units (pre-2006), and we used this information as a measure of relative effort. Each visit consists of a set of consecutive days or number of visits described by the surveyors. Two observers were assigned for each survey unit, where they conducted surveys on roadsides or well-established hiking trails. We reviewed and compiled bird survey reports from 1997-2005 to 2013-2019 (hereafter survey periods). And extracted bird lists for the smallest spatial units and gleaned information on the number of visits per species list. We selected most recent round of surveys (2013-2019) and earliest period (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) that cover the study area twice. Survey reports are available from the K-MOE Digital Library 1 , in Korean. We supplemented missing survey data in 2014-2019 from 2013 to provide better spatial coverage. For selecting breeding records of earlier (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) survey data, we used the information from the reports to filter out non-breeding records, using breeding status, list of survey dates, and survey environment. The data from the later periods came with dates of observation, so we selected all records from April to August, a period that represents the general breeding season in south Korea for our target species. Next, we removed all waterbird species that are covered in the Asian Waterbird Census (Mundkur et al., 2017) and removed rare and localized species that only occur in limited habitats (i.e., coastal islands, subalpine zones). We also excluded species pairs that can be easily misidentified species (i.e., Bush Warblers; Horornis diphone / Horornis canturians complex) and nocturnal species (owls and nightjars; order Strigiformes and Caprimulgiformes). Since having no false-positive identification is one of the most important assumptions in occupancy models (MacKenzie et al., 2002), we carefully removed species that are prone to misidentification. By excluding rarer species and only including conspicuous species that can be detected and identified easily with sound and sight, we tried to avoid the violation of this key assumption. Lastly, we removed species with <10% occurrence on all sites in both survey periods. We selected all migratory bird species (41 species) and an additional 11 species of selected sedentary passerine species from what remains from the above filtering process, leaving a total of 52 species for the analysis (Supplementary Table 1). We did not include woodpeckers as the primarily road-side survey design was not appropriate for these taxa that mainly use the forest interior and are underdetected in passive counts (Warren et al., 2005;Kumar and Singh, 2010;Saracco et al., 2011). Among the species selected, three species (Oriental Turtle Dove Streptopelia orientalis, Eurasian Magpie Pica pica, and Large-billed Crow Corvus macrorynchos) are game animals in South Korea; however, we did not consider the impact of direct human-induced mortality such as legal or illegal hunting on landbirds in this study, because the bird harvest, persecution, and poaching activities are believed to have been well-controlled in South Korea since the 1990s .
Because the spatial extent and grain size of each survey unit in the early surveys (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) do not match with the later years (2013-2019), we aggregated survey units in later years to the earlier survey period's spatial scale. To achieve this spatial matching, we overlaid the rectangular plots from the third and fourth phase on the second phase survey units, thereby creating a new sampling unit (hereafter "sampling unit"; aggregation of 1-43 third/fourth phase plots to second phase survey units, mean area = 13,840 hectares and standard deviation = 7,086 hectares). Whenever there was a grid unit from a later period that crossed the boundaries of an earlier period's survey unit, we used the preponderance of area included in each grid as a criterion for aggregation. Because the grids from 2013 to 2019 are fairly large, the boundaries did not match perfectly. This spatial error may introduce errors in parameter estimates, but it should not have a systematic bias, as the process is basically converting vector sampling units into raster cells. We selected survey units that had surveys from both periods (1997-2005 and 2013-2019), but excluded survey reports that did not provide the number of visits for sampling units for each survey period. Following these data filters, we obtained 464 sample units for both periods for use in our analysis ( Figure 1A).
We reviewed species accounts (Gore and Won, 1971;Billerman et al., 2020) to classify each species into simple categories based on their breeding season diet (arthropod, vertebrate, and others), general habitat association (i.e., farmland and forests), and migratory behavior (i.e., sedentary, shortdistance migrant, long-distant migrant). Most species were assigned to the arthropod diet group or the vertebrate group diet. We classified crow and magpies into the vertebrate diet group together with raptors, shrikes and kingfishers, and two remaining dove and grouse species into the "others" category (Supplementary Table 1). Given the lack of definite classification in this region, we considered short-distance migrants as those species that are partial migrants or migrating within Korea, Japan, and East China, whereas long-distance migrants were those that winter in tropical climate zones in Southeast Asia, South Asia, and South Africa (Supplementary Table 1).

Environmental Covariates: Land Use and Cover, Topography and Climate Data
We obtained land-use and cover maps produced by the Korean Ministry of Environment for the 1997-2000 period (classified from 5-m resolution imagery) and 2018 (classified from 1-m resolution imagery) from Environmental Geographic Information Services 2 and used these to assess land-use and cover amounts in each period. As two maps were created from satellite imagery with different resolutions, we rasterized the provided vector maps to 30-m resolution for final use. We grouped detailed land-use and cover classes to four general landcover variables. We reclassified all roads, buildings, other built structures, and landscaped artificial green spaces into "urban, " all dry, non-rice open agricultural lands to "dry-field, " all forest types to "forest, " and all rice paddies to "ricefield." We then extracted the area of each variable per survey unit for both survey periods (Figure 1). We used historical monthly climate data for South Korea from CRU-TS 4.03 (Harris et al., 2014), downscaled in WorldClim 2.1 (Fick and Hijmans, 2017). We chose two climate covariates, breeding season (May-June) average daily maximum temperatures and cumulative precipitation for the same period, matching each survey unit's survey years. We then extracted the average values per each survey unit from summarized 1-km resolution raster datasets of two variables. Also, elevation and topographic roughness indices were extracted from a 10-m resolution digital elevation model. All covariates were standardized with mean and variance (average is set to zero) for the model fitting. No pairwise relationships between covariates showed a high level of correlation (r < 0.70). We prepared all environmental covariate data using the following R packages: "raster" (Hijmans, 2020), "exactextractr" (Baston, 2020), "sf " (Pebesma, 2018), and "dplyr" (Wickham et al., 2020).

Modeling Changes in Bird Occurrences
We modeled bird occurrences and dynamics between two survey periods (1997-2005 and 2013-2019) in our study area using the single-visit dynamic occupancy model approach (Peach et al., 2017). This approach uses an effort variable to estimate detection probability instead of repeated observation data in robust multiple-visit surveys (MacKenzie et al., 2003). This allowed us to model each bird species occurrence while accounting for imperfect detection in bird survey data with no repeated visits (Lele et al., 2012). The model uses the combined probability of detection (P i,j * ) for detection probability given the species is present (p i,j ) at each observation per site j and season i, for unit effort (E i,j ) and covariates (x 1 ... x 1+k ) with intercept (β 0 ) and coefficients β 1+k for each x 1+k covariates. We modeled each survey period's detection probability independently as two separate sub-models (p 1j, and p 2j ), each with separate intercepts. For the effort variable, we used the number of visits to each sampling unit.
Dynamic occupancy models explicitly model initial occupancy (ψ 1,j ) and two dynamic parameters, the probability of colonization (γ i,j ) and probability of extinction (ε i,j ) at each site j, and primary periods i, while accounting for imperfect detection by adjusting observational data with estimated detection probability (MacKenzie et al., 2003). Note that since we have two primary periods (1997-2005 and 2013-2019), one parameter will be estimated per site for each probability of colonization and extinction. logit The initial occupancy model for each species includes the land-cover area variables of the first period (cover type_t1), breeding season average daily maximum temperatures (tmax_t1), and precipitation of the first period (precip_t1). The landcover variables included land-cover types of forest, dry-field, ricefield, and urban cover, as described above. Two dynamic parameters (probability of colonization: γ j and the probability of extinction: ε j ) were modeled with the land-cover change variables and climate change variables. We defined landcover as both proportional change in land cover between two periods and the remaining amount of land-cover at the second survey period. To incorporate this definition, we included both land-cover area the second period (land-cover area; cover type_t2), and their % change from the first period (landcover change; cover type_d) per each land-cover type. For climate variables, we used the difference between two periods (precip_d, tmax_d) for dynamic parameter submodels. Because our sampling units are not homogenous in size and shape, we always included the land-cover area variable and landcover change variable of each land-cover type together as a pair for dynamic parameter models (Supplementary Materials S2). By using this approach, we intended to reflect the overall land-cover type change, including both the retained amount of land-cover type after the change and the magnitude of land-cover change between the two periods. The detection probability sub-models included the sampling unit area, elevation, and terrain roughness index (TRI; Wilson et al., 2007) as covariates, accounting for accessibility and spatial extent of each sampling unit.
We built 64 candidate model sets for the probability of colonization, extinction, and initial occupancy, and eight candidate models for each detection probability models with these covariates. We fit models using maximum likelihood and selected the best subset using Akaike's Information Criteria (AIC) and evidence ratios for inference (Burnham and Anderson, 2004). Instead of comparing the full combination of these submodels sets, which would result in more than 16.7 million (64 3 × 8 2 ) possible models, we used a "secondary candidate set strategy" described by Morin et al. (2020). This ad hoc model selection strategy selects the best model sets for multilevel models in two stages of model selection. In the first stage, the candidate sets for each sub-models (ψ 1,j , γ i,j , ε i,j , p 1,j and p 2,j ) are fitted with all candidate target sub-model structures while non-target sub-models (e.g., first fitting and finding first stage candidate sets of initial occupancy while holding all other submodel structure to most general model). We then established a secondary model set using the selected model sets from the first step that had AIC < 5. Secondary model set included all possible combinations of sub-models selected from the first step. We used evidence ratios (ER) for model selection criteria for this secondary set evaluation, accepting models with a likelihood of support (ER) of 0.5 (<2 evidence ratio) than the best model. This approach is computationally efficient (208 models per species for first stage per each species) while giving comparable results to the "all-plausible combinations" approach (Morin et al., 2020).
Occupancy is often used as a surrogate for abundance (MacKenzie and Nichols, 2004) and tends to be positively correlated with abundance (Strayer, 1999;Gaston et al., 2000;Zuckerberg et al., 2009). Thus, we used our estimates on occupancy and occupancy trends as a proxy for population changes in each species (Steenweg et al., 2018). The correlation between occupancy and population size is imperfect and scale dependent. This is partly due to inter-species variation in home range sizes (Steenweg et al., 2018). In the case of the NES, the size of sampling unit varied over time, which likely further complicated the relationship between abundance and occupancy. Nevertheless, spatial matching of two surveys allowed us to infer the magnitude of change in occupancy between two periods across the study area (Zuckerberg et al., 2009). As our main interest was in identifying changes in occurrences that reflect trends in populations across South Korea, we derived mean occupancy in the second survey period (2013-2019) for each species (ψ t = 2 ), and mean trend parameter (rate of change in occupancy; λ = ψ t = 2 /ψ t = 1 ). We used parametric bootstrapping (bootstrap sample size n = 1,000) to provide uncertainty measures for all parameter estimates. We conducted goodnessof-fit tests for final models following MacKenzie and Bailey (2004) using Pearson's Chi-square statistics to check the fit of the selected models. Bootstrap goodness of fit using Pearson's Chi-square confirmed acceptable fit of all models from the final set (p > 0.58; Supplementary Table 2), though for species with very high or low prevalence, the overdispersion parameter (ĉ) was very low (<0.1). We also assessed for agreement within each species' model sets by visually assessing bootstrapped distributions and departures from the maximum likelihood estimates. With the final set of the model, we inferred each parameter from pooled bootstrapped parameters for evaluating each species trend, using confidence intervals (95% for strong evidence and 90% for moderate evidence). We considered a species' trend parameter (λ) to support a statistically significant decline if the upper confidence interval was below 1, and an increase when the lower confidence interval was above 1. For all species' model sets, we examined the effect of each land-cover change and climate change covariate using bootstrapped 95% confidence intervals of beta coefficients for these model terms in the dynamic parameter sub-models. We combined beta coefficients of land-cover area covariates for the 2013-2019 survey period (cover type_t2) and land-cover change covariate (cover type_d) to estimate the overall effect of the land-cover change on dynamic parameters. Single-visit dynamic occupancy models were fitted using the "optim" function in the base package of R, with the Newton-Raphson algorithm ("BFGS"; Nash, 1990;R Core Team, 2020). The code for specifying the maximum likelihood of single-visit occupancy models was modified from Peach et al. (2017Peach et al. ( , 2019.

Trait Group and Species Status Analysis
We estimated the magnitude of change in declining species for each ecological trait group and initial occurrence level (using the median of initial predicted occupancy = 0.63 as a cutoff for common versus uncommon species) to characterize ecological traits associated with declining species. We used linear mixed model analysis of variance to compare the differences between the probability of occupancy of 1997-2005 and 2013-2019 ( ψ = ψ t = 2 − ψ t = 1 ). We used bootstrapped estimates of the ψ to account for uncertainty in each point estimate. The model was specified with "species" as a random effect and each ecological trait (i.e., migratory behavior, diet, habitat association, and their initial abundance) as fixed effects. We conducted post hoc pairwise comparisons for the trait groups with some evidence of differences from Analysis of Variance. Alpha level was adjusted for multiple comparisons using the Tukey method. We estimated marginal means for each trait groups and plotted them to show differences among categories within each trait group. We fitted generalized linear models using package "lme4" (Bates et al., 2015) with the default maximum likelihood estimation option. Pairwise comparisons and estimation of marginal means were conducted using R package "emmeans" (Lenth, 2020).
Lastly, we used IUCN's red list criteria (IUCN Standards and Petitions Subcommittee, 2016) to evaluate the status of declining birds for the national level and compared these with multiple conservation status assessments, including the current IUCN Red List (IUCN, 2020), a current national red list prepared by the Korean Ministry of Environment (Suh et al., 2014). We also compared avian conservation status with legal protection status based on the "Wildlife Protection and Management Act (Act No. 15835)" and "Cultural Heritage Protection Act (Act No. 15827)" (Korea Legislation Research Institute, 2020). The former act mandates the legal protection of endangered species listed by the K-MoE, and the latter protects species that are listed as "natural monuments." We used program R version 4.0.2 (R Core Team, 2020) for all statistical analyses.

RESULTS
We selected 89 final sets of plausible models for all 52 species from the model selection procedure (Supplementary Table 2). Twenty species of the 52 that met our criteria for inclusion (38.4%) showed moderate levels of evidence for declining trends in occupancy, ten species (19.2%) increased during the study period, and for the remainder we did not detect statistically significant occupancy changes (Figure 2). Browneared Bulbuls (Hypsipetes amaurotis), Chinese Sparrowhawk (Accipiter soloensis), Daurian Redstart (Phoenicurus auroreus), Japanese Wagtail (Motacilla grandis), and Northern Hawk Cuckoos (Hierococcyx hyperythrus) exhibited moderate evidence for declines. The remaining 22 species showed no evidence of change or remained stable during the study period (Supplementary Table 4 and Figure 2). Most of these declining species (15 species) had initial occupancies over 50% across the study area in 1997-2005, while the other five species (Yellowrumped Flycatcher Ficedula zanthopygia, Ruddy Kingfisher Halcyon coromanda, Northern Hawk-Cuckoo, Lesser Cuckoo Cuculus poliocephalus, Brown Shrike Lanius cristatus) were less common than others. The magnitude of change varied greatly among these declining common species, from a 3% decline in Brown-eared Bulbuls to a drastic, 95% decline in Black-capped Kingfishers (Halcyon pileata). The probabilities of occupancy of all five less common declining species declined more than 50% between the two survey periods. Increasing species had moderate levels of initial occupancy overall (<80%) and increased from 4% (Bull-headed Shrike Lanius bucephalus, White's Thrush Zoothera aurea) to 81% in the Large-billed Crow. The distributions for each species' parameter estimates are provided in the Supplementary Materials ( Supplementary  Figures 1-4).
From 89 selected models for all species, only 31 models for 12 species retained local extinction covariates in the models, and 44 models for 16 species retained local colonization covariates. For two species, Chinese Sparrowhawk and Eurasian Hoopoe (Upupa epops), landscapes with increased amounts of ricefields and forest, and the remaining amount of each land cover types, respectively, had the lowest rates of local extinction. The other four species' (Gray Starling Spodiopsar cineraceus, Red-rumped Swallow Cecropis daurica, Eurasian Kestrel Falco tinnunculus, and Hazel Grouse Tetrastes bonasia) colonization probabilities were related to changes in breeding season mean daily maximum temperatures and change in breeding season precipitation. Precipitation effects on colonization were speciesdependent, with rainfall having a positive effect on Eurasian Hoopoe colonization, and a negative effect for Hazel Grouse. Aside from these six species, however, most species' change in occupancy was not explained by land-cover and climate covariates (Figure 3).
From our review of the IUCN Red List and the South Korean Red List, we learned that 14 species out of 20 species that were declining our study have no quantitative support for their population trend statements ( Table 1), even though the IUCN red list describes "trend unknown" for only one FIGURE 2 | Distribution of bootstrapped parameter estimates of the probability of occupancy (occurrence) of each species in 1997-2005 (dotted line) and 2013-2019 (solid line) period and rate of change between the two periods. Species are ordered by the magnitude of change from low to high. Each density distribution can be used to understand the uncertainty of estimates (more uncertainty when it is flat and wide, for average occupancies and trend parameters). Note that the rate of change is log-transformed with two base. Thus, one unit change refers to a two-fold change in the probability of occupancy.  Frontiers in Ecology and Evolution | www.frontiersin.org species (Gray Starling). Seven declining species are eligible for threatened categories (Endangered and Critically Endangered) at the national level when we apply IUCN's criteria to our findings. Black-capped Kingfisher, which declined 95%, is eligible for Critically Endangered, and other six species (Brown Shrike, Japanese Wagtail Motacilla grandis, Lesser Cuckoo, Northern Hawk-Cuckoo, Ruddy Kingfisher, and Yellowrumped Flycatcher) are eligible for the Endangered category at the national level based on A2 criteria in IUCN guidelines (population decline with the source of threat unidentified nor ceased). However, the current South Korean Red List published by the K-MoE has classified only the Chinese Sparrowhawk in the Vulnerable category, based on the D1 criteria (abundance less than 1,000 individuals). No species on our list were classified as threatened in the IUCN's global Red List. Among declining species from our study, only two species are legally protected in South Korea as endangered species class II (Chinese Sparrowhawk) and natural monument (Lesser Cuckoo and Chinese Sparrowhawk).

DISCUSSION
While occupancy estimates for 32 species were stable or increasing, about 38% of the species (20/52) we examined were declining, and the majority of those species were common (14/20) or migratory (14/20) species. For species that have a small breeding range in East Asia, such as Yellow-rumped Flycatcher, Japanese Wagtail, and Northern Hawk Cuckoos, our trends could reflect a high proportion of global population trends for these species. The general pattern of decline in South Korea follows similar studies in Europe (Inger et al., 2014;Gregory et al., 2019) and North America (Rosenberg et al., 2019), where common and widespread breeding bird species have declined rapidly in recent decades. In East Asia, studies of Japanese breeding bird data (Yamaura et al., 2009) using unadjusted encounter data have found similar results. The relationship between environmental drivers and dynamic parameters was cryptic for many species; few of the land-use or climate variables we identified were predictive of occupancy trends. This leaves us without clear suggestions as to the primary drivers of breeding bird declines in South Korea, despite apparent declines in 20 species.
The declines of common birds are a serious concern because similar percentage occupancy losses in these species equates to many more individual birds (and even higher losses to avian biomass) than declines in rare species (Gaston et al., 2018). Such species are much more likely to contribute to ecosystem processes and services (Smith and Knapp, 2003;Maas et al., 2016). For example, breeding landbird populations may play a key role as keystone predators of arthropod communities and populations (Terborgh, 2015;Nyffeler et al., 2018). In addition to their roles in trophic interactions, some of these declining common species are also known to provide other ecosystem services and functions such as seed dispersal (Brown-eared Bulbuls;Fukui, 1995;Kim et al., 2015) and scavenging (Eurasian Magpies; Inger et al., 2016). Thus, the observed decline of once-common breeding landbirds in this study suggests associated changes in the bird-related communities, ecosystem services, and functions in the study area.
Migratory species, especially long-distance migrants, have suffered greater declines between the 1997-2005 and 2013-2019 periods. Our results agree with the general patterns of sharper declines in long-distance migrants from studies in North America and Europe (Vickery et al., 2014;Rosenberg et al., 2019), as well as with long-distance migrant declines in Japanese breeding birds (Amano and Yamaura, 2007;Yamaura et al., 2009). Migratory bird species are affected by climate change and habitat losses across separate breeding, migration and wintering locations (Buehler and Piersma, 2008;Reudink et al., 2009;Faaborg et al., 2010), which amplifies the potential exposure to these and other stressors. For example, two declining species from our study, Common and Lesser Cuckoos, are extreme long-distance migrants -wintering in Africa and India (Erritzøe et al., 2012). Previous studies suggest that these cuckoo species are known to be exposed to threats on the non-breeding grounds and climate change-driven phenological mismatches with their host species on the breeding grounds (Saino et al., 2009;Hewson et al., 2016;Yun et al., 2020). These species and other declining long-distance migrants that winter in Southeast Asia also face stressors from similar sources.
We found that all four species with significant changes in the probability of colonization were affected by climate; both temperature and precipitation in May and June were the most common drivers of dynamic processes of bird occurrence change in our study. Three species had positive effects of climate variables on colonization probability, including breeding season maximum temperature (Eurasian Kestrel, Gray Starting) and increased precipitation (Eurasian Hoopoe). These three species are all short-distance migrants, which have previously been shown to respond more rapidly to changing climate, and are capable of adjusting migration and breeding phenology more effectively than long-distance migrants (Yamaura et al., 2009;Ollerton et al., 2011). Thus, we speculate that these species have potential to colonize vacant habitat as the climate changes. Eurasian Kestrels in Northern Europe have expanded their range and had greater breeding success as spring temperature have increased (Elmhagen et al., 2015;Huchler et al., 2020). Gray Starlings and Eurasian Hoopoe forage on the ground (Joo et al., 2016;Plard et al., 2020), and warmer temperatures and increased precipitation could lead to increased prey availability in agricultural landscapes (Arlettaz et al., 2010;Plard et al., 2020). In general, precipitation during the breeding season is decreasing in Korea (Figure 1C), so any positive effects of precipitation on the probability of colonization by Eurasian Hoopoes may not be widely observed across South Korea. On the other hand, the only species that showed strong evidence for negative effects of warming and increased precipitation colonization was Hazel Grouse. Hazel Grouse is a strictly sedentary species with limited dispersal capacity (Åberg et al., 1995) and a narrow dietary niche. Studies from Europe found that increased precipitation and temperature during the chick-rearing period had negative effects on chick survival (Klaus, 2007). In Galliform birds that have precocial chicks, local climatic conditions, especially precipitation, negatively affect the growth and survival of chicks and population expansion (Viterbi et al., 2015;Terhune et al., 2019). However, it is unclear whether temperature effects in our study were due to challenges associated with thermoregulation (Sunday et al., 2012), the indirect effects of climate change on trophic relationships (Both et al., 2006;Miller-Rushing et al., 2010), or alternative mechanisms. Given neither hunting or poaching have been reported for this species in South Korea, overexploitation can likely be excluded as a cause. Instead, small-scale and subtle habitat changes or food resources that are not afforded by the spatial scale of our study may be implicated. Three species had estimates for local extinction that were influenced by land-cover and climate covariates. Chinese Sparrowhawks almost exclusively forage in shallow freshwater wetlands, especially ricefields surrounded by forests (Kwon and Won, 1975;Choi et al., 2012), and Eurasian Hoopoes also prefer forest edges (Tagmann-Ioset et al., 2012); the covariates on extinction probabilities reflected these associations. Indeed, declines in Chinese Sparrowhawks are likely explained by increased local extinction rates in the many areas where ricefields are being lost (Figures 2, 3). However, it is unclear why Redrumped Swallows show reduced extinction probabilities in areas that warmed during the study period. We speculate that warming spring temperatures could have advanced the emergence of prey species for these swallows, increasing the availability of prey population temporally (Jonsson et al., 2015;Anderson et al., 2019), which could affect colonization of this declining species.
Given that we were only able to detect changes in occupancy with low sensitivity over large sample units, our findings are likely a conservative estimate of the proportion of species in decline and the magnitude of those declines (Strayer, 1999). Decreases in occupancy could only be observed in our study if complete local extinction occurred within a sample unit, so more subtle declines have been obscured. However, our proxy of population decline should be understood with caution when exact changes in abundance are needed, as the relationship between occupancy and abundance can vary among each species and for the density of individuals (Strayer, 1999;Steenweg et al., 2018). In addition, detailed information with diverse environmental metrics, such as vegetation composition, vegetation structure, soil conditions, and riparian management, would greatly improve explanatory climate and land-cover variables we used to model dynamic parameters. Our sampling units may partially obscure associations with these drivers; landcover variables in our study only quantify relatively broadscale changes and may not adequately describe specific habitat requirements for these species (Betts et al., 2014). Future analyses at finer scales of analysis will be possible as data from fineresolution survey units that have been accumulating since 2006. Also unlike our analysis, exactly matching sampling units in future sampling may improve the precision of the estimates, and enable estimation of abundance-occupancy relationships in survey units (Steenweg et al., 2018).
For example, our results indicate that Black-capped Kingfishers have declined severely in South Korea (a 95% reduction in occupancy), which meets the criteria for inclusion as a nationally critically endangered species. This species almost entirely nests in excavated cavities in riverbanks and exposed soils on the hillside, but in the past two decades, intensive riparian management throughout all rivers and streams were conducted for flood management, agricultural irrigation improvement, and as an economic boost (Normile, 2010;Woo, 2010). These projects have removed exposed riverbanks in almost every corner of lowland riparian areas, hardening the riverbanks and altering the flow, removing shallow wetlands on the sides of the rivers and creating deeper channels (Im et al., 2020). Loss of both foraging and nesting habitats for this kingfisher species is potentially related to the drastic decline in their breeding population in South Korea, yet their non-breeding period stressors or even their migration connectivity to wintering grounds have not been identified. To understand better the causes of population decline, it will be critical to understand both the migration ecology and fine-scale habitat requirements of this species.
In addition, the apparent weak effects of land-cover and climate on species' population changes (i.e., extinction and colonization rates) could be due to the presence of other unmeasured stressors. For instance, many declining migratory species spend most of their annual cycle outside their breeding ranges in South Korea. Clearly, we were unable to test for these migration and wintering-ground stressors in our analysis. Understanding demographic processes throughout the full annual life cycle in migratory species is essential for the conservation of migratory species (Marra et al., 2015). For example, direct mortality sources (Loss et al., 2012) are not incorporated into our models. Direct mortality from collision (Bing et al., 2012;Low et al., 2017), poisoning from pesticides and pollutants Barghi et al., 2018), and illegal trapping and consumption (Kamp et al., 2015;Yong et al., 2015) are reported in the region, and are possible cause of declines of these species. These sources of direct mortality, in addition to habitat loss and degradation from land-use change are still major threats on the breeding grounds in South Korea and non-breeding grounds.
Obtaining better knowledge of each species' population trends and potential causes are the first critical steps for biodiversity conservation. Based on our review of the conservation and legal protection status of declining landbird species (Table 1), it is clear that information was lacking for both national and international species status classifications. We recommend that at least seven species from our declining species list should have immediate status re-evaluation and appropriate conservation action to identify and protect remaining breeding habitats. This should be followed by the assessment of major threats, close monitoring of demographic rates, and identification of primary migration routes and key sites for the implementation of conservation measures (McComb et al., 2010). We strongly suggest the development of a transparent and robust plan for monitoring avian biodiversity. Point counts (Ralph et al., 1995), or constant effort surveys, accompanied by explicit recording and documentation of effort and key variables that affect bird detections will greatly improve our ability to infer on bird population changes (MacKenzie and Royle, 2005). Based on our results, we recommend more detailed follow-up studies on individual species as well as the creation and implementation of conservation plans in near the future.
Our study's scope is limited to common and widespread breeding landbird species occurring at relatively low elevations, and specifically the areas where road systems provided access. Even though such surveys are efficient, they are well known to poorly represent the status of bird species that are associated with less accessible habitats (Betts et al., 2007;Harris and Haskell, 2007). Due to the NES's survey protocol, we were not able to incorporate montane forest species, subalpine species, or forest interior species. To overcome this limit we suggest a long-term monitoring plan for this specific stratum of avian habitat in South Korea, which covers more than two thirds of the country.
The current focus on the "Asian Songbird Crisis" tends to be trapping for the cagebird trade in Southeast Asia (Marshall et al., 2020). However, widespread threats such as habitat loss and degradation, as well as climate change across the wintering grounds could be the hidden drivers of population declines in temperate-breeding migratory species that migrate through or winter in tropical regions of the East Asian Flyway (Yong et al., 2015). Without an effort to collect robust information that supports conservation actions, loss of landbird populations and their diversity could pass unnoticed. Even though some recent studies reported changes in avian biodiversity from this region (Amano and Yamaura, 2007;Ko et al., 2014;Tamada et al., 2017;, it still remains a small portion compared to the breeding and wintering range of landbirds in East Asia.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These datasets presented in this study can be found in online repositories. NES data can be accessed from the Environment Digital Library of the Ministry of Environment, South Korea https://library.me.go.kr/, using the keyword " ". The land use and cover maps can be accessed from the National Environmental Geographic Services (https://egis.me.go.kr/).

ETHICS STATEMENT
This study uses existing survey data of wild bird populations. Surveys were done passively, by observing and recording bird abundance and species. Hence no animals were used or manipulated in this study.

AUTHOR CONTRIBUTIONS
HK, C-YC, YM, BM, and MB contributed to the conception, design of the study, and wrote sections of the manuscript. HK and YM collated and prepared the data for analysis. HK performed the statistical analysis and wrote the first draft of the manuscript. All authors contributed to manuscript revision and approved the submitted version.

ACKNOWLEDGMENTS
We are grateful to participants of the surveys, managers, and policymakers who enabled this study for starting national surveys in the early years and maintaining it with substantial improvement in survey designs in the past 10 years. The national atlas data is currently managed by National Institute of Ecology, Seocheon, South Korea. Special thanks to Jungmoon Ha for helping the digitizing of data from printed reports in 2015. Sukhyun Joo and Joonghoon Shin helped to review the code for analysis. Minsu Jeong provided invaluable discussions in preparing the manuscript.