Feather samples from BWTE were collected during the autumn migration and wintering seasons in Cuba. During the autumn migration of 2021, samples were collected on ten different dates, from September to November, in western Cuba. The sampling area was Bajos de Santa Ana, situated west of Havana on the north coast of Cuba (23°04’06’’ N, 82°31’33’’ W). This location is a known stopover site for waterfowl during migration and consists of a mangrove forest covering approximately 21.9 hectares. The site experiences significant human activity due to its proximity to the town of Santa Fe. Presently, because of habitat degradation and hunting pressure, it serves as a stopover site for only a brief period. These characteristics ensured that the samples were obtained from individuals exclusively engaged in migration.

During the winter season, feathers were collected from the Southern Sancti Spíritus Wetland, located in the central region of Cuba (21°38’09’’ N, 79°12’20’’ W), in December 2020. This area is recognized as an Important Bird Area. The majority of this wetland region is occupied by rice fields, covering approximately 27,000 hectares (Acosta et al., 2002). Separated from the sea by a coastal belt of wetlands containing mangroves and several lagoons (BirdLife International, 2020) to the west, it adjoins a group of lagoons and intertidal mudflats that are part of the Tunas de Zaza Fauna Refuge. The entire coastal area serves as a vital refuge for migratory birds, including those from the Charadriiformes and Anseriformes orders. It hosts large congregations, sometimes numbering up to 100,000 individuals of BWTE (BirdLife International, 2020).

All sampled individuals were from carcasses obtained from the hunting season in Cuba. The right wing of each individual was taken, and their sex (Female/Male) and age (After Hatch Year: AHY and Hatch Year: HY) were identified based on wing plumage characteristics following the methods described by Pyle (2008). The primary flight feathers were added to the feather collection at the Felipe Poey y Aloy Museum, University of Havana. For isotopic analysis, the most distal primary flight feather (P10) was selected.

Primary feathers (P10; n = 276) were analyzed for δ ²H at the Laboratory for Stable-Isotope Science – Advance Facility for Avian Research (LSIS-AFAR, Western University, London, Ontario, Canada). Feathers were cleaned of surface oils using a 2:1 methanol:chloroform solvent (>12h) and allowed to air dry (24h). Feather vane was clipped from the distal end of the feather and a subsample (0.35 ± 0.02 mg) packed into silver cups and crushed. Samples were loaded, simultaneously with laboratory standards, in a Uni-Prep autosampler (Eurovector, Milan, Italy), heated to 60°C and flushed with dry helium and maintained under helium pressure (Wassenaar et al., 2015). Samples were combusted using flash pyrolysis (~1,350°C) on glassy carbon, separated via a Eurovector 3000 elemental analyzer interfaced with a Thermo Delta V Plus continuous-flow isotope-ratio mass spectrometer (CF-IRMS; Thermo Instruments, Bremen, Germany). Feather δ ²H values were derived using the comparative equilibration method of Wassenaar and Hobson (2003). Two keratin standards (CBS, δ ²H = −197‰; KHS, δ ²H = −54.1‰) were included at the beginning (n = 3, per standard) and end (n = 2, per standard) of each run (n = 38 samples per tray). All values were corrected for linear instrumental drift. At LSIS-AFAR, across-run error has been estimated to be ±2‰. All values are reported in delta (δ) notation in per mil (‰) deviations from the Vienna Mean Ocean Water standard (VSMOW).

The effect of season, date of migration, sex and age on feather δ ²H was investigated. We compared isotopic compositions between seasons using a Mann-Whitney test due to the non-normal distribution of the data. The result of the tests was interpreted based on the degree of evidence provided by the p-values (Muff et al., 2022). Furthermore, confidence intervals were calculated using the Monte Carlo resampling method (10,000 Iterations). Non-overlapping confidence intervals were considered as additional evidence of differences among samples.

To evaluate the effect of migration date, age and sex we employed linear regression models and the Information Theoretic framework with the Akaike index corrected for small sample sizes (AICc) (Dochtermann and Jenkins, 2011). We analyzed seven models for migration and four models for the winter season. That set of models included different combinations of age, sex, and migration date effects, along with a null model, to identify the most parsimonious models ( Table 1 ). To satisfy the assumption of data normality, we removed five extreme values from the dataset for the migration season. The ΔAICc value and the ranking relative to the null model served as evidence for the model(s) with the strongest support. In cases of model selection uncertainty, we used model averaged estimates, and we calculated 85% and 95% confidence intervals for parameter estimates. Predictive variables with 85 and 95% confidence intervals excluding zero, along with a relatively high AIC weight, were considered important (Arnold, 2010). In addition, we plotted the predicted values and 95% confidence intervals for variables in model set to visualize relationships. To assess the fit of models with stronger support, we examined residuals using graphical tools. All the analyses were conducted using the program R version 4.3.1 [R Core Team, 2023; R Project for Statistical Computing (RRID: SCR_001905)] and the Excel complement program Poptools V.3.2.3 [Hood, 2010; PopTools (RRID: SCR_022840)].

The transition probabilities from each breeding area to the region of Cuba were estimated using banding recovery data. These estimates were subsequently employed as prior information in the origin assignment analysis (Hobson et al., 2009). Banding data were requested from the North American Bird Banding Laboratory (BBL) (www.pwrc.usgs.gov/bbl). This database, which is freely accessible, includes capture and recapture information for thousands of bird species from 1955 to the present (Cohen et al., 2014).

The transition probabilities were estimated based on Burnham’s live-recapture dead-recovery modeling framework (Cohen et al., 2014). This model utilizes encounter histories of birds that were captured, marked, released, and subsequently either recaptured, resighted, recovered, or never re-encountered. The model provides estimates for four demographic parameters: survival probability (S), live recapture probability (p), dead recovery probability (r), and transition probability between strata (Psi). The encounter histories were constructed assuming a life span of 12 years (Rohwer et al., 2020), for individuals banded during the reproductive season (May to September) and with records during the migratory and wintering season (September to April).

For BWTE, the estimation of transition probabilities initially included four breeding areas in North America: Western Prairies (Alberta, Saskatchewan and Montana), Eastern Prairies (Manitoba, North and South Dakota), Upper Midwest (Wisconsin, Minnesota, Iowa) and Eastern Region (Ontario, Quebec, Michigan, Ohio, Pennsylvania, Maryland, Delaware, New York, New Jersey, Connecticut, Massachusetts, Rhode Island, New Hampshire, Vermont, Maine). These breeding areas are a modification of the proposals of Szymanski and Dubovsky (2013), based on the availability of banding information and their respective contributions to different harvest areas.

The delineation of populations and their distributions is a fundamental step in assessing migratory connectivity, making it important to establish breeding populations based on ecological criteria (Cohen et al., 2017). Subsequently, validation of these delineated areas was conducted utilizing population characteristics. We estimated population trajectories using the bbsBayes package (Edwards and Smith, 2021), which allows users to access breeding bird survey (BBS) data and to model annual indices and trends using hierarchical Bayesian models. We used a general additive model with a random effect of year (chains = 3, burn-in = 20,000, iterations = 10,000, thinning = 10, samples = 3,000). We used a t-distribution to model the extra-Poisson error distribution, as suggested by Edwards and Smith (2021). We modelled annual indices (1970–2018) stratified by the four regions: Western and Eastern Prairies, Upper Midwest, and Eastern. Based on the results, it was decided to consolidate the original number of breeding regions from four to three, merging the Western and Eastern Prairie regions due to their similarity. Details of these analyses can be found in Supplementary Material 1 .

The nine non-breeding regions corresponded to the flyways in which the recapture and recovery records were grouped in the original banding data (Canada, Pacific Migratory Flyways (MF), Central MF, Mississippi MF, Atlantic MF, and Central America) along with northern South America, the Caribbean, and Cuba (which was treated separately from the Caribbean). These nine regions satisfy the assumption that the recovery, resighting, and survival probabilities of banded individuals are uniform within the same area. This assumption is crucial for estimating transition probabilities (Cohen et al., 2014).

Reports of recapture and recovery often exhibit variations among non-breeding regions, with greater irregularities observed in the Neotropics. To reduce the error in estimating these parameters, covariates of effort for ‘r’ and ‘p’ were included in the model, as suggested by Cohen et al. (2014). Effort covariates for each region were calculated based on the number of recaptures and recoveries for all Neotropical migratory duck species in the BBL database relative to the area of available habitats. Available habitats were identified as wetlands and waterbodies, extracted from the Global Land Cover Map of 2015 (Kobayashi et al., 2017). Survival probability was assumed to be constant across all regions. For the analysis, we had data on 87,567 individuals banded in the breeding study areas, including 4,369 total records of recapture or recovery 897 of which were from Cuba. All analyses were conducted using the RMark package (Laake, 2013).

A Bayesian framework was used to establish and map the area of origin of individuals sampled (Vander Zanden et al., 2014; Campbell et al., 2020). Stable-isotope values within waterfowl flight feathers are representative of the natal site for hatch-year birds and the molting site for adults, as they are grown during a synchronous post-breeding molt.

We created an isotopic distribution map (isoscape) for North America, representing the expected δ ²H values for BWTE feathers (δ ²H_f). To accomplish this, a mean growing-season precipitation (δ ²H_p) isoscape (Bowen et al., 2005; Bowen, 2024) was calibrated using the equation based on mallard (Anas platyrhynchos) feathers (Equation 1) (van Dijk et al., 2014), as no calibration data of known-origin feathers were available for this species.

We estimated the likelihood that each cell (pixel) on the map represents a potential area of origin for each individual, using a normal probability density function. This function compares the individual feather δ ²H value against the expected (mean) δ ²H_f value at a given cell, while accounting for expected variation (i.e., calibration error and isoscape error). We used a residual standard deviation of 12.8‰ (Palumbo et al., 2019). The assignment of the area of origin was confined within the known breeding range (BirdLife International and Handbook of the Birds of the World, 2024) and the most relevant breeding area origin for Cuba based on the banding data analysis.

In order to encompass all possible origins of the sample and facilitate the interpretation of probable surface maps, a hierarchical clustering approach was applied to group individuals into similar geographic regions of origin for each season. Following Campbell et al. (2020), we generated a surface similarity matrix by conducting pairwise comparisons among the probability surfaces for all individuals using the Schoener’s D-metric. Subsequently, a hierarchical clustering analysis was performed on each similarity matrix using the average method. The approximately unbiased p-values (AU) of each cluster were computed via bootstrap resampling. To determine clusters, we cut the trees at a height of 0.5.

The estimated probabilities surfaces were normalized to sum to 1 and binarized by selecting cells with origin probability values exceeding 67%, corresponding to a 2:1 ratio (Van Wilgenburg and Hobson, 2011). Layers depicting the areas of probable origin for each individual were summed for each season and within each origin cluster for the season. This resulted in a representation of origin based on the number of individuals assigned to each pixel. All analyses were conducted using the isocat package version 0.2.6 (Campbell et al., 2020).

To comprehend the epidemiological implications of BWTE migratory connectivity in the context of Cuba, an analysis of the annual dynamics of Avian Influenza Virus (AIV) prevalence was conducted. This analysis encompassed both the molting and natal areas and the likely migratory flyways. To accomplish this, the predictive prevalence models by Kent et al. (2022) were used, which are freely accessible. Those authors provided a comprehensive spatiotemporal assessment of AIV prevalence in 30 duck species, utilizing data from two major surveillance databases in the continental United States: the U.S. Department of Agriculture (USDA) Animal and Plant Health Inspection Service (APHIS) National Wild Bird Avian Influenza Surveillance Program and the National Institutes of Health (NIH) National Institute of Allergy and Infectious Diseases (NIAID) Influenza Research Database.

For the analysis of the annual dynamics AIV prevalence, we worked with the average prevalence values for BWTE derived from the two aforementioned databases, covering all counties within each area of interest. With these data, a maximum and average weekly prevalence value was calculated for each area. These results are specific to the United States, as Kent et al. (2022) indicated that the data availability for Canada was insufficient to ensure robust predictions.

The most probable molting or natal area of individuals sampled during migration in western Cuba was extensive and located further north than those sampled in the central region of Cuba during winter. The substantial isotopic variation observed in the migratory season samples, along with strong evidence of differences from the wintering season samples, suggests a diverse range of potential origins for BWTE individuals using Cuba as a wintering or stopover site en route to the rest of the Caribbean, Central, and South America. However, during the migratory season, the probable area of origin for the samples corresponds to wetlands in the forests and prairies of the central-western United States and southern Canada, whereas during the winter season, it is confined mainly to prairie regions within the United States.

Both identified areas align with the regions characterized by the highest density of breeding pairs for the species, spanning from northeastern Iowa and the Dakotas to southern Manitoba, Saskatchewan, and Alberta (Rohwer et al., 2020). In this core breeding area, BWTE was one of the most abundant duck species, while in the remaining parts of its distribution range, pairs are found at lower densities and in irregular patterns.

During migration, BWTE individuals typically align their route with the migratory flyway closest to their breeding area; although there is approximately a 4.8% probability that some will cross to another flyway (Sharp, 1972; Garvon et al., 2011). The migration to wintering areas spans a duration of 6–8 weeks (Szymanski and Dubovsky, 2013). According to the analysis of Szymanski and Dubovsky (2013), a significant portion of individuals use the Central and Mississippi flyways with a smaller contingent migrating via the Pacific and Atlantic flyways. More specifically, birds originating from core breeding areas tend to follow the course of the Mississippi River Valley, ultimately reaching Mexico, Central America, and South America through Texas and Louisiana (Bellrose, 1980; Garvon et al., 2016).

Our results, which draw on banding data analysis, validate the notion that individuals predominantly utilize the migratory flyway closest to their breeding grounds. Specifically, for the Prairie breeding population, the primary migratory routes align with the Central and Mississippi flyways. The relative stability that these connectivity patterns support inferences of increased risk of AIV incursion in certain areas of Cuba based on publicly available information on virus detection in wild birds in the regions of origin at least in Canada and the United States (CWHC, 2024; USDA, 2024).

According to Kent et al. (2022) the prairie region of the United States, which coincides with the core breeding area of BWTE, is one of the regions with the highest incidence of AIV primarily in the months of October to November and mainly associated with the H5 subtype. These findings were exclusively derived from USDA data. However, we observed that, for birds from the Prairies, the peak starts earlier (August with peak in September). While Kent et al. (2022) focused on the U.S. region, it is likely that a similar pattern occurs in the core breeding areas of southern Canada, given their geographic continuity and shared ecoregion. Indeed, Papp et al. (2017) reported that the likelihood of a duck testing positive for AIV reached a peak in mid-August in the Prairie Provinces of Canada. On the other hand, for the U.S. segment of the Central and Mississippi Flyway, the proportion of positive cases remains relatively high until March during the spring migration. However, the elevated values observed during the spring season are likely associated with the H7 subtype (Kent et al., 2022).

We determined that BWTE sampled in Cuba during a period of likely stopover of birds en route to wintering sites further south (Central and northern South America) suggests the possibility that some of the BWTE individuals that use the island as a stopover location may ultimately reach one of these two regions. The highest probability of transition for these birds was attributed to the Upper Midwest and Prairie breeding population. This aligns with the analysis of band recoveries made by Fernández-Ordóñez and Albert (2022) for birds recovered in Venezuela for BWTE and suggests a fall leapfrog migration pattern for this species.

Distinctions in transition probabilities were identified between Cuba and the rest of the insular Caribbean, with the highest transition probabilities associated with the Upper Midwest and Eastern BWTE breeding areas. This suggests the presence of various migratory routes within the Caribbean. However, we are not aware of any similar connectivity pattern described for BWTE in the Caribbean and note a tendency for researchers to treat the Caribbean as a single region (Cohen et al., 2014). Possibly, the route passing through the smaller Caribbean islands connects with central and eastern regions of Cuba but remains undetected due to regional biases in band capture and reporting probabilities. This underscores the need to incorporate samples from other Cuban areas when investigating BWTE migration.

We found compelling evidence indicating that both migration date and individual age significantly influenced feather isotopic composition. Specifically, based on δ ²H_f, we observed that individuals tend to exhibit a more northerly origin as the migratory season progresses. Generally, there is a decreasing trend in the AIV prevalence with lower latitudes (Munster et al., 2007; Bevins et al., 2014). This gradient is associated with environmental factors such as temperature or host population density (Papp et al., 2017), as well as the temporal dynamics of virus circulation itself, with the highest prevalence typically occurring at the onset of the autumn migration (Olsen et al., 2006; Latorre-Margalef et al., 2009).

For BWTE, being among the first duck species to initiate migration, a decreased probability of infection by pathogens and parasites may result. Consequently, many individuals arrive in non-breeding regions without prior immunological experience (Garvon et al., 2016; Humphreys et al., 2020). This phenomenon would explain the species’ relatively high AIV prevalence later during the winter season in regions such as the Gulf Coast and Neotropics (Hanson et al., 2005; Ferro et al., 2010; González-Reiche et al., 2012; Hill et al., 2012). All of these findings underscore the critical importance of the initial phase of autumn migration for Cuba (Ferrer et al., 2014). During this phase, individuals tend to occupy molting areas to the south, increasing the likelihood of lacking prior immunological experience and, consequently, raising the possibility of acquiring the virus during the migration process.

There was a slight indication that adult BWTE individuals sampled during migration in Cuba occupied molting areas further south in North America than juveniles. Notably, higher AIV prevalence and viral shedding have been observed in juvenile (hatch-year) individuals within wild bird populations (Hill et al., 2012; Hoye et al., 2012; Avril et al., 2016). This outcome is linked to the fact that juveniles possess an immune system lacking immunity against AIV, in contrast to adults who may have acquired it in prior exposure. For the BWTE, adult males are among the first to migrate and can become abundant on the Gulf of Mexico coasts starting in late August (Bellrose, 1980; Thurber et al., 2020). However, as of the sampling initiation date in September, no discernible temporal variation pattern was observed in the age distribution of migrants.

Part of our study was carried out during the winter season in the central region of Cuba, specifically in one of the country’s most significant wetlands, the Southern Sancti Spíritus Wetland (BirdLife International, 2020). This wetland stands out not only for its high concentration of BWTE but also for its recognition as a crucial area in the circulation of other viruses that utilize birds as hosts, such as West Nile virus (Pupo-Antunez et al., 2018). Nevertheless, there are other wetlands in the southern region of Cuba that are used by this species during winter. For instance, the Southern Los Palacios Wetland in the western region is noted for its remarkable richness and abundance of aquatic species (Aguilar et al., 2020). Furthermore, due to its geographical position as one of the primary entry points for migration and its proximity to poultry commercial farms (Ferrer et al., 2014; Montano et al., 2020), this area is important in terms of risk in AIV transfer between wild and domestic birds. In the future, it would be advantageous for studies to encompass additional areas like these to enhance our comprehension of the migratory connectivity of species utilizing Cuba as their winter habitat.

Our establishment of migratory connectivity of BWTE using Cuba emphasizes opportunities for addressing various research issues associated with pathogen circulation (Altizer et al., 2011). Furthermore, our approach offers a pathway to scientifically reinforce risk-based surveillance solutions for early alert, prevention and control of these pathogens in Neotropical regions. The pattern of connectivity found in our study, together with knowledge of AIV infection rates at migration origin sites (Hobson et al., 2014), allows for a more accurate stratification of the risk of AIV incursion at destination and, consequently, better possibilities for risk-based planning towards resilience. The next phases of research should aim to integrate information on the infection dynamics of birds throughout their life cycle and their migratory behaviors, spanning both individual and population levels (Lebarbenchon et al., 2009; van Toor et al., 2018).

Through these combined analyses, we can identify critical areas along migratory flyways, including important stopover sites, which play a pivotal role in virus transmission and its potential movement across borders. Clearly, we have demonstrated the added value of risk assessment based on migratory connectivity deduced by stable isotope studies in BWTE.