Introduction

The Qinghai-Tibet Plateau (QTP) is the most extensively elevated surface on Earth, with an average elevation of ∼5 km over an area of 2.5 million km² (1). The uplift of the plateau exerts profound influences on the environment and determines the biogeographic boundaries both within and across continents (2). The unique geological developments of the QTP, especially its high elevation, are believed to have influenced various taxonomic groups continuously (3–5). However, the plateau’s uplift also brings up the Asian monsoon, one of the most vigorous phenomena in the global climate system (6, 7). The Asian monsoon dominates large areas extending from the Indian sub-continent eastwards to Southeast and East Asia (8). Its evolution and variability have caused significant variations in the redistribution of water and heat via a series of natural processes, such as drought, flood, and heat waves (6). Given the large impacts of Asian monsoons on the environment, they can reconfigure the spatial patterns of biodiversity and ecosystem processes. This entails movement patterns that shape the effects of the environment on organisms (9–11). However, owing to the difficulties in studying the complex effects caused by monsoons, most studies that explored the influence of the QTP just conceived the plateau as an orographic barrier (12–14). The role of monsoons in shaping species movement patterns remains poorly understood.

Animal movement sustains their spatial distributions and changes in multiple ways in response to environmental changes. A classical animal movement behaviour is migrating between breeding and wintering grounds to fulfil the annual cycle (15–17). Those migratory journeys that billions of animal individuals embark on land, in the oceans and skies, crossing national boundaries and continents, have fascinated humans for millennia (18). As a result, this has intrigued a body of different approaches and indicators to describe and model migration, including migratory direction, speed, timing, distance, and staging periods (18, 19). Amongst them, the migratory direction is one of the most prominent indicators for migration patterns, evidenced by a majority of animals migrating latitudinally between wintering and breeding areas (20). This can be explained by not only the fact that wintering sites are usually located in the warmer south (e.g. tropic) and breeding sites located in the cooler north (e.g. arctic), but also the earth’s magnetic fields that are arguably believed to affect the latitudinal migration of animals (21–24). However, the migratory direction can be changed from latitudinal to longitudinal when the animal faces environmental changes (20, 25–28).

Environmental fluctuations in the QTP are relatively small over the longue durée after the final-stage uplift (29), but few studies have evaluated how environmental heterogeneity across the QTP might influence the migratory behaviour of birds (but see migratory pattern descriptions, e.g., Zhao, Heim, Nussbaumer, van Toor, Zhang, Andersson, Bäckman, Liu, Song, Hellström, Roved, Liu, Bensch, Wertheim, Lei and Helm (30), Pu and Guo (31)). Yet it remains unclear whether and how these shifts systematically alter species’ migration patterns rather than a simple assumption that the QTP birds exploit resources according to their availability. Therefore, testing whether migration patterns vary consistently for birds that migrate across the QTP is key to our understanding of the processes that determine movement patterns and provides insights into how they may affect community organisation and functioning under the context of global environmental change.

In this work, we leverage the community-contributed and satellite-tracking data to explore the impacts of the QTP uplift in terms of both the development of its high elevation and Asian monsoons on the migratory strategies for the birds that migrated across the plateau. We do this by reconstructing the environments before the uplift and contrasting migratory directions of 50 bird species between breeding and wintering areas in environments before the uplift with those at present. We also calculate the migratory directions between adjacent stopover sites and assess the relationship between migratory directions and environmental stress. Our findings yield the most comprehensive picture to date of how the QTP uplift shapes migratory patterns of birds, revealing insights into the challenges and opportunities for migratory birds in a changing world.

Results and discussions

We find migratory birds tend to cluster both their wintering and breeding distribution with the uplift of QTP (Figure 1-E, G, and I). Specifically, before the uplift, migratory birds have a higher probability of breeding across a vast area at low and middle latitudes on the Eurasia continent, including West Asia, Siberia, QTP regions, and even Africa, whilst their most likely breeding areas move northeastward to the extreme north of Russia after the uplift. Different from the breeding area, the wintering area of migratory birds has a larger change in distribution probability. Birds that migrate across the QTP in the real scenario have a higher probability of wintering in Southwest Asia and North Africa, whereas they have a higher probability of moving southeast to winter in Southern China and more areas of Africa before the uplift (Figure 1-D F and H). Taken together, the distribution of migratory birds had extended in longitude and narrowed in latitude, suggesting that they are more likely to migrate along a longitudinal gradient in present environments as a result of the QTP uplift.

Figure 1.

Our results show that wind cost, temperature, and precipitation are three major factors that influence the overall migratory directions (both autumn and winter) of birds, despite the differences in autumn and spring migration across different geographic areas (Figure 2). During autumn migration, wind cost is the most important factor for birds’ migration direction (Figure 2). A higher wind cost is associated with migration, which suggests a higher opportunity for birds to use the wind to facilitate their longitudinal migration (Figures S1 and S9). They also choose to follow a flyway of relatively higher annual precipitation and temperature as they migrate from breeding areas to wintering areas during autumn migration (Figure S2, S3 and S9). Apart from those three factors, no evidence is found for strong impacts of elevation and vegetation on the direction of migration (Figure 2 and Figure S9).

Figure 2.

Birds diversify their preferences in terms of environmental conditions when they migrate in different geographic areas, i.e., areas west of (longitude < 73°E, West QTP), areas in the central (73°E ≤ longitude < 105°E, middle QTP), and areas east of the QTP (longitude ≥ 105°E, East QTP). Despite the fact that wind cost is the most important factor for the overall autumn migration, temperature is the most prominent factor in the areas east of the QTP and in the central QTP (Figure 2 and Figure S9). In addition, the average annual temperature in the central QTP is lower than that in the areas east of the QTP, but birds migrate across those two areas with increasing temperatures consistently (Figure S2, S6 and S9). Once they reach the regions west of the plateau (West QTP), low wind cost in the longitudinal direction and higher precipitation become priority choices for their migration (Figure 2 and Figure S9).

Compared with autumn migration, higher temperatures act as a major clue in the areas east and west of the plateau during spring migration, whereas the westerly outweighs temperature when birds migrate in the central plateau (Figure 1, Figure S1, S3 and S9). Besides temperature, precipitation also plays a role in all stages of spring migration. When birds migrate, they tend to follow a flyway of decreasing annual precipitation. Elevation has a slightly larger impact at the early stages of migration, i.e., areas east of the plateau during autumn migration and areas west of the plateau during spring migration, as birds migrate towards higher elevations during these stages.

It is commonly claimed that the initiation of migration is inherently inflexible in migratory birds (32), owing to the weak or insufficient responses by migratory birds to adjusting migration behaviour (e.g., migration timing and route) (33). This claim is particularly invoked for long-distance migrants, who may face greater temporal (e.g., migration timing) or physiological constraints given the varied phenologies en route (33). Our results show that a major change in avian migratory patterns in response to environmental change can be adjusting migration direction from the latitudinal to the longitudinal at the scale of their whole migration circle. This highlights substantial changes in migratory bird distribution and their biogeographic patterns as a result of the uplift of the Qinghai-Tibet Plateau (Figure 1).

One of the biggest climatic consequences of the uplift of the Qinghai-Tibet Plateau is the development of a unique monsoon system that has shaped environments across continents (34). One typical feature of Asian monsoons is the seasonal climatic change, comprising a dry cold winter phase and a wet warm summer phase (34). Asian monsoons also consist of several sub-systems, including the northeast monsoon and the East Asian winter monsoons that dominate the weather and climate in different parts of the plateau across different geographical periods. Our results showed that wind cost, temperature, and precipitation have more important impacts on avian migration than elevations in different geographic areas (Figure 2). This suggests that the monsoon system, rather than the high elevations of the plateau per se, is an important factor during avian migration on the plateau (Figure 2). Specifically, when birds begin their autumn migration in early September, the influence of the Siberian High on migration emerges as the East Asian winter monsoons start to reach the area east of the QTP, and their impacts at this stage are mainly reflected by varied temperatures and relatively less precipitation (35, 36). This can explain why higher temperatures and more precipitation play a more important role than wind cost in the area east of the plateau (Figure 2) since higher annual temperatures and more precipitations mean more food resources for migrants (37, 38), whereas wind during this period is less strong than that in winter (36). Whilst birds migrate westwards, less wind cost is becoming more important to determine their migration direction. This is, on the one hand, because the northeast monsoon begins to dominate the climate in the southwest of the QTP from around the end of September (39). The northeast monsoon brings cold wind to sweep the Qinghai-Tibet Plateau down towards the vast spans of the Indian Ocean (39), which could facilitate the westward migration of birds. On the other hand, in the northwest of the Qinghai-Tibet Plateau, the extended Siberia High and associated atmospheric systems that deliver cold and dry air masses to the Mediterranean surface also provide positive wind conditions for migrants (40).

When migrating toward breeding grounds in spring, birds adopt strategies different from their autumn migration. Temperature becomes more important than wind cost in spring migration (Figure 2 and S9). Given high temperatures usually mean relatively rich food resources for birds (41, 42), this suggests that birds that migrate across the Qinghai-Tibet Plateau may focus on energy accumulation during their spring migration rather than reducing flight costs in an effort to meet the energetic demands. Those birds also tend to follow a more ‘capital’ breeding strategy where birds rely on endogenous reserved energy gained prior to reproduction (43) rather than an ‘income’ breeding strategy where birds ingest nutrients mainly collected during the period of reproductive activity (44). This collaborates with studies on breeding strategies of migratory birds in Asian flyways (45). Another important reflection of the migration strategy is the role of elevation during different geographic areas of migration. As the elevation in the areas east and west of the plateau is much lower than that in the central, when migratory birds fly across the plateau, they need to follow a turned “u-shaped” elevation distribution and fly toward a higher elevation at the early and late stages of spring migration. This also implies that more energy is needed at the beginning of spring migration. Considering the need to balance energy accumulation and flight cost, areas with higher annual average temperatures and precipitation with high-level food resources are preferred during this migration stage to meet birds’ energy requests.

Our study provides a novel understanding of how QTP shapes migration patterns of birds. Albeit with the extensive influence of the plateau uplift on geology and geography, the resultant monsoon system, rather than its high elevation, is found to be a key factor shaping present avian migration patterns. Our study unveils shifts in avian migratory directions and their underlying mechanisms in the contexts of the QTP uplift, enhancing comprehension of the complex biogeographic effects on animal migration.

Methods and materials

Summary

We used two approaches to determine the migratory flyways of birds across the Qinghai-Tibet Plateau. First, we quantified the distributional change of each avian species by comparing the distribution range before and after the uplift of the plateau. For the present distribution, we used a dynamic spatiotemporal abundance model - Adaptive Spatiotemporal Model (AdaSTEM) that we have developed - to obtain the seasonal distribution of birds (19). We then used a species distribution model (i.e., MaxEnt) to measure the correlation between present distribution and environments (46). We calculated the distribution of migratory birds before the uplift of the plateau by projecting the correlation between their current distribution and environments onto environments before the uplift. Second, we obtained the specific migratory routes for each species by measuring the migratory directions (i.e., the azimuth angle between stopover sites) en route. Similarly, we used the relationships between present migratory directions and environments to predict the migratory directions pre-uplift of the plateau. Since our aim here was a prediction, we used random forest models, but we also used Bayesian multivariate regression modelling to measure the influence of environments on migratory directions of birds.

eBird checklist

We used a community-contributed database for the dynamic spatiotemporal abundance model to measure the seasonal distribution. Specifically, we first obtained the list of bird species that might migrate across the Qinghai-Tibet plateau based on Prins and Namgail (47). We then requested and downloaded the eBird Basic Dataset in Feb 2022 (48) for a sum of 64 species. We used data from the year 2019 to avoid the potential influence of the pandemic on bird observation (49, 50) and bird behaviour (51). For each bird observer, we calculated the historical cumulative species count for the bird observer throughout their historical eBird checklists prior to 2019 as a proxy to measure the expertise of bird observers. We then filtered the checklists as suggested by recent studies (52, 53):

1. Only checklists labelled as complete were included.

2. Only checklists with Traveling or Stationary protocol types were included. For checklists with the protocol type Traveling, only those with a travelling distance of less than 3km were included.

3. The observation duration should be longer than 5 min and shorter than 300 min.

4. Observers with expertise lower than 2.5% percentile were removed since they are less representative and may induce large bias.

To address spatiotemporal imbalances in data distribution and the potential overrepresentation of birding hotspots, we conducted spatiotemporal subsampling following the method proposed (52, 53). We first assigned each checklist with a global hexagonal hierarchical geospatial indexing system (H3 system; 54, 55), with a resolution of level 7 (∼5.16 km per cell). Then, to avoid biased sampling in rare species with unusual active temporal periods, we split the 24 hours of the day into 12 equal bins and assigned a checklist to each of the bins. We then randomly subsampled only one checklist for each year - day of the year - hour bin of day - cell combination. The subsampling resulted in 5,037,088 checklists for the year 2019.

Predictor variables for spatiotemporal abundance modelling

For each remaining checklist, we extracted six types of environmental variables based on their geographical coordinates:

1. Sampling effort variables, which include protocol type, travelling distance, observation duration minutes, number of observers, and observers’ expertise (measured in historical species count).

2. Temporal variables, which include day of year and observation started time of day.

3. Topographic variables, which include the mean and standard deviation of elevation, slope, north, and east aggregated within the 3 km × 3 km buffered area for each checklist. The Topographic data was downloaded from EarthEnv (56) in a 1 km resolution.

4. Land cover data. We used the Copernicus Climate Change Service (C3S) Global Land Cover data with a 300 m resolution (57). We calculated the landscape variables for each of the land cover types presented in the 3 km × 3 km buffered area for each checklist, including percentage cover, patch density, largest patch index, edge density, mean patch size, standard deviation of patch size for each land cover type, and entropy across heterogeneous land cover patches.

5. Bioclimate variables. We downloaded the ERA-5 hourly data at a 0.25° resolution (58). Hourly data of a 2 m temperature and total precipitation layer were firstly aggregated to daily level by taking the average. The day-level data were calculated using 19 bioclimate variables, which were then assigned to each checklist according to the geographical coordinates.

6. Normalized Difference Vegetation Index (NDVI). NDVI data were extracted from Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices 16-Day (MOD13A2) Version 6.1 product with a resolution of 16 days and 1 km (59). We further aggregated the data to hexagon level 5 based on the H3 indexing system (with edge length ∼9.85km). For each hexagon, we leveraged the pyGAM package (60) to apply a GAM model with 30 splines to interpolate the data to temporal ranges that were not provided by the original data. This resulted in a daily resolution dataset. We calculated six features based on NDVI and included them in subsequent modelling, i.e., the median, maximum, and minimum of NDVI, and the median, maximum, and minimum of the first derivative of NDVI against day of year (sometimes referred to “green wave”) for each hexagon throughout the year.

The feature engineering resulted in 106 predictor variables, including 6 sampling effort variables, 2 temporal variables, 8 topographic variables, 19 annual climatic variables, 65 land cover variables, and 6 vegetation index-related variables. All calculations are conducted in Python version 3.9.0.

Spatiotemporal abundance modelling

To adjust for sampling error and obtain the general migration pattern at the species level, we applied an Adaptive Spatio-Temporal Model (AdaSTEM) for each species to model abundance using stemflow package version 1.0.9.1, which we have recently reported (19).

AdaSTEM is a machine learning modelling framework that takes space, time, and sample size into consideration at different scales. It has been frequently used in modelling eBird data Fields (53, 61, 62) and has been proven to be efficient and advanced in multi-scale spatiotemporal data modelling. To briefly summarise the methodology, in the training procedure, the model recursively splits the input training data into smaller spatiotemporal grids (stixels) using the QuadTree algorithm (63). For each of the stixels, we trained a base model only using data contained by itself. Stixels were then aggregated and constituted an ensemble. In the prediction phase, stemflow queries stixels for the input data according to their spatial and temporal indexes, followed by the prediction of corresponding base models. Finally, we aggregated prediction results across ensembles to generate robust estimations (see Fink et al., 2013 (61) and stemflow documentation (19) for details).

We used XGBoost (64) as our classifier and regressor base model for its capability and balance between performance and computational efficiency. We set 10 ensemble folds, a maximum grid length threshold of 25 degrees, a minimum grid length threshold of 5 degrees, a temporal sliding window size of 50 DOY and a step of 20 DOY, and required at least 50 checklists for each stixel in model training. Trained models were then used to predict on prediction dataset with 0.1° spatial resolution and weekly temporal resolution, where the variables were annotated with the same methodology as that of the training dataset. Only spatiotemporal points with more than seven ensembles covered are predicted. In downstream analyses, we removed data points with abundance lower than 0.1 quantiles to obtain reliable predictions for each week.

Environmental variables for species distribution modelling

Given the challenges in simulating environmental and climatic conditions before the uplift of the Qinghai-Tibet Plateau, we modelled the environments before and after the uplift with five variables, i.e., monthly wind (speed and direction), annual temperature, annual precipitation, elevation and annual vegetation.

In detail, following Zhang, Jiang and Zhang (7), we used version 1.0.4 of the Community Earth System Model (CESM) coupled model with a dynamic atmosphere (CAM4), land (CLM4), ocean (POP2), and sea-ice (CICE4) components to simulate pre-uplift environments. CESM and its previous versions have been widely used in climate modelling, e.g., Meehl, Arblaster, Caron, Annamalai, Jochum, Chakraborty and Murtugudde (65) and are claimed to be capable of broadly reproducing the features of present-day climate (66). For CAM4, there is a horizontal resolution of ∼1.9° in latitude by 2.5° in longitude and 26 layers in the vertical direction. POP2 adopts a finer grid and has a nominal 1° horizontal resolution (320 × 384 grid points, latitude by longitude) and 60 layers in the vertical direction. The land and sea-ice components share the same horizontal grids as the atmosphere and ocean components, respectively. In CLM4, multiple land surface types and plant functional types (PFTs) are contained within one grid, and CLM4 can be run in a dynamic vegetation mode to simulate natural vegetation, including trees, grass, and shrub plant functional types, e.g., Yu, Wang, Parr and Ahmed (67) Qiu and Liu (68).

We initiated the modelling with two different scenarios, i.e., the actual elevation and a maximum elevation of 300m. We then used the same default preindustrial simulation for the two scenarios with a modern ice sheet, an atmospheric CO2 concentration of 280 ppmv, modern orbital parameters (the year 1950), modern solar constant (1,365 W/m²), other atmospheric greenhouse gas concentrations set to preindustrial values (CH4 and N2O set to 760 and 270 ppbv, respectively), and preindustrial aerosol conditions. We ran for 750 model years to ensure the combined atmospheric, ocean, and vegetation effects in response to the uplift of the plateau can be investigated.

Species distribution modelling

We used Maximum entropy (MaxEnt) models to compare the avian distributional change between pre- and post-uplift environments. MaxEnt compares the environmental features at presence points to those of pseudo absences to discriminate the suitable area (69). MaxEnt builds models using a generative approach and thus has an inherent advantage over a discriminative approach, especially when the amount of training data is small (69). Due to its good performance compared to other species distribution modelling techniques, MaxEnt is widely used in the study of biogeography and conservation biology.

We ran the MaxEnt model using default settings but with 1,000 iterations. For each model, we ran 20 bootstrap replications, and each time 75% of locations were selected at random as training samples, while the remaining 25% were used as validation samples.

Migratory direction

To obtain the species list of birds that migrate across the Qinghai-Tibet Plateau with available tracking data, we checked Movebank (movebank.org) together with literature reporting avian migratory routes across the plateau. For those who did not upload their data to Movebank, we digitalised the routes. Specifically, we built a new geographic layer with the same coordinate systems of each reported route and matched the layer with the images of routes. We then delineate migratory routes on the new geographic layer where the geographic information of the routes was achieved. This resulted in seven representative species that migrated across the plateau.

We used the same environmental variables, except for wind, for the species distribution model to investigate the potential influence of environments on migratory directions. We calculated wind connectivity to account for the influence of wind, considering wind connectivity has been identified as a key factor driving avian flying patterns (70). Since we aimed to investigate the migration patterns at large spatiotemporal scales, we measured the wind connectivity at a monthly resolution to enable analysis of seasonal differences. We adjusted the R package rWind for the computation. In detail, we replaced the default wind data from the Global Forecasting System atmospheric model with our monthly wind data from CESM as input. For both wind costs before and after the uplift of the plateau, we then calculate the movement cost from any starting cell to one of its eight neighbouring cells (Moore neighbourhood). This includes three parameters, i.e., wind speed at the starting cell, wind direction at the starting cell (azimuth), and the position of the target cell. A movement connectivity map was then determined after performing the default algorithms (71). We reversed the values of cells on the connectivity map, as we aimed to investigate the influence of wind cost, whereas the map showed the importance of the cell to maintain connectivity.

We used a random forest model and a multivariate linear regression model under the Bayesian framework to analyse the influence of environments on avian migratory directions. We first used the random forest model to measure the correlation between migratory directions and modern environments and predict the migratory direction before the uplift of the plateau. We then compared the influence between modern environments and environments before the uplift using a multivariate linear regression model under the Bayesian framework. We adopted two strategies for those two modelling approaches. First, we applied regression to different combinations of season-stage separately (seasons: spring, autumn; stages: overall, east QTP, central QTP, west QTP), resulting in eight regression models. Second, we additionally included species as random variables by applying hierarchical modelling, which also resulted in eight regression models.

All variables were standardised for comparison. All Bayesian models were conducted with PyMC version 5.5 (72) in Python version 3.9.0 environment. We used a NUTS sampler with a numpyro backend (jax.sample_numpyro_nuts) to run four chains, each with 30,00 tuning and 3,000 posterior chain sampling. We assessed the model convergence using potential scale reduction factor (Rhat) and effective sample size (ESS), where all parameters in all models met the criteria of Rhat < 1.03 and ESS > 400.

Data and materials availability

All data, code, and materials used in the analysis are available from GitHub (https://github.com/plmyann/QTPbirds).

Additional information

Funding

This work was supported by the CAS Project for Young Scientists in Basic Research (YSBR-097), the Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (No. 2019QZKK0501), and the National Natural Science Foundation of China (No. 31821001);

Author contributions

XZ, WZ, and ZG conceived the idea and designed the methodology of the study, WZ, ZG, YC and ZR performed the analyses and XZ shaped the study. WZ and XZ wrote the paper. All authors contributed critically to the drafts and gave final approval for publication;