Unravelling migration connectivity reveals unsustainable hunting of the declining ortolan bunting.

In France, illegal hunting of the endangered ortolan bunting Emberiza hortulana has been defended for the sake of tradition and gastronomy. Hunters argued that ortolan buntings trapped in southwest France originate from large and stable populations across the whole of Europe. Yet, the European Commission referred France to the Court of Justice of the European Union (EU) in December 2016 for infringements to legislation (IP/16/4213). To better assess the impact of hunting in France, we combined Pan-European data from archival light loggers, stable isotopes, and genetics to determine the migration strategy of the species across continents. Ortolan buntings migrating through France come from northern and western populations, which are small, fragmented and declining. Population viability modeling further revealed that harvesting in southwest France is far from sustainable and increases extinction risk. These results provide the sufficient scientific evidence for justifying the ban on ortolan harvesting in France.

INTRODUCTION

Unregulated hunting can be a major source of additive mortality in animal populations () and is considered a driver of biodiversity loss () that can lead to large-scale defaunation (), even driving superabundant species to extinction (, ). When illegal killings further occur on top of regulated exploitation, conservation risks can only be exacerbated (). In Europe, the Birds Directive (Directive 1979/409/EEC) is the legal instrument addressing the conservation of all native bird species occurring in the European Union (EU), covers the protection and management of these species and their habitats, and outlines rules for their potential exploitation, including hunting (). A limited number of hunting activities normally prohibited under the Birds Directive (Articles 5 to 8) are permissible by way of constrained derogations (Article 9). Derogations can concern the capture, keeping, and judicious use of certain birds in small numbers under strictly supervised conditions and on a selective basis. Small numbers are any sample of less than 1% of the total annual mortality rate of the concerned population, and the catches must be nondetrimental to the population’s survival ().

The ortolan bunting is a small Palearctic migrant songbird that breeds from Spain to Mongolia and from Iran to northern Finland and migrates to sub-Saharan African winter quarters. Numerous threats (habitat loss, climate change, and illegal captures) are believed to explain the continuous population decline across Europe (), with a drop of 88% since 1980 and a current European population size estimated at 4,750,000 breeding pairs (). Hunters in southwest France have harvested ortolan buntings for traditional uses for decades, taking up to 30,000 individuals annually during fall migration (), despite the species being listed on Appendix I of the EU’s Birds Directive since 1979 (Directive 1979/409/EEC), protected in France since 1999, and red-listed as endangered in many European countries including France (). Conservation organizations have denounced this practice, while local hunters, politicians, and leading French Michelin-starred chefs () defend the persistence of this cultural heritage and gastronomic tradition. French hunters have repeatedly asked for an exception to EU hunting regulations in order to restore legal ortolan bunting hunting in France, requesting an allowable harvest of up to 30,000 individuals, arguing that their catches are negligible compared to the ortolan bunting’s European population size ().

To inform stakeholders, a Pan-European team of researchers conducted a 5-year research project to determine the migratory patterns of the ortolan bunting across Europe (Fig. 1A) and to indisputably determine the origins, numbers, and temporal trends of migrant populations crossing southwest France in autumn. From a biological perspective, the challenge was to determine the migration strategy of a small songbird at the continental scale by combining multiple innovative methodologies. From a political perspective, the research aimed to provide scientific evidence to inform how the Birds Directive should be applied to the ortolan bunting.

Fig. 1

Range, genetic population structure, and migration flyways of European ortolan buntings Emberiza hortulana as revealed by light loggers and stable isotopes.

(A) Top left: Distribution range (orange) during the breeding season (Europe) and winter (Africa). Black triangles locate sites where autumn migrants were sampled. Dots indicate studied breeding populations: dark blue, western genetic cluster; red, eastern genetic cluster; green, northern genetic cluster. (B) Top right: Tracks of 46 migrant buntings obtained with light loggers. Blue, western flyway; green, eastern flyway; orange, putative central flyway. Black dots, INTIGEO loggers; open dots, SOI loggers. (C) Bottom left: Assignment to origin of ortolan bunting feathers using δ2H measurements in feathers collected in populations using the western flyway. Europe: 74 migrating individuals sampled in autumn 2012–2015 in France. Africa: 238 individuals sampled from various breeding populations of the western and eastern genetic clusters in spring 2013–2015. Assignment probabilities were rescaled from 0 to 1 according to the maximum value obtained in a pixel during the assignation process, for each continent, separately. (D) Bottom right: Assignments to origin of ortolan bunting feathers using δ2H measurements in feathers collected in populations using the eastern flyway. Europe: 78 migrating individuals sampled in Kuwait in fall 2015. Africa: 297 individuals sampled in various breeding populations of the eastern genetic cluster in spring 2013–2015. Assignment probabilities were rescaled to range from 0 to 1 according to the maximum value obtained in a pixel during the assignment process, for each continent, separately.

RESULTS AND DISCUSSION

We used three main approaches to investigate migration strategies in ortolan buntings: archival light loggers, stable isotopes, and population genetics. Light loggers are small electronic devices that are affixed to the back of a bird and record light intensity, which enables calculation of the approximate daily positions of a bird (). We retrieved 61 loggers from returning breeding ortolan buntings (table S1), providing 60 complete or partial migration tracks including routes and wintering ground locations (Fig. 1B). Deuterium is a stable hydrogen isotope, which concentrations in feathers (δ2Hf) can be used to probabilistically determine areas where the feathers were grown (). We used δ2Hf of ortolan buntings to identify the probable wintering grounds of breeding populations sampled across Europe (n = 535 individuals) and the probable breeding grounds of migrants captured along fall migration flyways (n = 152 individuals). Last, we genotyped 266 migrant buntings captured in France along the Atlantic flyway to assign these migrants to potential breeding populations, owing to the mapped genetic structure of breeding populations across the continent performed with 24 variable microsatellite loci (n = 575 individuals from 26 breeding sites from Spain throughout Europe to eastern European Russia) ().

The combination of the different complementary approaches identified two migration flyways for the ortolan bunting separating eastern from western European populations. Migration tracks obtained from light loggers illustrate this migratory divide, running from Belarus to the Balkans (Fig. 1B). Bayesian clustering and discriminant analyses uncovered a weak but definitive genetic structure consisting of two main clusters: a northern cluster containing individuals from Fenno-Scandinavia and part of the Baltic states and a southern cluster, further longitudinally subdivided into an eastern cluster (Russia, Belarus, Serbia, and Greece) and an admixed western cluster weakly differentiated from the eastern one (from Poland to Spain; Fig. 1A). Deuterium concentrations in scapulars/body coverts (molted in Africa during the winter) of breeding/spring migrating birds revealed two distinct groups of countries representing the flyways: an eastern group (Belarus, Russia, and Serbia, together with spring migrants captured in Israel and Kuwait) and a western group (breeding birds from Finland, France, Lithuania, Poland, and Spain) (fig. S1). These two groups correspond to individuals having molted in similar isotopic environments, probably in the same geographical areas. Geographical assignments to the winter range (Fig. 1, C and D) show that these two groups correspond to a segregation between birds wintering in East Africa (eastern flyway) and birds wintering in West Africa (western flyway).

Within the restricted estimate of 469,000 pairs (range, 354,500 to 619,500) breeding west of the migratory divide (), only ~81,000 pairs (range, 46,000 to 116,000) and their offspring are estimated to use the Atlantic route of the western flyway in autumn, and thus visit southwest France (data S1), while the other western ortolan buntings migrate along the Mediterranean coasts. With an average fledging success rate of 1.72 juveniles per pair (), the fall flux in southwest France (the Atlantic route) should consist of an average of 300,000 migrants, including 138,000 fledged juveniles and 162,000 adults. Stable isotope and genetic assignments to origin of migrants revealed that, among these numbers, ~65 to 75% originate from Western Europe, mainly Poland and Germany, while the rest originates from the Baltic countries and Fennoscandia (Fig. 1C and figs. S3 and S4). We estimated the recent trends (2000–2014) of the corresponding breeding populations as a decline of 20 to 30% (Supplementary Materials and data file S1). This decline represents an annual decrease of 1.7 to 2.7%, translating to 2750 to 4400 mature birds failing to return to breed in the spring every year. Identifying the genetic structure of European breeding populations confirmed that ~30% of the migrant ortolans sampled in France during autumn migration belong to the northern genetic cluster (fig. S4). This northern cluster is genetically isolated with almost no emigration or immigration from Western or Eastern Europe (), and the breeding populations there have declined by almost 5% annually (53% from 2000 to 2014; data file S1).

High-quality and updated scientific information on population size () and demographic parameters including mortality are a prerequisite for reliable calculations of the “small numbers” (1%) concerned by a potential constrained derogation allowable under the Birds Directive (). Monitoring of a ringed breeding population in Norway provided estimates of adult and first-calendar-year survival probabilities and reproductive success (table S2) (, ). One percent of the annual mortality of postbreeding numbers visiting southwest France represents ~1700 (1518 to 1922) individuals.

Population viability analyses confirmed that current northern populations of ortolan buntings are directly threatened with extinction and could not persist without marked increases in survivorship, with prospective decreased extinction risk associated with enhanced habitat quality improving reproductive success (Fig. 2). When parameterizing models with the demographic parameters observed in Norway (table S2), the population is deterministically decreasing (growth rate, λ = 0.83) and goes rapidly to extinction (median time to extinction, 23 years). On the basis of criterion E (quantitative estimate of viability) of the International Union for Conservation of Nature (IUCN) Red List applied at the population scale (), the population would qualify as “Endangered,” facing a high risk of extinction in the near future. Assuming most optimistic scenario in terms of operational sex ratio (three females for four males, SR = 0.75) and clutch size (CS = 4.9) does not result in a viable population (λ = 0.96; time to extinction, 52 years) and qualifies the species for a “Vulnerable” status, with a high risk of endangerment in the medium term. Improvements of only the sex ratio and the habitat quality/clutch size are not sufficient to improve population viability, and all scenarios indicating a reduced extinction risk require substantial improvements in survival.

Fig. 2

The 100-year extinction risk of northern ortolan bunting breeding populations with respect to various scenarios of increased survival probability (decreased illegal harvesting).

Different scenarios of increased habitat quality, and hence breeding performance (increased clutch size and more balanced sex ratio) are shown with different symbols, based on population viability analyses conducted with vital rates drawn from population dynamic research carried out in Norway (triangles; clutch size of 4.25 eggs, sex ratio with one male out of two finding a mate) and in Finland (squares; clutch size of 4.57 eggs, sex ratio with two males out of three successfully pairing), while further more optimistic scenarios are shown with circles (first with a larger clutch size of 4.9 eggs and then also with a more balanced sex ratio, with three of four males successfully pairing, as would happen in populations hypothetically well connected within a metapopulation system). Increased survival probability is shown in terms of the number of additional individuals (first and second values of the x axis legend) surviving the autumn migration (second value) or until the next spring (first value). The red rectangle encompasses the range of known and/or suspected harvested numbers in southwest France (estimated between 15,000 and 30,000 ortolans per autumn). The colors depict the IUCN Red List status of the population as estimated from maximum extinction probability after 10, 20, and 100 years and three and five generations (see Materials and Methods for details; green, Least Concern; orange, Vulnerable; red, Endangered). In the best demographic scenario, the extinction risk is predicted to drop from 100 to 66% if 15,000 additional buntings survive the fall or to only 15% if 30,000 additional buntings survive the fall migration. Photo credits: S. Minkevicius.

Assuming some absolute decrease in mortality probabilities leads to notable reductions of extinction probabilities and improvements of the population conservation status in all cases (Fig. 2 and fig. S5). However, with the default Norwegian demographic parameters (table S2), only strong reductions in mortality can significantly decrease the risk of extinction at a 100-year time horizon. A decrease of 10 percentage points in mortality probabilities, i.e., ~70,000 additional individuals surviving the autumn migration, can make the 100-year extinction probability drop from 100 to 17% only. For populations in better fate, such as Finnish populations (larger clutch size and better balanced sex ratio; Fig. 2 and fig. S5), mild to moderate reductions in mortality probability lead to strong changes in extinction risk and improvements in conservation status. An increase in survival representing an additional 30,000 postbreeding live birds (i.e., a reduction of approximately four points of the mortality percentage) would induce a 23% decrease in the 100-year extinction probability, while this decrease would reach 50 to 85% in the most favorable scenarios of reproductive output, mimicking conditions that occurred in northern Europe before the attested large-scale decline of farmland birds (). In these favorable scenarios, a moderate reduction of mortality (five points) is sufficient to improve the Red List conservation status of the population from Endangered to “Least Concern.” The models illustrate that harvesting in France potentially greatly affected the fate of northern populations when breeding conditions were better and are still responsible for nearly a quarter of current extinction risk of Finnish-like populations. French hunting is therefore identified as being partly responsible for the continuous decline of the ortolan populations using the western Atlantic flyway.

Together, this information diminishes the credibility for a regularized and sustainable harvest of the endangered ortolan bunting. Moreover, the 15,000 to 30,000 individuals recurrently harvested yearly in France () largely overpass the 1% of the annual mortality of postbreeding numbers visiting southwest France (estimated as 1518 to 1922 individuals). Hence, with this study, French authorities now have rigorous scientific data to make an informed decision to conclusively ban ortolan bunting hunting, actively police poaching, and increase the chances of the ortolan bunting to survive global change.

MATERIALS AND METHODS

Breeding populations

We used updated national population sizes and their recent trends (), benefitting also from the recent reporting of EU members on Article 12 of the Birds Directive. Further updates were necessary not only for non-EU members but also for Poland, which previously published erroneous figures [see the Supplementary Materials and (–)]. We combined these national population sizes with data on isotopic and genetic flyway assignments to estimate the total numbers of individuals using each flyway. For routes along the western flyway, we proposed three scenarios based on the proportions of logger-equipped males of each country using the Atlantic or the Mediterranean route (see the Supplementary Materials and data file S1). Field work included capturing males on their breeding territories, sampling one tail feather and three to six body coverts (scapulars), and equipping a sample of males with light loggers. We used DNA extracted from bases of tail feathers to characterize the genetic structure of the breeding European populations. Body coverts molt twice a year, first during the complete summer premigration molt and then during the winter before spring migration. Hence, breeding birds have body coverts grown on their winter molting grounds so that the deuterium concentration in these spring body feathers was used to assign breeding individuals to potential wintering grounds.

Sampling along migration flyways

We also sampled fall migrants along the western and eastern flyways to collect tail feathers. In autumn, tail feathers either have grown in the nest for first-calendar-year birds or have been molted on the breeding grounds after reproduction and before fall migration for older birds. Hence, fall migrants have tail feathers grown on their breeding grounds, so we used the deuterium concentration in these summer-grown tail feathers to assign fall migrants to potential breeding grounds. We sampled migrants in France in August to September of 2012 to 2015 in two ways: first, wild migrants were captured, ringed, and released by ringers, and second, wild migrants were captured by hunters but seized by police, placed in a wildlife rescue center, and sampled at ringing before release into the wild. Seized birds comprised both wild birds recently captured and live dummies kept in captivity for 1 year or more. Seized birds were sampled (one tail feather; n = 274) when being ringed at the care center just before being released in the wild. We divided feather data from these seized birds in three groups: (i) first-calendar-year individuals (n = 34; aged using a combination of molt contrasts, body feather streaking, and characteristic pointed tail feathers; these birds hatched in the year, so in September, they have retained feathers grown in the wild in the summer at their hatching site), (ii) attested decoys [n = 21, feathers grown in captivity, including all individuals with coloration anomalies (one or more white secondary or primary, fuliginous body coverts; this sample also includes one male kept in a care center for more than 1 year)], (iii) and unknown status (n = 179; which should be a mix of nonobvious decoys and wild recently captured adults). Whatever the status of all individuals, they all have hatched in the wild (there is no captive breeding of the species), but their feathers have grown either in the wild on their hatching/breeding grounds (if recently caught) or in captivity (for nonobvious decoys).

Along the eastern flyway, we captured fall migrants in Kuwait in September 2014 and spring migrants in Kuwait and Israel in April 2015. We collected a tail feather from each bird, first to obtain DNA to assign migrants to a geographically defined genetic breeding cluster. We also used deuterium concentration in the summer-grown tail feathers to assign fall migrants to potential breeding grounds. Spring migrants have body coverts molted on the wintering grounds, so spring migrants captured in Israel and Kuwait also provided additional cover feathers for identifying potential wintering grounds of birds following the eastern flyway.

Stable isotopes

Deuterium concentrations (δ2H) have a latitudinal structure across some continents that are linked to their concentrations in precipitation (). Bird feathers grown at a given place have δ2H values (δ2Hf) reflecting amount-weighted mean growing season precipitation at sites where feathers were grown. Hence, by measuring δ2Hf, it is possible to infer a probabilistic geographical space where this feather might have grown. We used a spatially explicit likelihood assignment method to delineate probable origins for ortolan buntings by converting an amount-weighted growing season precipitation surface () to a feather isoscape using a calibration equation developed for Eurasian reed warblers (Acrocephalus scirpaceus; δ2Hf = −10.29 + 1.28* δ2Hp) (). The residual SD [SD = 10.36 per mil (‰)] from the linear regression model used to calibrate the precipitation surface for Eurasian reed warbler feathers was included in the assignments as an estimate of error. As the ortolan bunting also forages in agricultural habitats, we did not consider feather carbon (δ13C) or nitrogen (δ15N) isotopes because their concentrations may be linked to local agricultural inputs.

Feathers were cleaned in 2:1 chloroform:methanol solvent rinse and prepared for δ2H analysis at the Stable Isotope Laboratory of Environment Canada, Saskatoon, Canada. The δ2H of the nonexchangeable hydrogen of feathers was determined using a method based on two calibrated keratin hydrogen isotope reference materials [CBS (Caribou Hoof Standard) and KHS (Kudu Horn Standard)]. We performed hydrogen isotopic measurements on H2 gas derived from high-temperature (1350°C) flash pyrolysis (EuroVector 3000; Milan, Italy) of feather subsamples (350 ± 10 μg) and keratin standards loaded into silver capsules. We analyzed resultant separated H2 on an interfaced Isoprime (Crewe, UK) continuous-flow isotope-ratio mass spectrometer. Measurement of the two-keratin laboratory reference materials corrected for linear instrumental drift were both accurate and precise with typical within-run mean δ2H ± SD values of −197 ± 0.79‰ (n = 5) for CBS and −54.1 ± 0.33‰ (n = 5) for KHS.We report all results for nonexchangeable H expressed in the typical delta notation, in units of per mil, and normalized on the Vienna Standard Mean Ocean Water–Standard Light Antarctic Precipitation standard scale.

We limited assignments to origin to the known range of the species, where feathers could have grown or molted. We obtained the digital distribution map of ortolan buntings from the BirdLife International and NatureServe () and modified it to include known breeding areas in Kazakhstan (see map at www.birds.kz/v2taxon.php?s=577&l=en). Ortolan buntings molt their feathers either on the breeding or wintering grounds; therefore, assignment to origin analyses only included the species’ breeding or wintering range, respectively.

To determine the wintering grounds of ortolan buntings from various breeding populations, we analyzed δ2Hf values from birds captured on their breeding territories. The ortolan bunting molts body coverts twice each year: before the autumn migration on the breeding grounds (August) and before the spring migration on the wintering grounds. Hence, breeding birds in May to June have body coverts molted on their African wintering grounds.

We analyzed samples (scapulars or body coverts) obtained from territorial males captured on their breeding grounds in May to June of 2013 to 2015 in various countries (sample sizes given in brackets): Belarus (n = 43), Finland (n = 139), France (n = 34), Lithuania (n = 43), Poland (n = 31), Serbia (n = 7), Spain (n = 19), and Russia (n = 192). We also analyzed body feathers of migrants captured in spring (April 2015) in Kuwait (n = 45) and Israel (n = 10) to similarly assign individuals to potential wintering grounds. Deuterium concentrations in these scapulars/body coverts (molted in Africa during the winter) revealed two distinct groups of countries representing an eastern flyway: (i) Belarus, Russia, and Serbia, together with spring migrants captured in Israel and Kuwait, and a western flyway: (ii) breeding birds of Finland, France, Lithuania, Poland, and Spain (fig. S1). Each group corresponds to individuals having molted in similar isotopic environments, probably in the same geographical areas.

To determine the potential geographical origins of ortolan buntings captured across France during autumn migration, we measured δ2Hf from live birds captured during migration and released in the wild. In this species, autumn birds either have tail feathers grown in the nest (for juveniles) or molted on the breeding grounds before the fall migration (for adults). We therefore used samples collected from migrating individuals in France in August to September (in 4 years, 2012–2015; n = 40) and samples collected from seized first-calendar-year individuals (n = 34), while we also did the same analysis with samples collected along the eastern flyway in Kuwait in September 2014 (n = 78).

In separate assignments to origin using δ2Hf, potential origins for birds caught along the western (n = 74) and eastern (n = 45) flyways were restricted to their respective breeding areas on either side of the longitudinal migratory divide following political borders, based on results from geolocation data and genetic structure of breeding populations. We used an odds ratio of 2:1 to assign each feather sample to potential origins, where cells in the isoscape in the upper 67% of probabilities were considered as likely (1) origins and all others were considered unlikely (0). We first conducted assignments by country in which birds were captured, and then, we grouped samples based on statistical analyses, grouping national populations with similar δ2Hf values [country effect in General Linear Models (GLMs); see, e.g., fig. S1]. Assignments resulted in a spatially explicit binary surface for individual birds, which we subsequently summed across assignments for all individuals within a group to represent potential origins for that group. For each pixel in the assignment raster, the presented value is the proportion of all sampled individuals, which could potentially have grown the analyzed feathers within the pixel, rescaled to range from 0 to 1.

To determine the latitudinal origin of migrants along the eastern flyway, we performed an assignment to origin analysis on samples collected in September 2015 in Kuwait from wild birds captured with mist nets while migrating. A preliminary assignment was conducted using the 78 individuals captured in Kuwait to the whole breeding range of the species. When restricting the assignment to countries with populations using the eastern flyway (i.e., wintering in East Africa, as defined by the genetics and light logger studies; Fig. 1D), the possible geographical origin of these migrants appears as being mainly southern Russia, southern Ukraine, and northern Romania, as well as Turkey. Given the respective population sizes breeding in these countries, we can consider that most migrants captured in Kuwait likely have origins in Russia, as almost all Russian ortolans (2 million to 4.3 million pairs) breed within the pale to dark blue areas depicted in Fig. 1D. From the assignment surface analysis, Turkish (0.5 million to 1 million pairs) and Romanian (225,000 to 550,000 pairs) populations probably also use this flyway, while Ukraine holds only 58,000 to 67,000 pairs ().

To determine the latitudinal origin of migrant ortolan buntings captured in August to September of 2012 to 2015 along the western Atlantic flyway (in France), we performed an assignment to origin analysis on samples collected by ringers on wild birds captured in France. We derived a preliminary assignment depiction from the 74 individuals captured by ringers or seized by the police from caged birds and aged as first-calendar-year birds. When restricting the assignment to countries with populations using the western flyway (i.e., wintering in West Africa; Fig. 1C), the possible geographical origin of these migrants appears as being mainly Poland and Germany and, to a lesser extent, France, the Baltic States, and southern Sweden. Very few individuals appear to come from Norway (which has a restricted breeding population), Finland (which may use a more continental route), or northern Russia (which uses the eastern flyway and has a very low population density). Again, it is possible that such a map—restricted to the countries with populations using the western flyway—mirrors local population densities, as the largest populations in the area of concern are located in Poland. The only exceptions here are Finland, which holds 7000 to 19,000 pairs, and Sweden, where population sizes are mostly smaller than in Poland, with a very low probabilistic contribution to migrants captured in France.

Archival light loggers

Light loggers are small electronic devices that record light intensity, which, if retrieved after a migration cycle, enables the calculation of the approximate position of the logger given the duration of the day and the time of solar noon, both depending on latitude and longitude for a given calendar day. Geolocators do not remotely transmit positions; therefore, tags must be retrieved to download data. This necessitates animal recapture and, thus, site fidelity at wintering, breeding, or stopover locations. We used data collected by SOI-GDL 3.0 (Swiss Ornithological Institute, 0.65 g) and INTIGEO P65C2-7 (Migrate Technology, 0.74 g) loggers, fixed with an ultraviolet-proof leg-loop string harness to the bird’s back. The device and associated attachment material represent 3.1 to 3.5% of a bird’s weight, a male ortolan bunting weighing on average (±SD) 21 ± 2 g during the breeding season. Colleagues from Sweden also successfully deployed Intigeo P65C2-7 loggers and previously published a large part of their tracking results (). We reanalyzed the raw data of tracks from nine individuals, corresponding to 11 migration tracks (two individuals were tracked during two migration cycles); the tracks obtained in 2015 were not published previously. We reanalyzed these data with similar parameterizations as for the other data. We also successfully deployed 18 Intigeo P65A9-10 longpin loggers in Germany. They provided tracks of 17 individuals, corresponding to 19 migration tracks (one individual fit in 2013 and in 2014 and another individual fit in 2013 and retrieved in 2015). Two loggers stopped recording light data during the winter, two failed during the spring migration, but all provided autumn tracks.

Overall, we retrieved 61 loggers (table S1), 57 of which recorded data: Belarus (n = 1), France (n = 1), Russia (n = 14), Lithuania (n = 1), Finland (n = 8), Norway (n = 1), Poland (n = 2), Spain (n = 2), Sweden (n = 9), and Germany (n = 18). Three loggers (two from Sweden and one from Germany) provided data for two annual migration cycles, so we obtained 60 available autumn migration tracks. All retrieved INTIGEO loggers contained complete data for the period carried. Eleven SOI loggers were recorded from May until early August to early October 2015, so they provided partial migration tracks. Four more SOI loggers failed within a couple of weeks. We used the GeoLight 2.0 package in R (https://CRAN.R-project.org/package=GeoLight) to analyze light data using the threshold method (). Light-level geolocation data of the loggers used in Germany were analyzed using the freely available R packages SGAT [Solar/Satellite Geolocation for Animal Tracking; GitHub Repository (http://github.com/swotherspoon/sgat)] and “BAStag” (basic data processing for light based geolocation archival tags; R package).

Genetics

The aim of the genetic part of the research program was to determine the genetic structure of European breeding populations [see () for additional details] and to assign autumn migrants captured on the western flyway (southwest France) and the eastern flyway (Kuwait and Israel) to genetic breeding groups/areas.

We used a shotgun sequencing approach on an ion PGM (Personal Genome Machine) platform (Life Technologies) to develop 24-microsatellite loci from muscle tissue obtained from one specimen found dead in Kuwait in September 2014. We extracted DNA from the tail feather collected on all individuals captured on breeding sites or along migration flyways. A total of 1127 samples, including 143 duplicates, were genotyped for the 24 loci. The resulting breeding dataset consisted of 575 individuals sampled at 26 breeding sites for population-level analyses. We carried out further analyses without three loci that displayed a high frequency of null alleles and deviated from the Hardy-Weinberg equilibrium (HWE). We performed a randomized G test with 1000 replicates to test for genetic differentiation among sites.

We used two methods to uncover genetic population structure. First, the Bayesian clustering program STRUCTURE () allowed the assignment of the 575 breeding individuals to K populations by minimizing deviations from the HWE. We obtained the optimal number of clusters K from ΔK, based on the rate of change in the log probability of data in successive K values (). Second, we applied a discriminant analysis of principal components (DAPC) () to the breeding dataset. DAPC is free from population genetic assumptions with inferences made on allelic similarity. It summarizes genetic variability of individuals within groups while optimizing group discrimination. We used sampling sites as the grouping variable. We carried out analyses out using the ADEGENET 2.0.1 package in R 3.3.1. We retained the first 115 principal components in the data transformation step, corresponding to 84.2% of genetic variance, and we saved three discriminant functions for further analyses. We reached consensus on population structure based on results from these two methods. The magnitude and direction of contemporary gene flow occurring between the consensus populations were estimated using the program BAYESASS 3.0.1 ().

We then performed a cross-validation of the population structure. We randomly split the full breeding dataset (575 individuals) into a training and a validation dataset by randomly assigning 70% of individuals from each sampling site to the training set (402 individuals) and the remaining 30% (173 individuals) to the validation set. The training set defined the genetic makeup of the clustering to be tested, and we assigned individuals from the validation set to one of these populations by the program GENECLASS () using the Bayesian method described in (). We repeated the process 10 times. The cross-validation on the three populations supported the strength of the northern and eastern clusters with an average of 85 and 79% of individuals, respectively, correctly assigned (fig. S2). Moderate correct assignment to the western population (52%) with a large contribution of the eastern population (36%) suggests low differentiation between these populations.

We finally assigned the 396 individuals captured during migration along the western flyway (in southwest France) and the eastern flyway (Kuwait and Israel) to one of the clusters defined previously by the program GENECLASS using a Bayesian method (). We first performed assignments of migrants to the three defined genetic clusters (western, eastern, and northern) to reflect both the higher level, and thus stronger, genetic structure and our knowledge of migratory flyways.

We grouped individuals captured in southwest France into three categories, as for isotopic analyses: (i) wild migrants (captured in the wild with mist nets by ringers and first-calendar-year individuals seized by the hunting police, hatched and grown the same year in the wild), (ii) dummies (adults obviously kept in captivity during their last molt because of aberrant colors, e.g., blackish plumage, or one or more totally white remiges), and (iii) status unknown (adult seized birds that could be dummies or wild migrants). Chi-square tests were used to compare the distributions of origin assignment and especially to compare wild individuals (n = 73), captive individuals (n = 21), and individuals of unknown status (n = 172). The number of wild and first-calendar-year individuals here is 73 and in isotopic analyses is 74, because we did not succeed in sequencing one of the 74 samples. We found no significant difference in the distribution of assigned individuals to breeding populations among the different categories of French migratory birds (χ2 test = 16.601, P = 0.165). For genetic analyses, we thus pooled all birds captured in France to represent the western flyway.

We found significant differences in assignments between individuals captured in the eastern and western flyways (χ2 test = 30.565, P < 0.001). Hence, the procedure assigned 67.69% of migratory birds using the eastern flyway to the eastern population and 23.84% to the western population (fig. S3). In contrast, the assignment of individuals from the western flyway was more equally distributed among populations, with 39.10% assigned to the western population, 38.72% to the eastern one, and 21.80% to the northern one.

Stable isotopes and light loggers data indicated that individuals caught in southwest France did not originate from the eastern population. We therefore reran the assignment for the western flyway and removed the eastern population as a potential breeding origin. Individuals originated mostly from the western (66.54%) and northern (33.08%) populations (fig. S4). Further details on the genetic study can be found in (), including pairwise FST (fixation index) between populations and estimates of contemporary gene flow.

Demographic parameters

Demographic data are based on the Norwegian ortolan bunting population monitored during the breeding period (May and June of each year) in central Hedmark County (60°29′ to 60°53′N and 11°40′ to 12°18′E) and Akershus County between 1996 and 2005 [e.g., (, )]. This is the last breeding population in Norway, which decreased from approximately 250 singing males in the 1990s to only 14 males in 2016. Juvenile and adult male birds were captured and fitted with a metal ring and colored rings, allowing individual identification of adult males (). In each year, field observers recorded recaptures (mainly visual) of already marked birds, including information on their breeding status (paired or unpaired). We did not estimate survival for females because the sample size was too small.

We analyzed 419 adult male capture-recapture histories using multievent capture-recapture models, with the program E-SURGE 1.9.0 (). We estimated three basic types of parameters: the initial state probabilities (from the Π-vector, following E-SURGE’s notations), the between-state transition probabilities (probabilities of transition between states and survival from the Ψ and Φ matrices respectively), and the event probabilities P (from the B matrix). We subdivided the B matrix into two matrices, the P matrix for recapture probabilities and the A matrix for probabilities to confirm the uncertain breeding states when detected. We modeled the between-state transitions (Ψ) conditional on survival (Φ). We estimated parameters simultaneously by standard maximum likelihood procedures, and models were ranked using Akaike’s information criterion (AIC) adjusted for small sample size AICc ().

Using this framework, we examined the effects of reproductive status (paired versus unpaired males) on male ortolan bunting survival Φ, breeding transition Ψ, and recapture P and state assignment A probabilities. On these parameters, we assessed the effect of time (t), as well as the effects of the breeding state of the individual before the current recapture (f) and during the current recapture (to). Parameter constancy is denoted with “i.” The combination of these effects in a given parameter could be either additive (+) or nonadditive (.). We achieved the model selection by starting from model Φ(f.t) Ψ(f.to.t) P(f.t) A(f) and using a sequential approach to model simplification. We assumed in all cases that the probabilities of assignment A of recaptured individuals to the various states were state dependent [i.e., A(f)]. We began the model selection by testing temporal and state effects and their interactions on P; then, we used the most parsimonious P model to model Ψ and, finally, Φ.

Because there is currently no general procedure for goodness-of-fit (GOF) tests of multievent models, we estimated the GOF of a simplified dataset with known breeding states. We considered the Jolly-Move (JM) model as the starting general model (). We used the program U-CARE 2.3.2 () to perform GOF tests and assess data conformity to JM.

GOF results indicated that data conformed to the JM model (global test, P value = 0.565 and c-hat = 0.96). On the basis of the model selection procedure described in methods, the most parsimonious model was Φ(i) Ψ (f) P(f + t) A(f), in which survival was constant over time and similar among breeding statuses (table S3). The survival probability averaged 0.64 (SE = 0.02) for adult males regardless of their breeding status.

Population viability models

As Norwegian-breeding ortolan buntings migrate via the Atlantic route of the western flyway, much of this population is likely to be exposed to hunting pressure in southwest France, a finding further supported by multiple ring recoveries (see the Supplementary Materials). By developing a population viability framework based on the Norwegian demographic default parameters (table S2), we illustrated the relative impacts of habitat quality (clutch size and sex ratio) and harvest on the predicted extinction risk of this population within the next 100 years. We used our estimate of 0.64 ± 0.02 for the annual survival probability of adult males; we assumed female adult survival to be equal to male survival. We retrieved juvenile survival and all breeding parameters from (). Table S2 reports all default demographic parameters implemented in population viability models.

In Norway, the ortolan bunting population is suffering from a male-biased sex ratio in adults, which is one of the factors limiting its growth (). This deficit in females is unlikely to result from a skewed sex ratio at birth or from sex-biased mortality. Rather, it is presumably due to female-biased natal dispersal () and thus largely results from the isolation of the Norwegian population (absence of female immigration). We modeled the deficit of females as the proportion of breeders (PB) in a given year () using the empirical estimate of 0.52 (i.e., 52% of adult males become paired and can breed). Because this low PB likely results from the isolation of the Norwegian population, we considered two more optimistic scenarios with PB = 0.66 and PB = 0.75 to model potential dynamics in less or nonisolated populations. In the main simulations, we assumed a mean clutch size of CS = 4.25 based on the empirical distribution provided in () for Norway. We further considered two additional scenarios of larger clutch sizes, corresponding to more optimistic scenarios of a higher-quality habitat yielding a higher clutch size: CS = 4.57 as observed in Finland () and CS = 4.9 as observed in Finland in the 1960s (), both implemented as truncated Gaussian distributions. By developing scenarios using these former enhanced demographic parameters, we aimed at testing (i) the potential former impact of harvesting on population fate or (ii) the potential impact that harvesting could have if the necessary changes in agricultural policies succeed in restoring favorable habitats and breeding conditions for currently declining farmland birds.

To assess the effects of some reduction of harvest related mortality, we assumed scenarios of absolute (not relative) reductions of mortality probabilities of both adult and juvenile birds. We implemented the enhanced survival as a reduction from 0.25 to 10% of mortality, corresponding to ca. 750 to 30,000 more buntings surviving to the next spring (in other terms, 1750 to 70,000 more buntings that will not die during the first autumn). These figures should include putative illegal autumn kills, estimated between 15,000 and 30,000 (), although numbers should be moderated by the proportion of the harvest mortality, which is not additive.

We conducted a population viability analysis using Vortex 10.2.7.0 (). Starting from a population of 400 individuals (which was the size of the ringed Norwegian population when the ringing study started) at the stable age distribution and with balanced sex ratio, we estimated the probability of population extinction (PE) under various ecological scenarios (see below) and over several time horizons. We considered some reduction in demographic parameters (juvenile survival and fecundity) due to inbreeding, assuming that 50% of inbreeding depression was due to recessive lethal mutations () and based on the generic estimate of 6.29 haploid lethal equivalents (). We included both environmental stochasticity and catastrophes. We considered environmental stochasticity based on specific estimates of interannual SD in demographic rates provided by (). We included the effects of demographic and environmental stochasticities, catastrophic events, and inbreeding depression. We considered catastrophes based on generic estimates (), assuming a 90% reduction in demographic parameters (survival and breeding success) occurring randomly with a 0.01 per generation frequency (). We transformed the per generation frequency of catastrophic events to an annual frequency based on the generation time and implemented in Vortex. We set the population carrying capacity to 2000 individuals in all scenarios considered.

In all scenarios, we computed extinction probabilities after 10, 20, and 100 years and after three and five generations to estimate the conservation status of populations based on the E criterion (quantitative estimate of viability) of the IUCN Red List framework (). In figures, the 100-year extinction probability is provided for each scenario, and the associated Red List status is indicated (by the color of the symbol; Fig. 2 and fig. S5). Statuses were determined according to the following rules:

1) Populations with PE < 0.1 after 100 years were categorized as Least Concern (green symbols).

2) Populations with PE ≥ 0.1 after 100 years and with PE < 0.2 after 20 years or five generations (whichever was the longer) were categorized as Vulnerable (orange symbols).

3) Populations with PE ≥ 0.2 after 20 years or five generations and PE < 0.5 after 10 years or three generations (whichever was the longer) were categorized as Endangered (red symbols).

In all scenarios, extinction probabilities were computed from Monte Carlo simulations based on 1000 iterations.

Converting mortality rates to bunting numbers

Scenarios of mortality reduction

Applied to all breeding populations using the Atlantic route of the western flyway, our scenarios of mortality reduction would represent 0 to 10% of 300,000 individuals, hence 0 to 30,000 individuals. The postbreeding population size visiting southwest France during the fall migration is indeed estimated at 300,000 individuals. This corresponds to 81,000 pairs (see the Supplementary Materials and data file S1), hence 162,000 adults with an annual survival probability of 0.64 ± 0.071 (table S2). The remaining birds are 138,000 first-calendar-year birds (deduced from the average of 1.7-fledged chicks per pair) () with a survival probability of 0.176 ± 0.063 (table S2).

If 30,000 birds are hunted in France, and if the capture probability to be hunted is not age dependent, then these should be, on average, 16,200 adults—of which 16,200 × 0.64 = 10,368 would have survived until the next spring—and 13,800 juveniles—of which 13,800 × 0.18 = 2484 would have survived until the next spring. Hence, the taking of 30,000 buntings in autumn would induce a reduction in size of the breeding population by 12,852 individuals the next spring. This is the way we estimated the potential impact of numbers harvested in autumn on extinction risk for each step of reduction in mortality probabilities and reported as harvested numbers in Fig. 2 (with two numbers reported on the x axis legend: the upper value is the number of harvested buntings and the lower value is the number of those harvested buntings, which would have survived until the next spring). These numbers are certainly overestimates as they are based on the hypothesis that mortality would not affect numbers between fledging and migration over France, which is obviously not the case for juveniles, although also for mature birds. The true numbers are therefore probably between those estimates and the total numbers represented by the reduction in mortality probabilities.

Small numbers

As the average annual mortality probability is 0.36 ± 0.07 for adult ortolan buntings and 0.82 ± 0.06 during the first year of life (), we deduced that 1% of the annual mortality (lower-upper values) of postbreeding numbers (300,000 individuals, i.e., 162,000 adults and 138,000 first-calendar-year individuals) would represent ~583 (468 to 698) mature and ~1137 (1050 to 1224) first-calendar-year individuals, for a total of ~1720 individuals (1518 to 1922).