Ecology

Fossil record of insectsWe compiled all species-level fossil occurrences of insects using https://paleobiodb.org/ (PBDB) as a starting point (downloaded October 12, 2021). The dataset obtained from PBDB contained initially 5808 occurrences for a period ranging from the Asselian to the Rhaetian. The dataset was cleaned of synonyms, outdated combinations, nomina dubia, and other erroneous and doubtful records, based on revisions provided in the literature and/or on the expertise of the authors. After correction, including data addition from the literature, our dataset was composed of 3636 species (1784 genera, and 418 families) for 17,250 occurrences resulting from an in-depth study and curation of the entire bibliography of fossil insects, spanning from the Asselian (lowermost Permian) to the Rhaetian (uppermost Triassic). Although most of the taxa included in the datasets are nominal taxa (published and named), a few unnamed taxa (genera or species) that are considered separate from others were also included, although not formally named in the literature or not published yet. These unpublished taxa are identifiable by the notation ‘fam. nov.’ or ‘gen. nov.’ following their names.Occurrences used here are specimens originating from a given stratigraphic horizon assigned to a given taxon. The age of each occurrence is based on data from PBDB, corrected with a more precise age (generally stage, sometimes substage), and the age of each time bin boundaries relies on the stratigraphic framework proposed in the International Chronostratigraphic Chart (updated to correspond with the ICS 2022/0295). Similarly, the ages of some species assigned to the wrong stage were corrected. In fact, some species from the French Permian deposit of Lodève were initially considered to be of Artinskian age in PBDB but most species from this deposit originate from the Merifons member, which is of Kungurian age96.Our data compilation allows a robust integration of data before and after our period of interest (i.e. the lower Permian and all geologic stages after the Carnian) to encompass occurrences of genera that may survive until the Late Triassic and to generate a sufficient background for the model to correctly estimate the extinction events around the P/T boundary. Since we used different datasets, the differences between genus-level or family-level occurrence numbers are explained by the systematic placement of some specimens that can only be placed confidently in a family but not in a genus (Supplementary Table 1). Tentative species identifications originally placed with uncertainty (reported as ‘aff.’ or ‘?’) were always included at a higher taxonomic level. Uncertain generic attributions were integrated as occurrences at the family level (e.g. a fossil initially considered Tupus? is recorded as an occurrence of Meganeuridae). Our total dataset was subdivided into smaller datasets, which represent orders or other subclades of insects (e.g. Mecoptera, Holometabola and Polyneoptera). Note that all the ichnospecies—a species name assigned to trace fossils (e.g. resting trace, nest and leaf damage)— and insect eggs (e.g. Clavapartus latus, Furcapartus exilis and Monilipartus tenuis) were not included in the analyses97. To prevent potential issues regarding the diversification estimates for clades with poor delineation, we refrained from analysing several orders that serve as taxonomic ‘wastebaskets’ (e.g. Grylloblattodea). These groups are poorly defined, likely polyphyletic or paraphyletic, and not supported by apomorphic characters—e.g. the monophyly of the ‘Grylloblattodea’ (Grylloblattida Walker, 1914 plus numerous fossil families and genera of uncertain affinities) is not supported by any synapomorphy, nor the relationships within this group. The occurrences assigned to these orders were rather included in analyses conducted at a higher taxonomic level (at the Polyneoptera level in the case of the ‘Grylloblattodea’). The detail of the composition of all the datasets is given in Supplementary Table 14, and each dataset is available in Supplementary Data 1.Studying extinction should, when possible, rely on species-level diversity to better circumscribe extinction events at this taxonomic rank, which is primarily affected by extinction98,99,100. However, in palaeoentomology, species-level occurrence data may contain less information than genus-level data, mainly because species are most of the time only known from one deposit, resulting in reduced life span, and are also sometimes poorly defined. Insects are also less prone to long-lasting genera or species than other lineages, maybe because of the relatively short time between generations (allowing for rapid evolution) or because morphological characters are better preserved or more diagnostic than in other lineages (i.e. wing venation), allowing easier differentiation. Another argument for the use of genus-level datasets is the possibility to add occurrences represented by fossils that cannot be assigned at the species level because of poor preservation or an insufficient number of specimens/available characters. By extension, the genus life span provides clues as to survivor taxa and times of origination during periods of post-extinction or recovery. A genus encompassing extinction events indicates that at least one species of this genus crossed the extinction. To get the best signal and infer a robust pattern of insect dynamics around the P/T events, we have chosen to analyse our dataset at different taxonomic ranks (e.g. genus, family and order levels) to extract as much evidence as possible.To further support our choice to work at these different levels, most recent works aiming to decipher the diversification and extinction in insect lineages have worked using a combination of analyses21,22,26; this also applies to non-insect clades51,101,102. This multi-level approach should maximise our understanding of the Permo–Triassic events.Assessing optimal parameters and preliminary testsPrior to choosing the settings for the final analyses (see detail in Dynamics of origination and extinction), a series of tests were carried out to better evaluate the convergence of our analyses. First, we analysed our genus-level dataset with PyRate36 running for 10 million generations and sampling every 10,000 generations, on ten randomly replicated datasets using the reversible-jump Markov Chain Monte Carlo (RJMCMC) model37 and the parameters of PyRate set by default. As the convergence was too low, new settings were used, notably increasing the number of generations to 50 million generations and monitoring the MCMC mixing and effective sample size (ESS) each 10 million generations. We modified the minimal interval between two shifts (-min_dt option, testing 0.5, 1.5 and 2), and found no major difference in diversification patterns between our tests. We have opted for 50 million generations with a predefined time frame set for bins corresponding to the Permian and Triassic stages, and a minimum interval between two shifts of two Ma. These parameters allow for maintaining a short bin frame and high convergence values while correctly identifying periods of diversification and extinction. For each analysis, ten datasets were generated using the extract.ages function to randomly resample the age of fossil occurrences within their respective temporal ranges (i.e. resampled ages are randomly drawn between the minimum and the maximum ages of the geological stratum). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 after excluding the first 10% of the samples as a burn-in period. The parameters are considered convergent when their ESS are greater than 200.Dynamics of origination and extinctionWe carried out the analyses of the fossil datasets based on the Bayesian framework implemented in the programme PyRate36. We analysed the fossil datasets under two models: the birth–death model with constrained shifts (BDCS38) and the RJMCMC (-A 4 option37). These models allow for a simultaneous estimate for each taxon: (1) the parameters of the preservation process (Supplementary Fig. 17), (2) the times of origination (Ts) and extinction (Te) of each taxon, (3) the origination and extinction rates and their variation through time for each stage and (4) the number and magnitude of shifts in origination and extinction rates.All analyses were set with the best-fit preservation process after comparing (-PPmodeltest option) the homogeneous Poisson process (-mHPP option), the non-homogeneous Poisson process (default option), and the time-variable Poisson process (-qShift option). The preservation process infers the individual origination and extinction times of each taxon based on all fossil occurrences and on an estimated preservation rate, denoted q, expressed as expected occurrences per taxon per Ma. The time-variable Poisson process assumes that preservation rates are constant within a predefined time frame but may vary over time (here, set for bins corresponding to stages). This model is thus appropriate when rates over time are heterogeneous.We ran PyRate for 50 million MCMC generations and a sampling every 50,000 generations for the BDCS and RJMCMC models with time bins corresponding to Permian and Triassic stages (-fixShift option). All analyses were set with a time-variable Poisson process (-qShift option) of preservation and accounted for varying preservation rates across taxa using the Gamma model (-mG option), that is, with gamma-distributed rate heterogeneity with four rate categories36. As explained above, the minimal interval between two shifts (-min_dt option) was modified and a value of 2 was used. The default prior to the vector of preservation rates is a single gamma distribution with shape = 1.5 and rate = 1.5. We reduced the subjectivity of this parameter, and favoured a better adequation to the data, allowing PyRate to estimate the rate parameter of the prior from the data by setting the rate parameter to 0 (-pP option). Therefore, PyRate assigns a vague exponential hyper-prior to the rate and samples the rate along with all other model parameters. Similarly, because our dataset does not encompass the entire fossil record of insects, we assumed that a possible edge effect may interfere with our analyses, with a strong diversification during the lowermost Permian and, conversely a strong extinction during the uppermost Triassic. Because the RJMCMC and BDCS algorithms look for rate shifts, we constrained the algorithm to only search for shifts (-edgeShift option) within the following time range 295.0 to 204.5 Ma. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 after excluding the first 10% of the samples as a burn-in period. The parameters are considered convergent when their ESS are greater than 200.We then combined the posterior estimates of the origination and extinction rates across all replicates to generate rates through-time plots (origination, extinction, and net diversification). Shifts of diversification were considered significant when log Bayes factors were >6 in the RJMCMC model, while we considered shifts to be significant in the BDCS model when mean rates in a time bin did not overlap with the 95% credibility interval (CI) of the rates of adjacent time bins.We replicated all the analyses on ten randomly generated datasets of each clade and calculated estimates of the Ts and the Te as the average of the posterior samples from each replicate. Thus, we obtained ten posterior estimates of the Ts and Te for all taxa and we used these values to estimate the past diversity dynamics by calculating the number of living taxa at each time point. For all the subsequent analyses, we used the estimated Ts and Te of all taxa to test whether or not the origination and the extinction rate dynamics were correlated with particular abiotic factors, as suggested by the drastic changes in environmental conditions known during the Permo–Triassic. We used proxies for abiotic factors, such as global continental fragmentation or the dynamic of major clades of plants, and for biotic factors via species interaction within and between ecological guilds. This approach avoids re-modelling preservation and re-estimating times of origination and extinction, which reduces drastically the computational burden, while still allowing to account for the preservation process and the uncertainties associated with fossil ages. Similarly, the times of origination and extinction used in all the subsequent analyses were obtained while accounting for the heterogeneity of preservation, origination and extinction rates. To discuss the magnitude of the periods of extinction and diversification, we compared the magnitude of these events to the background origination and extinction rates (i.e. not during extinction or diversification peaks).The PyRate approach has proven to be robust following a series of tests and simulations that reflect commonly observed biases when modelling past diversity dynamics31,38. These simulations were based on datasets simulated under a range of potential biases (i.e. violations of the sampling assumptions, variable preservation rates, and incomplete taxon sampling) and reflecting the limitations of the fossil record. Simulation results showed that PyRate is able to correctly estimate the dynamics of origination and extinction rates, including sudden rate changes and mass extinction, even if the preservation levels are low (down to 1–3 fossil occurrences per species on average), the taxon sampling is partial (up to 80% missing) or if the datasets have a high proportion of singletons (exceeding 30% of the taxa in some cases). The strongest bias in birth–death rate estimates is caused by incomplete data (i.e. missing lineages) altering the distribution of taxa; a pervasive effect often mentioned for phylogeny-based models104,105,106. However, in the case of PyRate, the simulations confirm the absence of consistent biases due to an incomplete fossil record36. Finally, the recently implemented RJMCMC model was shown to be very accurate for estimating origination and extinction rates (i.e. more accurate than the BDCS model, the boundary-crossing and three-time methods) and is able to recover sudden extinction events regardless of the biases in the fossil dataset37.The severity of extinctions and survivorsFor each event—the Roadian–Wordian, the LPME, and the Ladinian–Carnian—we quantified the percentage of extinctions and survivors at the genus level. We used the Te and Ts from our RJMCMC analysis and computed the mean for the Te (Tem) and for the Ts (Tsm) of each genus. We then filtered our dataset to keep only the genera with a Tsm older than the upper boundary of the focal event, i.e., we only kept the genera that appeared before the end of the event. Then, we discarded the genera that have disappeared before the lower boundary of the focal event, i.e. Tem older that the lower boundary of the event. The remaining genera, which corresponds to all the genera (total) present during the crisis (Ttgen), can be classified into two categories, ‘survivor genera’ (Sgen), i.e. those that survived the crisis, and those that died: ‘extinct genera’ (Egen). The survivors have a Tem younger than the upper boundary of the focal event, while the ‘extinct genera’ died out during the event and have a Tem between the lower and upper boundaries of the event of interest. To obtain the percentage of survivors, we used the following formula: (Sgen/Ttgen) × 100. Similarly, the percentage of extinction is calculated as: (Egen/Ttgen) × 100.Age-dependent extinction modelWe assessed the effect of taxon age on the extinction probability by fitting the age-dependent extinction (ADE; -ADE 1 option) model50. This model estimates the probability for a lineage to become extinct as a function of its age, also named longevity, which is the elapsed time since its origination. It is recommended to run the ADE model over time windows with roughly constant origination and extinction rates, as convergence is difficult—but not impossible—to reach in extinction or diversification contexts50. We ran PyRate for 50 million MCMC generations with a sampling every 50,000 generations, with a time-variable Poisson process of preservation (-qShift option), while accounting for varying preservation rates across taxa using the Gamma model (-mG option). We replicated the analyses on ten randomised datasets and combined the posterior estimates across all replicates. We estimated the shape (Φ) and scale (Ψ) parameters of the Weibull distribution, and the taxon longevity in a million years. According to ref. 50, there is no evidence of age-dependent extinction rates if Φ = 1. However, the extinction rate is higher for young species and decreases with species age if Φ  1. Although ADE models are prone to high error rates when origination and extinction rates increase or decrease through time, simulations with PyRate have shown that fossil-based inferences are robust50. We investigated the effect of ADE during three different periods (-filter option) as follows: (1) between 264.28 Ma and 255 Ma (pre-decline), (2) between 254.5 Ma and 251.5 Ma (decline) and (3) between 234 Ma and 212 Ma (post-crisis). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Selection of abiotic and biotic variablesTo test correlations of insect diversification with environmental changes, we examined the link between a series of environmental variables and origination/extinction rates over a period encompassing the GEE, the LPME and the CPE but also for each extinction event. We focused on the role of nine variables, also called proxies, which have been demonstrated or assumed to be linked to extinctions and changes in insect diversity26,67.The variations in the atmospheric CO2 and O2 concentrations are thought to be correlated with the diversification of several insect lineages, including the charismatic giant Meganeuridae65,66,67. Because the increase of O2 concentration has likely driven the diversification of some insects, its diminution may have resulted in the extinction or decline of some lineages. Therefore, we investigated the potential correlation of the variations of this variable with insect dynamics using data from ref. 55. We extracted the data, with 1-million-year time intervals, spanning the Permo–Triassic.Similarly, the modification of CO2 concentration, notably its increase, is known to promote speciation in some modern insect groups107. Therefore, a similar effect may have occurred during the Permian and Triassic but remains to be tested. We based our analyses on the dataset of ref. 108. We used their cleaned dataset and extracted all verified values for the Permo–Triassic interval. Because the initial data (i.e. independent estimates) were made in various locations for the same age, different values of the CO2 concentration are provided. We incorporated all these values in our analysis, allowing PyRate to search for a correlation for each value of the CO2 concentration. We obtain a final correlation independent of the sampling location, in line with our large-scale analysis.The continental fragmentation, as approximated by plate tectonic change over time, has recently been proposed as a driver of Plecoptera dynamics26. Because the period studied encompasses a major geological event, the fragmentation of the supercontinent Pangea, we investigated the effect of continental fragmentation on insect diversification dynamics. We retrieved the index of continental fragmentation developed by ref. 69 using paleogeographic reconstructions for 1-million-year time intervals. This index approaches 1 when all plates are disjoined (complete plate fragmentation) and approaches 0 when the continental aggregation is maximal.Climate change (variations in warming and cooling periods) is a probable driver of diversification changes over the history of insects21,109. Temperature is likely directly linked with insect dynamics109 but also with their food sources, notably plants110. Because it was demonstrated that modification of temperature impacted floral assemblages110, we tested the correlation between temperature variations and the diversification dynamic of insects. Major trends in global climate change through time are typically estimated from relative proportions of different oxygen isotopes (δ18O) in samples of benthic foraminiferan shells111. We used the data from ref. 112, converted to absolute temperatures following the methodology described in Condamine et al.113 (see their section Global temperature variations through time). The resulting temperature data reflects planetary-scale climatic trends, with time intervals inferior to 1-million-year, which can be expected to have led to temporally coordinated diversification changes in several clades rather than local or seasonal fluctuations.The fluctuation in relative diversity of gymnosperms, non-Polypodiales ferns, Polypodiales ferns, spore-plants, and later the rise of angiosperms has likely driven the diversification of numerous insects57,60,61,114. Close interactions between insects and plants are well-recorded during the Permian and Triassic57,60,61. In fact, herbivorous insects are known to experience high selection pressure from bottom-up forces, resulting from interactions with their hosts or feeding plants30,72. Therefore, it appears crucial to investigate the effect of these modifications on the insects’ past dynamics. We used the data from ref. 38 for the different plant lineages (all with 1-million-year time intervals). All the datasets for these variables are available in the publications cited aside from each variable or in Supplementary Data 1.Multivariate birth–death modelWe used the multivariate birth–death (MBD) model to assess to what extent biotic and abiotic factors can explain temporal variation in origination and extinction rates55. The model is described in ref. 55, where origination and extinction rates can change through time in relation to environmental variables so that origination and extinction rates depend on the temporal variations of each factor. The strength and sign (positive or negative) of the correlations are jointly estimated for each variable. The sign of the correlation parameters indicates the sign of the resulting correlation. When their value is estimated around zero, no correlation is estimated. An MCMC algorithm combined with a horseshoe prior, controlling for over-parameterisation and for the potential effects of multiple testing, jointly estimates the baseline origination (λ0) and extinction (µ0) rates and all correlation parameters (Gλ and Gµ)55. The horseshoe prior is used to discriminate which correlation parameters should be treated as noise (shrunk around 0) and which represent a true signal (i.e. significantly different from 0). In the MBD model, a correlation parameter is estimated to quantify independently the role of each variable on the origination and the extinction.We ran the MBD model using 20 (for short intervals) or 50 million MCMC generations and sampling every 20,000 or 50,000 to approximate the posterior distribution of all parameters (λ0, µ0, nine Gλ, nine Gµ and the shrinkage weights of each correlation parameter, ωG). The MBD analyses used the Ts and the Te derived from our previous analyses under the RJMCMC model. The results of the MBD analyses were summarised by calculating the posterior mean and 95% CI of all correlation parameters and the mean of the respective shrinkage weights (across ten replicates), as well as the mean and 95% CI of the baseline origination and extinction rates. We carried out six analyses, over: (1) the Permo–Triassic (between 298.9 and 201.3 Ma); (2) the Roadian–Wordian (R/W) boundary (between 270 and 265 Ma), (3) the LPME (between 254.5 and 250 Ma), (4) the Ladinian–Carnian (L/C) boundary (between 240 to 234 Ma), (5) the Permian period (between 298.9 and 251.902 Ma) and (6) the Triassic period (between 251.902 and 201.3 Ma). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Multiple clade diversity-dependence modelTo assess the potential effect of diversity-dependence on the diversity dynamics of three or four insect guilds, we used the multiple clade diversity-dependence (MCDD) model in which origination and extinction rates are correlated with the diversity trajectory of other clades31. This model postulates that competitive interactions linked with an increase in diversity results in decreasing origination rates and/or increasing extinction rates. The MCDD model allows for testing diversity-dependence between genera of a given clade or between genera of distinct clades sharing a similar ecology.We estimated the past diversity dynamics for three (i.e. herbivores, predators, and a guild composed of generalists + detritivores/fungivores dubbed ‘others’) or four insect groups or guilds (i.e. herbivores, predators, generalists and detritivores/fungivores) by calculating the number of living species at every point in time based on the times of origination (Ts) and extinction (Te) estimated under the RJMCMC model (see above) (Supplementary Figs. 19–24). We defined our four insect groups with a cautious approach i.e. insect genera, families or orders for which nothing is known about the ecology or about the ecology of their close relatives were not considered for the analysis. For example, no diet was assigned to Diptera, Mecoptera or Glosselytrodea. The ecology of the Triassic Diptera and Permo–Triassic Mecoptera is difficult to establish because extant Diptera and Mecoptera have a wide diversity of ecology. Fossil Mecoptera are also putatively involved in numerous interactions with plants (species with elongated mouthparts), suggesting a placement in the herbivore group, while other species were likely predators. Therefore, we cannot decide to which group each species belongs. Similarly, nothing is known about the body and mouthparts of the Glosselytrodea, most of the time described based on isolated wings; we did not assign the order to any group. The definition and delineation of insect clades have also challenged the placement of several orders (e.g. ‘Grylloblattodea’) in one of our four groups. The order ‘Grylloblattodea’ is poorly delineated and mostly serves as a taxonomic ‘wastebasket’ to which it is impossible to assign a particular ecology. Finally, genera, species, or families not placed in a higher clade (e.g. Meshemipteron, Perielytridae) were not included in the analysis. Oppositely, the guilds ‘herbivores’ and ‘predators’ are well defined, and their ecology is evidenced by the morphology of their representatives and the principle of actualism. For example, the ecology of Meganeurites gracilipes (Meganeuridae) has been deeply studied, and its enlarged compound eyes, its sturdy mandibles with acute teeth, its tarsi and tibiae bearing strong spines, and the presence of a pronounced thoracic skewness are specialisations today found in dragonflies that capture their prey while in flight115. All Odonatoptera are well-known predator insects. The raptorial forelegs of the representatives of the order Titanoptera and their mouthparts with strong mandibles are linked with predatory habits81. The Palaeodictyopteroidea were herbivorous insects with long, beak-like, piercing mouthparts, and probably a sucking organ81,82. Most Hemiptera are confidently considered herbivorous insects by comparison with their extant representatives. For example, the Cicadomorpha or Sternorrhyncha are known to feed on plants and their fossil representatives likely possessed the same ecology because of similar morphologies116. Some hemipteran families (e.g. Nabidae) are predators and we cautiously distinguished herbivorous and carnivorous taxa among Hemiptera. The detail of the ecological assignations for the 1009 genera included in our analyses can be found in Supplementary Data 1 (Table MCCD).We calculated ten diversity trajectories from the ten replicated analyses under the RJMCMC model. The estimation of past species diversity might be biased by low preservation rates or taxonomic uncertainties. However, such trajectory curves are likely to provide a reasonably accurate representation of the past diversity changes in the studied clades, notably because the preservation during the Permian and Triassic period is relatively good for insects (i.e. no gaps).Our MCDD analyses comprise all the insect genera spanning from the lowermost Permian to the uppermost Triassic and were run and repeated on ten replicates (using the Te and Ts estimated under the RJMCMC model) with 50 million MCMC generations and a sampling frequency of 50,000. For each of the four insect groups, we computed the median and the 95% CI of the baseline origination and extinction rates (λi and µi), the within-group diversity-dependence parameters gλi and gµi, and the between-groups diversity dependence parameters gλij and gµij. The mean of the sampled diversity dependence parameters (e.g. gλij) was used as a measure of the intensity of the negative (if positive) or positive interactions (if negative) between each pair of groups. The interactions were considered significant when their median was different from 0 and the 95% CI did not overlap with 0. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.We cross-validated the result of the MCDD model using the MBD model. The MBD model can be used to run a multiple clade diversity-dependence analysis by providing the diversity trajectories of insect guilds as a continuous variable. These data are directly generated by PyRate using the lineages-through-time generated by the RJMCMC analyses (-ltt option). We ran the MBD model using 50 million MCMC generations and sampling every 50,000 to approximate the posterior distribution of all parameters (λ0, µ0, four Gλ, four Gµ and the shrinkage weights of each correlation parameter, ωG). We carried out three analyses, over the period encompassing the three extinction events (between 275 and 230 Ma): (1) for herbivores; (2) for predators; and (3) for ‘others’. For each analysis, the lineages-through-time data of the two other guilds are used as continuous variables to investigate a diversity dependence effect. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

SamplesWe examined 1,022 avian species (~ 10% recorded species) in this study, representing 418 genera, from 91 families (37% recorded families) and 29 orders (73% of all orders). Specimens were from the skin collection of the World Museum Liverpool, Tring Natural History Museum, Manchester Museum and Wollaton Hall Museum, all situated in the United Kingdom. All work was carried out in accordance with ethical regulations at Manchester Metropolitan University and with the permission of all aforementioned museums. Only the best-preserved adult specimens (no signs of cut off feathers or holes in the skin near the beak) were chosen for this study to ensure accurate measurements of bristle length, shape and presence, which should not be affected by the process of skin removal and specimen conservation. Species were randomly chosen, without targeting our sampling towards species known a priori to have bristles. Where possible, two specimens per species were measured (occurring in 82% of all species examined). Specimens of each sex were measured when present; however, this was not always possible since labelling was often inaccurate or missing. In total, the sample included 508 males, 412 females and 374 individuals of unknown sex. Both sexes were examined in 274 species and there was no difference whatsoever between the presence of bristles on male or female species (n = 97 with bristles present and n = 180 with bristles absent for both males and females). Length (Mann–Whitney U test, W = 37,962, N = 552, P = 0.94) and shape (Chi-square test, χ2 = 0, N = 552, df = 3, P = 1) of rictal bristles also did not significantly differ between males and females. Therefore, rictal bristles are likely to be sexually monomorphic and data for males and females was pooled for further analyses. Overall, rictal bristles were absent in 64% of species examined (n = 656) and just over a third of species (n = 366) had bristles present.Bristle descriptionsFacial bristles were initially identified by sight and touch in each specimen. Bristles were recorded as either present or absent from the upper rictal, lorial, lower rictal, narial and interramal regions (Fig. 1a). We use the term ‘rictal bristle’ here for bristles on both the upper rictal and/or the lorial region, since there was no clear differentiation and morphological differences between the bristles found in these regions forming a continuum of bristles above the edge of the beak. When present, rictal bristle shape was recorded as: (i) unbranched rictal bristles, (ii) rictal bristles with barbs only at the base (“Base”) and (iii) branched rictal bristles (“Branched”), i.e. barbs and barbules present along the bristle rachis (Fig. 1b). The three longest rictal bristles were measured on both sides of the head of each specimen using digital callipers, and these lengths were averaged to provide a mean length of rictal bristles per species. In species lacking rictal bristles, a length of “0” and a shape category of “Absent” was recorded.Ancestral reconstruction of facial bristle presenceFollowing Felice et al.19, a single consensus phylogenetic tree was generated from the Hackett posterior distribution of trees from Birdtree.org20 with a sample size of 10,000 post burn-in, using the TreeAnnotator utility in BEAST software21 with a burn-in of 0. Maximum Clade Credibility (MCC) with the option “-heights ca” was selected as the method of reconstruction. The common ancestor trees option (-heights ca) builds a consensus tree by summarising clade ages across all posterior trees. Both the consensus tree and posterior distribution of 10,000 trees were imported into RStudio v. 1.2.5 for R22,23 and pruned so that only species present in the dataset of this study remained in the phylogeny. Taxon names were modified where necessary to match those from the Birdtree.org (http://birdtree.org) species record. Negative terminal branches in our consensus tree were slightly lengthened to be positive using ‘edge.length[tree$edge.length  More

Khan, M. et al. Plant species and communities assessment in interaction with edaphic and topographic factors; an ecological study of the mount Eelum District Swat Pakistan. Saudi J. Biol. Sci. 24(4), 778–786 (2017).Article 

Google Scholar 
Ur Rahman, A. et al. Impact of multiple environmental factors on species abundance in various forest layers using an integrative modeling approach. Global Ecol. Conserv. 29, e01712 (2021).Article 

Google Scholar 
Arneth, A., Uncertain future for vegetation cover. Nature 524(7563), 44–45.Goldsmith, F., Description and analysis of vegetation. Methods Plant Ecol. (1976).Rahman, I. U. et al. First insights into the floristic diversity, biological spectra and phenology of Manoor Valley Pakistan. Pak. J. Bot 50(3), 1113–1124 (2018).
Google Scholar 
Khan, S.M., Plant communities and vegetation ecosystem services in the Naran Valley, Western Himalaya, 2012, University of Leicester.Haq, F., Ahmad, H. & Iqbal, Z. Vegetation description and phytoclimatic gradients of subtropical forests of Nandiar Khuwar catchment District Battagram. Pak. J. Bot 47(4), 1399–1405 (2015).
Google Scholar 
Iqbal, M. et al. A novel approach to phytosociological classification of weeds flora of an agro-ecological system through Cluster, two way cluster and indicator species analyses. Ecol. Ind. 84, 590–606 (2018).Article 

Google Scholar 
Shaw, M. R. et al. Grassland responses to global environmental changes suppressed by elevated CO2. Science 298(5600), 1987–1990 (2002).Article 

Google Scholar 
Drenovsky, R.E., Effects of mineral nutrient deficiencies on plant performance in the desert shrubs Chrysothamnus nauseosus ssp. consimilis and Sarcobatus vermiculatus2002: University of California, Davis.Iqbal, M. et al. Vegetation classification of the Margalla Foothills, Islamabad under the influence of edaphic factors and anthropogenic activities using modern ecological tools. Pak. J. Bot 53(5), 1831–1843 (2021).Article 

Google Scholar 
Bai, Y. et al. Landscape-level dynamics of grassland-forest transitions in British Columbia. J. Range Manag. 57(1), 66–75 (2004).Article 

Google Scholar 
Zhao, T. et al. Retrievals of soil moisture and vegetation optical depth using a multi-channel collaborative algorithm. Remote Sens. Environ. 257, 112321 (2021).Article 

Google Scholar 
Austin, M., Chapter 2: Vegetation and environment: discontinuities and continuities. IN VAN DER MAAREL, E.(Ed.) Végétation ecology. Etats‐Unis, 2005, Blackwell Publishing.Peña-Claros, M. et al. Soil effects on forest structure and diversity in a moist and a dry tropical forest. Biotropica 44(3), 276–283 (2012).Article 

Google Scholar 
Miao, R. et al. Effects of long-term grazing exclusion on plant and soil properties vary with position in dune systems in the Horqin Sandy Land. CATENA 209, 105860 (2022).Article 

Google Scholar 
Abbas, Z. et al. Plant communities and anthropo-natural threats in the Shigar valley,(Central Karakorum) Baltistan-Pakistan. Pak. J. Bot. 52, 987–994 (2020).Article 

Google Scholar 
Anwar, S., et al., Plant diversity and communities pattern with special emphasis on the indicator species of a dry temperate forest: A case study from Liakot area of the Hindu Kush mountains, Pakistan. Trop. Ecol. 1–16 (2022).Mumshad, M. et al. Phyto-ecological studies and distribution pattern of plant species and communities of Dhirkot, Azad Jammu and Kashmir, Pakistan. PLoS ONE 16(10), e0257493 (2021).Article 

Google Scholar 
Baldeck, C. A. et al. Soil resources and topography shape local tree community structure in tropical forests. Proc. R. Soc. B Biol. Sci. 280(1753), 20122532 (2013).Article 

Google Scholar 
Guerra, T. N. F. et al. Influence of edge and topography on the vegetation in an Atlantic Forest remnant in northeastern Brazil. J. For. Res. 18(2), 200–208 (2013).Article 

Google Scholar 
Townsend, A. R., Asner, G. P. & Cleveland, C. C. The biogeochemical heterogeneity of tropical forests. Trends Ecol. Evol. 23(8), 424–431 (2008).Article 

Google Scholar 
Becknell, J. M. & Powers, J. S. Stand age and soils as drivers of plant functional traits and aboveground biomass in secondary tropical dry forest. Can. J. For. Res. 44(6), 604–613 (2014).Article 

Google Scholar 
Geri, F., Rocchini, D. & Chiarucci, A. Landscape metrics and topographical determinants of large-scale forest dynamics in a Mediterranean landscape. Landsc. Urban Plan. 95(1–2), 46–53 (2010).Article 

Google Scholar 
Lomolino, M. V. Elevation gradients of species-density: historical and prospective views. Glob. Ecol. Biogeogr. 10(1), 3–13 (2001).Article 

Google Scholar 
Zhang, K. et al. An integrated flood risk assessment approach based on coupled hydrological-hydraulic modeling and bottom-up hazard vulnerability analysis. Environ. Model. Softw. 148, 105279 (2022).Article 

Google Scholar 
Liu, Y. et al. A hybrid runoff generation modelling framework based on spatial combination of three runoff generation schemes for semi-humid and semi-arid watersheds. J. Hydrol. 590, 125440 (2020).Article 

Google Scholar 
Mir, A. Y. et al. Ethnopharmacology and phenology of high-altitude medicinal plants in Kashmir Northern Himalaya. Ethnobot. Res. Appl. 22, 1–15 (2021).
Google Scholar 
Vetaas, O. R. & Grytnes, J. A. Distribution of vascular plant species richness and endemic richness along the Himalayan elevation gradient in Nepal. Glob. Ecol. Biogeogr. 11(4), 291–301 (2002).Article 

Google Scholar 
Li, W. et al. Fine root biomass and morphology in a temperate forest are influenced more by the nitrogen treatment approach than the rate. Ecol. Ind. 130, 108031 (2021).Article 

Google Scholar 
Su, N. et al. Landscape context determines soil fungal diversity in a fragmented habitat. CATENA 213, 106163 (2022).Article 

Google Scholar 
Yang, Y., et al., Nitrogen fertilization weakens the linkage between soil carbon and microbial diversity: a global meta‐analysis. Global Change Biol. (2022).Ahmad, Z. et al. Weed species composition and distribution pattern in the maize crop under the influence of edaphic factors and farming practices: A case study from Mardan Pakistan. Saudi J. Biol. Sci. 23(6), 741–748 (2016).Article 

Google Scholar 
Rahman, A. U. et al. Ecological assessment of plant communities and associated edaphic and topographic variables in the Peochar Valley of the Hindu Kush mountains. Mt. Res. Dev. 36(3), 332–341 (2016).Article 

Google Scholar 
Ashton, P. S. A contribution of rain forest research to evolutionary theory. Ann. Mo. Bot. Gard. 64(4), 694–705 (1977).Article 

Google Scholar 
Yang, Y. et al. Negative effects of multiple global change factors on soil microbial diversity. Soil Biol. Biochem. 156, 108229 (2021).Article 

Google Scholar 
Pärtel, M. Local plant diversity patterns and evolutionary history at the regional scale. Ecology 83(9), 2361–2366 (2002).Article 

Google Scholar 
Taylor, D.R., Aarssen, L.W., & Loehle, C. On the relationship between r/K selection and environmental carrying capacity: A new habitat templet for plant life history strategies. Oikos 239–250 (1990).Knapp, A. K. et al. Rainfall variability, carbon cycling, and plant species diversity in a mesic grassland. Science 298(5601), 2202–2205 (2002).Article 

Google Scholar 
Zscheischler, J. et al. Short-term favorable weather conditions are an important control of interannual variability in carbon and water fluxes. J. Geophys. Res. Biogeosci. 121(8), 2186–2198 (2016).Article 

Google Scholar 
Gao, C. et al. Simulation and design of joint distribution of rainfall and tide level in Wuchengxiyu Region China. Urban Clim. 40, 101005 (2021).Article 

Google Scholar 
Wang, S. et al. Exploring the utility of radar and satellite-sensed precipitation and their dynamic bias correction for integrated prediction of flood and landslide hazards. J. Hydrol. 603, 126964 (2021).Article 

Google Scholar 
Zhang, K. et al. The sensitivity of North American terrestrial carbon fluxes to spatial and temporal variation in soil moisture: An analysis using radar-derived estimates of root-zone soil moisture. J. Geophys. Res. Biogeosci. 124(11), 3208–3231 (2019).Article 

Google Scholar 
Yang, Y. et al. Increasing contribution of microbial residues to soil organic carbon in grassland restoration chronosequence. Soil Biol. Biochem. 170, 108688 (2022).Article 

Google Scholar 
Li, J. et al. Differential mechanisms drive species loss under artificial shade and fertilization in the Alpine Meadow of the Tibetan Plateau. Front. Plant Sci. 13, 832473–832473 (2022).Article 

Google Scholar 
Fischer, C. et al. How do earthworms, soil texture and plant composition affect infiltration along an experimental plant diversity gradient in grassland?. PLoS ONE 9(6), e98987 (2014).Article 

Google Scholar 
Zhao, T. et al. Soil moisture experiment in the Luan River supporting new satellite mission opportunities. Remote Sens. Environ. 240, 111680 (2020).Article 

Google Scholar 
Marandi, A., Polikarpus, M. & Jõeleht, A. A new approach for describing the relationship between electrical conductivity and major anion concentration in natural waters. Appl. Geochem. 38, 103–109 (2013).Article 

Google Scholar 
Xu, J. et al. Modeling of coupled transfer of water, heat and solute in saline loess considering sodium sulfate crystallization. Cold Reg. Sci. Technol. 189, 103335 (2021).Article 

Google Scholar 
Chen, X. et al. Spatiotemporal characteristics and attribution of dry/wet conditions in the Weihe River Basin within a typical monsoon transition zone of East Asia over the recent 547 years. Environ. Model. Softw. 143, 105116 (2021).Article 

Google Scholar 
Ali, G., Siddique, S. & Suliman, M. Effect of canopy cover on natural regeneration of pinus wallichiana in moist temperate forest of Yakh Tangay, District Shangla Swat Pakistan. FUUAST J. Biol. 8(2), 193–201 (2018).
Google Scholar 
Zhang, K. et al. Characteristics and influencing factors of rainfall-induced landslide and debris flow hazards in Shaanxi Province, China. Nat. Hazard. 19(1), 93–105 (2019).Article 

Google Scholar 
Khan, W. et al. Vegetation mapping and multivariate approach to indicator species of a forest ecosystem: A case study from the Thandiani sub Forests Division (TsFD) in the Western Himalayas. Ecol. Ind. 71, 336–351 (2016).Article 

Google Scholar 
Iqbal, J. & Ahmed, M. Vegetation description of some pine forests of Shangla district of Khyber Pakhtunkhwa Pakistan: A preliminary study. FUUAST J. Biol. 4(1), 83–88 (2014).
Google Scholar 
Sparrow, B. D. et al. A vegetation and soil survey method for surveillance monitoring of rangeland environments. Front. Ecol. Evol. 8, 157 (2020).Article 

Google Scholar 
Esri, R., ArcGIS desktop: release 10. Environmental Systems Research Institute, CA (2011).Salzer, D., & Willoughby, J. Standardize this! The futility of attempting to apply a standard quadrat size and shape to rare plant monitoring. in Proceedings of the symposium of the North Coast Chapter of the California Native Plant Society: the ecology and management of rare plants of northwestern California. Arcata, CA. Sacramento, CA: The California Native Plant Society (2004).Bano, S. et al. Eco-Floristic studies of native plants of the Beer Hills along the Indus River in the districts Haripur and Abbottabad Pakistan. Saudi J. Biol. Sci. 25(4), 801–810 (2018).Article 

Google Scholar 
Perveen, A. & Qaiser, M. Pollen flora of Pakistan–XXXI Betulaceae. Pak. J. Bot. 31, 243–246 (1999).
Google Scholar 
Raunkiaer, C., The life forms of plants and statistical plant geography; being the collected papers of C. Raunkiaer. The life forms of plants and statistical plant geography; being the collected papers of C. Raunkiaer. (1934).Hussain, S.S., Pakistan manual of plant ecology1984: National Book Foundation.Kamran, S. et al. The role of graveyards in species conservation and beta diversity: A vegetation appraisal of sacred habitats from Bannu Pakistan. J. For. Res. 31(4), 1147–1158 (2020).Article 

Google Scholar 
Manan, F. et al. Environmental determinants of plant associations and evaluation of the conservation status of Parrotiopsis jacquemontiana in Dir, the Hindu Kush Range of Mountains. Trop. Ecol. 61(4), 509–526 (2020).Article 

Google Scholar 
Tfaily, M. M. et al. Sequential extraction protocol for organic matter from soils and sediments using high resolution mass spectrometry. Anal. Chim. Acta 972, 54–61 (2017).Article 

Google Scholar 
Chaney, R., Slonim, S., & Slonim, S. Determination of calcium carbonate content in soils, in Geotechnical properties, behavior, and performance of calcareous soils1982, ASTM International.McCune, B., & Mefford, M. PC-ORD, Multivariate analysis of ecological data, Version 5 for Windows edition. MjM Software Design, Gleneden Beach, Oregon USA (2005).Lepš, J., & Šmilauer, P. Multivariate analysis of ecological data using CANOCO2003: Cambridge university press.Xie, W. et al. A novel hybrid method for landslide susceptibility mapping-based geodetector and machine learning cluster: A case of Xiaojin county, China. ISPRS Int. J. Geo Inf. 10(2), 93 (2021).Article 

Google Scholar 
Li, L., Lei, Y. & Pan, D. Economic and environmental evaluation of coal production in China and policy implications. Nat. Hazards 77(2), 1125–1141 (2015).Article 

Google Scholar 
Team, R.C., R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/ (2013).Anwar, S. et al. Floristic composition and ecological gradient analyses of the Liakot Forests in the Kalam region of District Swat Pakistan. J. For. Res. 30(4), 1407–1416 (2019).Article 

Google Scholar 
Haq, S. M. et al. Exploring and understanding the floristic richness, life-form, leaf-size spectra and phenology of plants in protected forests: A case study of Dachigam National Park in Himalaya Asia. Acta Ecol. Sin. 41(5), 479–490 (2021).Article 

Google Scholar 
Ilyas, M. et al. A Preliminary checklist of the vascular flora of Kabal Valley, Swat Pakistan. Pak. J. Bot 45(2), 605–615 (2013).
Google Scholar 
Amjad, M. S. et al. Floristic composition, biological spectrum and phenological pattern of vegetation in the subtropical forest of Kotli District, AJK Pakistan. Pure Appl. Biol. (PAB) 6(2), 426–447 (2017).
Google Scholar 
Shaheen, H. et al. Species diversity, community structure, and distribution patterns in western Himalayan alpine pastures of Kashmir Pakistan. Mount. Res. Dev. 31(2), 153–159 (2011).Article 

Google Scholar 
Abbas, Z. et al. Ethnobotany of the balti community, tormik valley, karakorum range, baltistan, pakistan. J. Ethnobiol. Ethnomed. 12(1), 1–16 (2016).Article 

Google Scholar 
Ahmed, M. et al. Phytosociology and structure of Himalayan forests from different climatic zones of Pakistan. Pak. J. Bot. 38(2), 361 (2006).MathSciNet 

Google Scholar 
Shehzadi, S. et al. Floristic compositions along an 18-Km long transect in Ayubia National Park District Abbottabad Pakistan. Pak. J. Bot. 41(5), 2115–2127 (2009).
Google Scholar 
Khan, W., et al., Life forms, leaf size spectra and diversity indices of plant species grown in the Thandiani forests, district Abbottabad, Khyber Pakhtunkhwa, Pakistan. Saudi J. Biol. Sci.Kharkwal, G. et al. Phytodiversity and growth form in relation to altitudinal gradient in the Central Himalayan (Kumaun) region of India. Curr. Sci. 1, 873–878 (2005).
Google Scholar 
Bennie, J. et al. Slope, aspect and climate: Spatially explicit and implicit models of topographic microclimate in chalk grassland. Ecol. Model. 216(1), 47–59 (2008).Article 

Google Scholar 
Choudhary, K. & Nama, K. S. Phyto-diversity of Mukundara hills national park of Kota district, Rajasthan India. Adv. Appl. Sci. Res. 5(1), 18–23 (2014).
Google Scholar 
Shimwell, D.W., Description and classification of vegetation (1971).Malik, Z.H., Comparative study of vegetation of GungaChotti and Bedori Hills, Distric Bagh, Azad Jammu and Kashmir with special reference to range conditions, 2005, University of Peshawar, Pakistan.Khan, W. et al. Life forms, leaf size spectra, regeneration capacity and diversity of plant species grown in the Thandiani forests, district Abbottabad, Khyber Pakhtunkhwa Pakistan. Saudi J. Biol. Sci. 25(1), 94–100 (2018).Article 

Google Scholar 
Grytnes, J. A. & Vetaas, O. R. Species richness and altitude: A comparison between null models and interpolated plant species richness along the Himalayan altitudinal gradient Nepal. Am. Nat. 159(3), 294–304 (2002).Article 

Google Scholar 
Majid, A., Khan, M. & Calixto, E. Ecological assessment of plant communities along the edaphic and topographic gradients of biha valley, District Swat Pakistan. Appl. Ecol. Environ. Res. 16(5), 5611–5631 (2018).Article 

Google Scholar 
Khan, S.M., et al., Vegetation dynamics in the Western Himalayas, diversity indices and climate change. Sci. Tech. Dev. 31(3), 232–243 (2012).Khan, S. M. et al. Identifying plant species and communities across environmental gradients in the Western Himalayas: Method development and conservation use. Eco. Inform. 14, 99–103 (2013).Article 

Google Scholar 
Shaheen, H. & Shinwari, Z. K. Phyto diversity and endemic richness of Karambar lake vegetation from Chitral Hindukush-Himalayas. Pak. J. Bot 44(1), 17–21 (2012).
Google Scholar 
Wana, D., Plant communities and diversity along altitudinal gradients from Lake Abaya to Chencha Highlands, 2002, MA Thesis, School of Graduate Studies, Addis Ababa University. Addis Ababa.Canfora, L. et al. Is soil microbial diversity affected by soil and groundwater salinity? Evidences from a coastal system in central Italy. Environ. Monit. Assess. 189(7), 1–15 (2017).Article 

Google Scholar 
Liu, S., et al., The distribution characteristics and human health risks of high-fluorine groundwater in coastal plain: A case study in Southern Laizhou Bay, China. Front. Environ. Sci. 568 (2022).Niu, Y. et al. Vegetation distribution along mountain environmental gradient predicts shifts in plant community response to climate change in alpine meadow on the Tibetan Plateau. Sci. Total Environ. 650, 505–514 (2019).Article 

Google Scholar 
Nadal-Romero, E. et al. Effects of slope angle and aspect on plant cover and species richness in a humid Mediterranean badland. Earth Surf. Proc. Land. 39(13), 1705–1716 (2014).Article 

Google Scholar  More

Huijbers, P. M. C., Flach, C.-F. & Larsson, D. G. J. A conceptual framework for the environmental surveillance of antibiotics and antibiotic resistance. Environ. Int. 130, 104880 (2019).Article 

Google Scholar 
Aarestrup, F. M. & Woolhouse, M. E. J. Using sewage for surveillance of antimicrobial resistance. Science 367, 630–632 (2020).Article 

Google Scholar 
European Commission. Proposal for a revised Urban Wastewater Treatment Directive. European Commission https://environment.ec.europa.eu/publications/proposal-revised-urban-wastewater-treatment-directive_en (2022).US Centres for Disease Control and Prevention. COVID-19 impacts on environment (e.g., water, soil) and sanitation: addressing antimicrobials and antimicrobial resistant threats in the environment. US Centres for Disease Control and Prevention https://www.cdc.gov/drugresistance/pdf/covid19/COVID19-Impacts-AR-Environment-Sanitation-508.pdf (2021).Flach, C.-F., Hutinel, M., Razavi, M., Åhrén, C. & Larsson, D. G. J. Monitoring of hospital sewage shows both promise and limitations as an early-warning system for carbapenemase-producing Enterobacterales in a low-prevalence setting. Water Res. 200, 117261 (2021).Article 

Google Scholar 
Larsson, D. G. J. & Flach, C.-F. Antibiotic resistance in the environment. Nat. Rev. Microbiol. 20, 257–269 (2022).Article 

Google Scholar 
Newton, R. J. et al. Sewage reflects the microbiomes of human populations. mBio 6, e02574 (2015).Article 

Google Scholar 
Huijbers, P. M. C., Larsson, D. G. J. & Flach, C. F. Surveillance of antibiotic resistant Escherichia coli in human populations through urban wastewater in ten European countries. Environ. Pollut. 261, 114200 (2020).Article 

Google Scholar 
Laxminarayan, R. & Macauley, M. K. The Value of Infromation: Methodological Frontiers and New Applications in Environment and Health 1st edn (Springer Dordrecht, 2012).Munk, P. et al. Genomic analysis of sewage from 101 countries reveals global landscape of antimicrobial resistance. Nat. Commun. 13, 7251 (2022).Article 

Google Scholar  More

Global statistics on declining biodiversity can give the impression that every population of every species is in a downward spiral. In fact, many populations are stable or growing, while a small number of species faces truly existential challenges. These photos capture some specific crises. They are images of threats unfolding, of desperate attempts at species defence and of the beautiful living world that is at stake.
The 15th United Nations Biodiversity Conference, COP15, opens in Montreal, Canada, on 7 December. At the meeting, delegates will attempt to agree on goals for stabilizing species’ declines by 2030 and reverse them by mid-century. The current draft framework agreement promises nothing less than a “transformation in society’s relationship with biodiversity”.
Help for the kelp. Tasmania’s forests of giant kelp (Macrocystis pyrifera) are dying as climate change shifts ocean currents, bringing warm water to the east coast of the temperate Australian island. The kelp forests host an entire ecosystem, including abalone and crayfish — both economically important species and part of local food culture. Now, researchers at the Institute for Marine and Antarctic Studies in Hobart are breeding kelp plants that can tolerate warmer conditions, and replanting them along the coast — a trial for what they hope will become a landscape-scale restoration. More

Ah-King, M. Flexible mate choice in Encyclopedia of Animal Behavior, 2nd edn Vol. 4 (ed Jae Chun Choe) 421–431 (Academic Press, 2019).Harestad, A. S. & Bunnell, F. L. Home range and body weight—A reevaluation. Ecology 60, 389–402 (1979).Article 

Google Scholar 
O’Neill, R. V., Milne, B. T., Turner, M. G. & Gardner, R. H. Resource utilization scales and landscape pattern. Landsc. Ecol. 2, 63–69 (1988).Article 

Google Scholar 
Tricas, T. C. Determinants of feeding territory size in the corallivorous butterflyfish, Chaetodon multicinctus. Anim. Behav. 37, 830–841. https://doi.org/10.1016/0003-3472(89)90067-5 (1989).Article 

Google Scholar 
Tremblay, I., Thomas, D., Blondel, J., Perret, P. & Lambrechts, M. M. The effect of habitat quality on foraging patterns, provisioning rate and nestling growth in Corsican Blue Tits Parus caeruleus. Ibis 147, 17–24. https://doi.org/10.1111/j.1474-919x.2004.00312.x (2005).Article 

Google Scholar 
Watts, D. P. The influence of male mating tactics on habitat use in Mountain Gorillas (Gorilla gorilla beringei). Primates 35, 35–47. https://doi.org/10.1007/BF02381484 (1994).Article 

Google Scholar 
Lescroël, A. et al. Working less to gain more: when breeding quality relates to foraging efficiency. Ecology 91, 2044–2055. https://doi.org/10.1890/09-0766.1 (2010).Article 
PubMed 

Google Scholar 
Tufto, J., Anderson, R. & Linnell, J. Habitat use and ecological correlates of home range size in a small cervid: the roe deer. J. Anim. Ecol. 65, 715–724. https://doi.org/10.2307/5670 (1996).Article 

Google Scholar 
Morellet, N. et al. Seasonality, weather and climate affect home range size in roe deer across a wide latitudinal gradient within Europe. J. Anim. Ecol. 82, 1326–1339. https://doi.org/10.1111/1365-2656.12105 (2013).Article 
PubMed 

Google Scholar 
Anderson, D. P. et al. Scale-dependent summer resource selection by reintroduced elk in Wisconsin, USA. J. Wildl. Manag. 69, 298–310. https://doi.org/10.2193/0022-541X(2005)069%3c0298:SSRSBR%3e2.0.CO;2 (2005).Article 

Google Scholar 
Olsson, P. M. O. et al. Movement and activity patterns of translocated elk (Cervus elaphus nelsoni) on an active coal mine in Kentucky. Wildl. Biol. Pract. 3, 1–8. https://doi.org/10.2461/wbp.2007.3.1 (2007).Article 

Google Scholar 
Porter, W. P., Sabo, J. L., Tracy, C. R., Reichman, O. J. & Ramankutty, N. Physiology on a landscape scale: plant–animal interactions. Integr. Comp. Biol. 42, 431–453. https://doi.org/10.1093/icb/42.3.431 (2002).Article 
PubMed 

Google Scholar 
Berg, J. E. et al. Mothers’ movements: shifts in calving area selection by partially migratory elk. J. Wildl. Manag. 85, 1476–1489. https://doi.org/10.1002/jwmg.22099 (2021).Article 

Google Scholar 
Lehman, C. P. et al. Elk resource selection at parturition sites, Black Hills, South Dakota. J. Wildl. Manag. 80, 465–478. https://doi.org/10.1002/jwmg.1017 (2016).Article 

Google Scholar 
Johnson, B. K., Kern, J. W., Wisdom, M. J., Findholt, S. L. & Kie, J. G. Resource selection and spatial separation of mule deer and elk during spring. J. Wildl. Manag. 64, 685–697. https://doi.org/10.2307/3802738 (2000).Article 

Google Scholar 
Grace, J. & Easterbee, N. The natural shelter for red deer (Cervus elaphus) in a Scottish glen. J. Appl. Ecol. 16, 37–48. https://doi.org/10.2307/2402726 (1979).Article 

Google Scholar 
Demarchi, M. W. & Bunnell, F. L. Estimating forest canopy effects on summer thermal cover for Cervidae (deer family). Can. J. For. Res. 23, 2419–2426. https://doi.org/10.1139/x93-299 (1993).Article 

Google Scholar 
Parker, K. L. & Gillingham, M. P. Estimates of critical thermal environments for mule deer. J. Range. Manag. 43, 73–81 (1990).Article 

Google Scholar 
Proffitt, K. M. et al. Changes in elk resource selection and distributions associated with a late-season elk hunt. J. Wildl. Manag. 74, 210–218. https://doi.org/10.2193/2008-593 (2010).Article 

Google Scholar 
Webb, S. L., Dzialak, M. R., Harju, S. M., Hayden-Wing, L. D. & Winstead, J. B. Influence of land development on home range use dynamics of female elk. Wildl. Res. 38, 163–167. https://doi.org/10.1071/WR10101 (2011).Article 

Google Scholar 
Rumble, M. A., Benkobi, L. & Gamo, R. S. Elk responses to humans in a densely roaded area. Intermt. J. Sci. 11, 10–24 (2005).
Google Scholar 
McCorquodale, S. M. Sex-specific movements and habitat use by elk in the Cascade Range of Washington. J. Wildl. Manag. 67, 729–741. https://doi.org/10.1890/15-1607.1 (2003).Article 

Google Scholar 
Saïd, S. & Servanty, S. The influence of landscape structure on female roe deer home-range size. Landsc. Ecol. 20, 1003–1012. https://doi.org/10.1007/s10980-005-7518-8 (2005).Article 

Google Scholar 
Seddon, P. J., Armstrong, D. P. & Maloney, R. F. Developing the science of reintroduction biology. Conserv. Biol. 21, 303–312. https://doi.org/10.1111/j.1523-1739.2006.00627.x (2007).Article 
PubMed 

Google Scholar 
Hale, S. L. & Koprowski, J. L. Ecosystem-level effects of keystone species reintroduction: a literature review. Restor. Ecol. 26, 439–445. https://doi.org/10.1111/rec.12684 (2018).Article 

Google Scholar 
Cheyne, S. M. Wildlife reintroduction: considerations of habitat quality at the release site. BMC Ecol. 6, 5. https://doi.org/10.1186/1472-6785-6-5 (2006).Article 
PubMed 
PubMed Central 

Google Scholar 
Hegel, T. M., Gates, C. C. & Eslinger, D. The geography of conflict between elk and agricultural values in the Cypress Hills, Canada. J. Eniron. Manag. 90, 222–235. https://doi.org/10.1016/j.jenvman.2007.09.005 (2009).Article 

Google Scholar 
Walter, W. D. et al. Management of damage by elk (Cervus elaphus) in North America: a review. Wildl. Res. 37, 630–646. https://doi.org/10.1071/WR10021 (2010).Article 

Google Scholar 
Jung, T. S. Extralimital movements of reintroduced bison (Bison bison): implications for potential range expansion and human–wildlife conflict. Eur. J. Wildl. Res. 63, 35. https://doi.org/10.1007/s10344-017-1094-5 (2017).Article 

Google Scholar 
Buchholtz, E. K., Stronza, A., Songhurst, A., McCulloch, G. & Fitzgerald, L. A. Using landscape connectivity to predict human-wildlife conflict. Biol. Conserv. 248, 108677. https://doi.org/10.1016/j.biocon.2020.108677 (2020).Article 

Google Scholar 
Hodgson, J. A., Moilanen, A., Wintle, B. A. & Thomas, C. D. Habitat area, quality and connectivity: striking the balance for efficient conservation. J. Appl. Ecol. 48, 148–152. https://doi.org/10.1111/j.1365-2664.2010.01919.x (2011).Article 

Google Scholar 
Murie, O. The Elk of North America (Stackpole Co., 1951).
Google Scholar 
VDWR. Virginia elk management plan 2019–2028 (ed Virginia Department of Wildlife Resources) (Virginia Department of Wildlife Resources, 2019).Lituma, C. M. et al. Terrestrial wildlife in the post-mined Appalachian landscape: status and opportunities in Appalachia’s Coal-Mined Landscapes (eds Carl E. Zipper & Jeff Skousen) 135–166 (Springer, 2021).Lupardus, J. L., Muller, L. I. & Kindall, J. L. Seasonal forage availability and diet for reintroduced elk in the Cumberland Mountains, Tennessee. Southeast. Nat. 10, 53–74. https://doi.org/10.1656/058.010.0105 (2011).Article 

Google Scholar 
Schneider, J. et al. Food habits of reintroduced elk in southeastern Kentucky. Southeast. Nat. 5, 535–546. https://doi.org/10.1656/1528-7092(2006)5[535:Fhorei]2.0.Co;2 (2006).Article 

Google Scholar 
Smith, T. N., Keller, B. J., Chitwood, M. C., Hansen, L. P. & Millspaugh, J. J. Diet composition and selection of recently reintroduced elk in Missouri. Am. Midl. Nat. 180, 143–159. https://doi.org/10.1674/0003-0031-180.1.143 (2018).Article 

Google Scholar 
Franklin, J. A., Zipper, C. E., Burger, J. A., Skousen, J. G. & Jacobs, D. F. Influence of herbaceous ground cover on forest restoration of eastern US coal surface mines. New. For. 43, 905–924. https://doi.org/10.1007/s11056-012-9342-8 (2012).Article 

Google Scholar 
Popp, J. N., Toman, T., Mallory, F. F. & Hamr, J. A century of elk restoration in eastern North America. Restor. Ecol. 22, 723–730. https://doi.org/10.1111/rec.12150 (2014).Article 

Google Scholar 
Cook, J. G., Irwin, L. L., Bryant, L. D., Riggs, R. A. & Thomas, J. W. Relations of forest cover and condition of elk: a test of the thermal cover hypothesis in the summer and winter. Wildl. Monogr. 141, 3–61 (1998).
Google Scholar 
Parker, K. L. & Robbins, C. T. Thermoregulation in mule deer and elk. Can. J. Zool. 62, 1409–1422. https://doi.org/10.1139/z84-202 (1984).Article 

Google Scholar 
Mao, J. S. et al. Habitat selection by elk before and after wolf reintroduction in Yellowstone National Park. J. Wildl. Manag. 69, 1691–1707. https://doi.org/10.2193/0022-541X (2005).Article 

Google Scholar 
Wolff, J. O. & Van Horn, T. Vigilance and foraging patterns of American elk during the rut in habitats with and without predators. Can. J. Zool. 81, 266–271. https://doi.org/10.1139/z03-011 (2003).Article 

Google Scholar 
Beck, J. L. & Peek, J. M. Diet composition, forage selection, and potential for forage competition among elk, deer, and livestock on aspen–sagebrush summer range. Rangel. Ecol. Manag. 58, 135–147. https://doi.org/10.2111/03-13.1 (2005).Article 

Google Scholar 
Ford, W. M., Johnson, A. S. & Hale, P. E. Nutritional quality of deer browse in southern Appalachian clearcuts and mature forests. For. Ecol. Manag. 67, 149–157. https://doi.org/10.1016/0378-1127(94)90013-2 (1994).Article 

Google Scholar 
Sikes, R. S., Gannon, W. L. & The American Care and Use Committee of the American Society of Mammalogists. Guidelines of the American Society of Mammalogists for the use of wild mammals in research. J. Mammal. 92, 235–253. https://doi.org/10.1644/10-mamm-f-355.1 (2011).Percie du Sert, N. et al. The ARRIVE guidelines 2.0: updated guidelines for reporting animal research. PLoS Biol. 18, e3000410. https://doi.org/10.1371/journal.pbio.3000410 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Powell, J. W. Physiographic Regions of the United States. (American Book Company, 1895).Braun, E. L. Forests of the Cumberland Mountains. Ecol. Monogr. 12, 413–447. https://doi.org/10.2307/1943039 (1942).Article 

Google Scholar 
Clark, J. B. The Vascular Flora of Breaks Interstate Park, Pike County, Kentucky, and Dickenson County, Virginia Master of Science thesis, Eastern Kentucky University (2012).Pericak, A. A. et al. Mapping the yearly extent of surface coal mining in Central Appalachia using Landsat and Google Earth Engine. PLoS ONE 13, e0197758. https://doi.org/10.1371/journal.pone.0197758 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Boettner, F. et al. An assessment of the natural assets in the Appalachian Region: forest resources (ed Appalachian Regional Commission Report) 97 (Washington, DC, 2014).NOAA. Summary of monthly normals Grundy, VA 1991 – 2020 data (National Oceanic and Atmospheric Administration (2022).U.S. Geological Survey (USGS) Gap Analysis Project (GAP). GAP/LANDFIRE national terrestrial ecosystems 2011: U.S. Geological Survey data release (2016).Clark, M. The Nature Conservancy Eastern Division & North Atlantic Landscape Conservation Cooperative. Terrestrial habitat, Northeast data (2017).ESRI. ArcGIS desktop version 10.8.1 (Environmental Systems Research Institute, 2020).Ford, W. M. et al. Influence of elevation and forest type on community assemblage and species distribution of shrews in the central and southern Appalachians in Advances in the Biology of the Shrews II Vol. 1(eds. J.F. Merritt, S. Churchfield, R. Hutterer and B.A. Sheftel) 303–315(Special Publication of the International Society of Shrew Biologists, 2006).Kniowski, A. B. & Ford, W. M. Predicting intensity of white-tailed deer herbivory in the Central Appalachian Mountains. J. For. Res. 29, 841–850. https://doi.org/10.1007/s11676-017-0476-6 (2018).Article 

Google Scholar 
Fleming, C. H. & Calabrese, J. M. ctmm: continuous-time movement modeling. R package version 0.6.0 (2021).R Core Team. R: a language and environment for statistical computing (R Foundation for Statistical Computing, 2020).Fleming, C. H. et al. Estimating where and how animals travel: an optimal framework for path reconstruction from autocorrelated tracking data. Ecology 97, 576–582. https://doi.org/10.1890/15-1607.1 (2016).Article 
PubMed 

Google Scholar 
Hijmans, R. J. raster: geographic data analysis and modeling. R package version 3.4-5 (2020).Becker, R. A., Chambers, J. M. & Wilks, A. R. The New S Language (Wadsworth and Brooks/Cole, 1988).Kuznetsova, A., Brockhoff, P. B. & Christensen, R. B. H. lmerTest package: tests in linear mixed effects models. J. Stat. Softw. 82, 1–26. https://doi.org/10.18637/jss.v082.i13 (2017).Article 

Google Scholar 
Burnham, K. P. & Anderson, D. R. Model Selection and Inference: A Practical Use of the Information-Theoretic Approach (Springer, 1998).Turner, M. G., Wu, Y., Romme, W. H. & Wallace, L. L. A landscape simulation model of winter foraging by large ungulates. Ecol. Modell. 69, 163–184. https://doi.org/10.1016/0304-3800(93)90026-O (1993).Article 

Google Scholar 
Taper, M. L. & Gogan, P. J. P. The northern Yellowstone elk: density dependence and climatic conditions. J. Wildl. Manag. 66, 106–122. https://doi.org/10.2307/3802877 (2002).Article 

Google Scholar 
Green, R. A. & Bear, G. D. Seasonal cycles and daily activity patterns of Rocky Mountain elk. J. Wildl. Manag. 54, 272–279. https://doi.org/10.2307/3809041 (1990).Article 

Google Scholar 
Craighead, J. J., Craighead, F. C. J., Ruff, R. L. & O’Gara, B. W. Home ranges and activity patterns of nonmigratory elk of the Madison Drainage herd as determined by biotelemetry. Wildl. Monogr. 33, 3–50 (1973).
Google Scholar 
Gittleman, J. L. & Thompson, S. D. Energy allocation in mammalian reproduction. Am. Zool. 28, 863–875. https://doi.org/10.1093/icb/28.3.863 (1988).Article 

Google Scholar 
Beier, P. & McCullough, D. R. Factors influencing white-tailed deer activity patterns and habitat use. Wildl. Monogr. 109, 3–51 (1990).
Google Scholar 
Ciuti, S., Davini, S., Luccarini, S. & Apollonio, M. Variation in home range size of female fallow deer inhabiting a sub-Mediterranean habitat. Rev. Ecol. 58, 381–395 (2003).
Google Scholar 
Vore, J. M. & Schmidt, E. M. Movements of female elk during calving season in northwest Montana. Wildl. Soc. Bull. 29, 720–725 (2001).
Google Scholar 
Wickstrom, M. L., Robbins, C. T., Hanley, T. A., Spalinger, D. E. & Parish, S. M. Food intake and foraging energetics of elk and mule deer. J. Wildl. Manag. 48, 1285–1301. https://doi.org/10.2307/3801789 (1984).Article 

Google Scholar 
Van Soest, P. J. Allometry and ecology of feeding behavior and digestive capacity in herbivores: a review. Zoo. Biol. 15, 455–479 (1996). https://doi.org/10.1002/(SICI)1098-2361(1996)15:53.0.CO;2-AEsmaeili, S. et al. Body size and digestive system shape resource selection by ungulates: a cross-taxa test of the forage maturation hypothesis. Ecol. Lett. 24, 2178–2191. https://doi.org/10.1111/ele.13848 (2021).Article 
PubMed 

Google Scholar 
Demment, M. W. & Van Soest, P. J. A nutritional explanation for body-size patterns of ruminant and nonruminant herbivores. Am. Nat. 125, 641–672 (1985).Article 

Google Scholar 
Anderson, D. P. et al. Factors influencing female home range sizes in elk (Cervus elaphus) in North American landscapes. Landsc. Ecol. 20, 257–271. https://doi.org/10.1007/s10980-005-0062-8 (2005).Article 

Google Scholar 
Maigret, T. A., Cox, J. J. & Yang, J. Persistent geophysical effects of mining threaten ridgetop biota of Appalachian forests. Front. Ecol. Environ. 17, 85–91. https://doi.org/10.1002/fee.1992 (2019).Article 

Google Scholar 
Beier, P. Sex differences in quality of white-tailed deer diets. J. Mammal. 68, 323–329. https://doi.org/10.2307/1381471 (1987).Article 

Google Scholar 
Parker, K. L., Barboza, P. S. & Gillingham, M. P. Nutrition integrates environmental responses of ungulates. Funct. Ecol. 23, 57–69. https://doi.org/10.1111/j.1365-2435.2008.01528.x (2009).Article 

Google Scholar 
Wichrowski, M. W., Maehr, D. S., Larkin, J. L., Cox, J. J. & Olsson, M. P. O. Activity and movements of reintroduced elk in southeastern Kentucky. Southeast. Nat. 4, 365–374. https://doi.org/10.1656/1528-7092(2005)004[0365:Aamore]2.0.Co;2 (2005).Article 

Google Scholar 
Relyea, R. A., Lawrence, R. K. & Demarais, S. Home range of desert mule deer: testing the body-size and habitat-productivity hypotheses. J. Wildl. Manag. 64, 146–153. https://doi.org/10.2307/3802984 (2000).Article 

Google Scholar  More

It is well known that parrots are excellent vocal learners; here we quantified that ability across a wide variety of species, using human mimicry as a proxy for vocal learning of natural repertoires. Results confirm that parrot vocal mimicry varies substantially both within and among species22. Parrot age, social interactions, and sex do not appear to be universal drivers of vocal learning ability within the order Psittaciformes, but all of these factors may have effects within individual species.Vocal learning variation by speciesWithin species, mimicry sound repertoires are extremely variable bird to bird; for example, our data indicate that a grey parrot may mimic anywhere from 0 to 600 different human words. Many other species showed smaller repertoires but similar variability. It is not entirely clear whether this range of variation would be present in natural sounds within wild parrot populations, but research has demonstrated intraspecific repertoire size variation in multiple species of parrots30,31.The vast majority of parrots presented a pattern in which their repertoire size was largest for words, intermediate for phrases (composed of the reported words), and smallest for non-linguistic sounds (Fig. 2). In the wild, parrots mimic the most socially relevant vocalizations, and presumably do so in captivity as well15. Thus, the spoken word and phrase interactions with their human “flock” likely reflect the most socially relevant cues. The interesting exceptions to this pattern were Fischer’s lovebirds, cockatiels, and Senegal parrots who all used more sounds than phrases. Cockatiels are well-known in the pet world to be excellent whistlers, and thus it was satisfying to see that our data support that informal information. We suspect that deviations from the typical patterns may represent acoustic learning preferences, templates, or limitations32.Although individual variation was substantial, we nevertheless saw strong evidence that overall vocal learning abilities differed by species. Pacific parrotlets and sun parakeets showed very limited human mimicry, while grey parrots, Amazona parrots, cockatoos, and macaws were generally very accomplished mimics. The patterns that we documented appeas to reflect natural vocal repertoire variation across species. The documented calls of wild parrots generally range from 5 to 15 calls25,33,34,35,36. Several species, however, present additional complexity: yellow-naped parrots (Amazona auropalliata), palm cockatoos (Probosciger aterrimus), and grey parrots all have natural repertoires of more than 25 discrete elements, with additional elements given in duets13,27,37 Members of these three groups, grey parrots, Amazona parrots and cockatoos also had relatively large repertoires in our study. In several of these species (particularly grey parrots) our measure of mimicked “words” (60) was higher than estimates of natural call “elements” (39) in the literature27. This discrepancy suggests that parrots are capable of learning vocalizations with more than 25 elements and, simultaneously, might reflect a sampling bias wherein survey-takers are more likely to report on individuals with high mimicry ability.Parrot species varied in their tendency to improvise new combinations of elements, although most species did rearrange words to some degree. Research shows that parrot vocalization length and structure carry signal content, so there may be selective pressures favoring this ability24,33. If so, then our data suggest that those pressures are strongest in some cockatoos and weakest in sun parakeets and green-cheeked parakeets. In general, species with larger repertoires also showed more vocal flexibility (Fig. 2, Appendix 6). Additionally, wild birds typically use particular vocalizations in set contexts, so the ability to do so is likely to be adaptive24. Previous studies of captive parrots have demonstrated contextual use of mimicked words, both in tutored lab settings and in home-raised birds28,38. In our sample, contextual use of learned sounds was supported across 89% of individuals and most species. Survey-taker responses on this topic are necessarily subjective, so we emphasize that this rate of contextual use should be interpreted as a general estimate. Nevertheless, the data indicated that parrots frequently associated mimicked human sounds with appropriate human contexts. This finding is particularly revealing because the relevant human contexts are, by their nature, outside the range of typical wild parrot experiences. Contextual vocalization use must, therefore, rely on extremely flexible vocal learning mechanisms.Vocal learning variation by ageOn average, birds aged with high confidence were younger than those aged with low or medium confidence. This pattern might indicate that people tend to overestimate the age of captive birds of uncertain age. This pattern might also reflect the facts that older birds are more likely to be wild-caught and that younger birds are more likely to have good hatch-date documentation. In either case, there are few ramifications of inaccurate age estimates relating to vocal behavior because our data gave no evidence that adult vocal mimicry repertoires varied with age. Our analyses of grey parrots confirmed that repertoires expanded through the juvenile phase, but did not show reliable expansion among adults. Studies of wild birds indicate that parrots can learn vocalizations throughout life; such open-ended learning is limited to a subset of vocal learning species, and can generate different outcomes as animals age15. In some species, animals can add new vocal features over the course of a lifetime, leading to repertoire expansion39,40. In other species, animals may replace parts of their repertoire with newly-learned vocalizations, leading to stable vocal production repertoire sizes across age groups39,41. Our data suggest that parrots fit the second pattern; although they are open-ended vocal learners, their adult repertoires change more by element replacement, than by expansion. This does not necessarily imply that vocalizations are “forgotten” through time, but merely that some sounds are no longer used as conditions change42. Many parrot vocalizations function in social coordination with flock-mates22. The fission–fusion nature of parrot flocks creates changing social conditions for each individual over its lifetime43. A vocal replacement model for repertoire learning would allow individuals to adjust their vocal signatures to match new social situations and stop producing vocalizations that are no longer socially relevant11,44.Vocal learning variation by sexOur analyses of the full data set confirmed the generally held understanding that males and females in most species of parrots have similar vocal learning abilities15. We did, however see sex differences in some species that merit future study. First, we found a substantial overrepresentation of males in our sample. This could be interpreted several ways; (1) there are legitimately more males in the parrot pet trade, (2) pet owners are giving us accurate data but are more likely to give us data on males or (3) some bias exists in which pet owners assume their talking parrots are males, rather than females. Possibilities 1 and 2 seem unlikely because after we eliminated all parrots sexed with low confidence, we were left with a nearly 1:1 ratio of males:females in the subset of parrots that were sexed with high confidence. That trend suggests that the male bias in our data comes (at least in part) from a human tendency to label their pet parrots as male when the sex is not clear. Among songbirds, there is a strong tendency to assume that singing birds are male, and a similar bias may hold true for parrots45. It is unclear whether parrots in this study were mislabeled as male because they vocalize or, more simply, because that is the default human tendency for any animal.Although we conclude that some of the male bias in our data is human error, we also saw patterns that suggest real sex differences in vocal learning some species. For example, Pacific parrotlets are a dimorphic species, and all of our sampled birds were sexed by plumage46. Thus, we expect sexing in this species to be fairly accurate. Our data set included 10 males and no females, a bias unlikely to result purely from sampling error. We saw a similar trend in cockatiels for which there was a large overabundance of males in the data set, even among the 17 birds sexed with high confidence. Humans may be more likely to report on parrots that are good mimics. Therefore, the results likely reflect a real-world tendency for male cockatiels to mimic more human sounds than females. Figure 3 suggests that the same might be true for galahs, sulphur-crested cockatoos, rose-ringed parakeets, Senegal parrots, and budgerigars. Existing research supports the idea that sex differences in vocal behavior are important in several of these species. Among galahs, male and female calls evoke different responses47, and patterns of call adjustment vary by sex among budgerigars20. We also note that several of these species (Pacific parrotlets, rose-ringed parakeets, budgerigars, and cockatiels; Appendix 2b) show sex-based differences in both plumage and vocal learning, raising questions about whether those traits co-evolve.In addition to sex-based differences in the tendency to mimic humans, several well-sampled species showed evidence of sex-based differences in repertoire sizes. Particularly interesting are the blue-and-yellow macaws, in which repertoire size was significantly male-biased. We had more females (15) than males (9) in the data set, but males used on average 3–4 times as many mimicry sounds, phrases and words as females did. Galahs and budgerigars showed a similar male-bias in repertoire sizes, matching the trend of males being overrepresented in our data set for those two species. Prior research on galahs and budgerigars has found that males can be more vocal and more flexible with their vocalizations; perhaps these abilities translate to learning more call types20,47. A similar, but weaker, male mimicry increase occurred in rose-ringed parakeets. In only one species, yellow-headed parrots, did females show a significantly larger mimicry repertoire than males in any category (Appendix 5). Interestingly, the tendency to mimic humans (measured as sampling in the data set) and repertoire sizes did not always show the same patterns. Among sulphur-crested cockatoos, cockatiels, and Senegal parrots, males were more likely to show human mimicry, but their repertoires were not larger than the repertoires of females. This suggests that in some species, females may be less likely to mimic vocalizations, but when they do so they have just as large a vocabulary as males.The reported sex differences in parrot vocal mimicry repertoires are intriguing, but also are tentative conclusions. In many species, including our best sampled species, grey parrots, we saw no evidence of sex-differences in repertoire size. The sex-biases that we did document lose statistical significance after controlling for the many comparisons that we conducted. Nevertheless, we expect that some of our data represent true biological differences, especially because studies of wild birds have shown similar trends47,48. Thus, we offer our data as a starting point for additional research. Taken together, the analyses by sex provide interesting points of comparison to other vocal learning animals. Our combined analyses suggest that sex differences in vocal learning are vastly smaller and less common among parrots than they are among oscine passerines and hummingbirds45,49,50. Sex-based patterns of vocal learning in parrots appear more similar to those of vocal learning mammals than to those of other vocal learning birds51. Overall, parrots and songbirds present excellent comparative study systems for all aspects of sex differences in song learning, from the mechanistic to the functional17,51.Vocal learning variation by social contextMany parrot vocalizations function in social organization for individuals within flocks, and the ability to learn from conspecifics is essential to parrot familial and social integration12,15,52. Although our study specifically examined vocal learning of human sounds, we thought it possible that the presence of other parrots would increase mimicry rates if parrots learned human vocalizations from their parrot companions. Anecdotal stories of parrots teaching words to other parrots abound53, and studies of grey parrot cognition show that vocal modeling by multiple tutors can lead to better learning of human words54. Most existing results, however, are based on human tutoring, with controlled studies of parrot-parrot word transmission lacking. Here we tested whether social interactions with other parrots correlated with more vocal learning of human sounds. Our data gave no evidence that parrot-parrot social interactions drive human vocal mimicry. This was true across the full sample (controlling for species identity), and for our best sampled species, grey parrots. Although companion parrots are known to learn from conspecifics, that learning does not appear to shape repertoire sizes53. Open questions remain about whether signal complexity, repertoire size, or aspects of vocal learning covary with social complexity at a larger scale among parrots55. Follow up studies should address these questions using phylogenetically-controlled methods56. More

To estimate the total number of non-broadcasting vessels, including those that were not detected by SAR, we: (1) obtained SAR detections of vessels from RADARSAT-2 and the corresponding vessel lengths as estimated from the SAR image; (2) processed a global feed of AIS data to identify every broadcasting vessel that should have appeared in the SAR images at the moment the images were taken; (3) developed a novel technique to determine which vessels in AIS matched to detections in SAR, which AIS vessels were not detected by SAR, and which SAR detections represented non-broadcasting vessels; (4) after matching SAR to AIS, we could then (a) model the relationship between a vessel’s actual length and the length as estimated by the SAR image (Fig. 3b) and (b) model the relationship between the likelihood that a vessel is detected and its length (Fig. 3a); and (5) finally, we combined these relationships to develop an estimate of the number and lengths of non-broadcasting vessels in the region.SAR imagery and vessel detectionsWorking with the satellite company Kongsberg Satellite Services (KSAT), we tasked the Canadian Space Agency’s satellite RADARSAT-2 to acquire SAR images from its ship detection mode (DVWF mode, GRD product), with a pixel size of about 40 m and a swath width over 400 km (19). These images were processed following standard procedures for GRD products (e.g. applying radiometric calibration and geometric corrections)29,30. Vessel locations were extracted from the images with the widely used ship detection algorithms, which discriminates objects at sea based on the backscatter difference (pixel values) between the sea clutter and the targets31. Vessel lengths were estimated by measuring distances directly on the images with the aid of a graphical user interface tool31.Identifying Vessels using AISIn each region, AIS data, obtained from satellite providers ORBCOMM and Spire, were processed using Global Fishing Watch’s data pipeline1. The identities and lengths of all AIS devices that operated near the SAR scenes in both space and time were first obtained using Global Fishing Watch’s database1. To be sure vessels were identified correctly, two analysts reviewed the tracks of every AIS device in each region.In both regions, it is common practice for fishers to put AIS beacons on their longlines, likely to aid in retrieving them, meaning that many AIS devices were longline gear and not vessels. Because gear outnumbered vessels by several-fold, it was critical to differentiate gear and fishing vessels. In the Indian Ocean, 521 unique AIS devices associated with gear were detected that were likely within the SAR scenes, and 390 unique AIS devices associated with gear in the Pacific that were likely within the SAR scenes. Transponders were determined to be associated with gear by inspecting the name broadcast in the AIS messages (gear frequently broadcasts one of several standard names and/or a voltage reading) and classification using the Global Fishing Watch vessel classification algorithm1. Most gear also had an MMSI number (unique identifier number for AIS) that started with 1, 8, or 9 or broadcast names that signified gear. We eliminated all gear from the analysis because (1) these gear buoys have reflectors that are only ~ 1 m in size, and they should not be visible in ~ 40 m resolution SAR images, and (2) we found that gear matched to SAR detections only when traveling faster than 2 knots (and thus was on the deck of a boat); of 159 instances of gear in scenes where the gear was traveling slower than two knots, zero matched to a radar detection (Fig. S9).Generating probability rasters for matching AIS to SARMost AIS positions did not correspond to the exact time when the SAR images were taken. Hence, to determine the likelihood that a vessel broadcasting AIS corresponded to a specific SAR detection, we first developed probability rasters of where a vessel was likely to be minutes before or after a GPS position was recorded (Figs. S1,S2). We mined one year of global AIS data, including roughly 10 billion GPS positions, and computed these rasters for six different vessel classes (trawlers, purse seines, tug, cargo or tanker, drifting longlines, and others) and considered six different speeds (1, 3, 5, 7, 9, and 12.5 knots) and 36 time intervals (− 448, − 320, − 224, − 160, − 112, − 80, − 56, − 40, − 28, − 20, − 14, − 10, − 7, − 5, − 3.5, − 2.5, − 1.5, − 0.5, 0.5, 1.5, 2.5, 3.5, 5, 7, 10, 14, 20, 28, 40, 56, 80, 112, 160, 224, 320, and 448 min).For example, we queried a year of AIS data to find every example of where a tugboat had two positions that were 10 min apart from one another when the vessel had been traveling at 10 knots at the first position. We then recorded each of these locations relative to the location the vessel would have been if it traveled in a straight line, with x coordinates being in the direction of travel and the y coordinates being perpendicular to the direction of travel. When collected for hundreds of thousands of examples across the AIS dataset, the result is a heatmap of where tug boats are located 10 min after a position when it was traveling at 10 knots. The raster is centered on a point that is the extrapolated position of the vessel based on its speed. For instance, the purse seine raster that corresponds to a vessel traveling between 6 and 8 knots between 96 and 128 min after the most recent position is centered at a point that is 13.1 km (7 knots × 112 min) straight ahead of the direction the vessel was traveling. Figure S1 shows samples of these rasters for different vessels.We built rasters of 1000 by 1000 pixels for each vessel class and time interval, with the area covered by the raster dependent on the time interval (longer time intervals imply longer traveled distances, covering more area). The scale of each pixel was given by:$${text{pixel}};{text{width = max(1, }}Delta {text{m) / 1000}}$$
(1)
where Δm is the time interval in minutes, and pixel width is measured in km. Thus, if the Δm is under one minute, the entire raster is one kilometer wide with each pixel one meter by one meter. If the time is 10 min, then each pixel is 10 m wide, and the entire raster is 10 km by 10 km.Since the pixel width varies between rasters, the units of the rasters are probability per km2, thus summing the area of each pixel times its value equals one. Six vessel classes with 36 time intervals for each and six speeds led to 1296 different rasters. This probability raster approach could be seen as a utilization distribution32—for each vessel class, speed and time interval—where the space is relative to the position of the individual.Combining probability rasters to produce a matching scoreFor a few vessels (~ 4%) there was only one AIS position available before or after the scene. This resulted from a long gap in the AIS data due to poor reception, a weak AIS device, or cases where the vessels disabled their AIS. For these vessels, we used the raster values for a single position. For the vast majority of vessels, however, there was a GPS position right before and after the scene, and thus two probability rasters. We used two methods to combine these probability rasters to obtain information about the most likely location:Multiply and renormalize the rastersTo multiply the rasters, we interpolated the raster values, using bilinear interpolation, to a constant grid at the highest resolution between the before and after rasters. Then, we multiplied the values at each point and renormalized the resulting raster (Fig. S2):$$p_{i} = frac{{p_{ai} cdot p_{bi} }}{{mathop sum nolimits_{k = 0}^{N} p_{ak} cdot p_{bk} cdot da}}$$
(2)
where pi is the probability in vessel density per km2 at location i, pai is the value of the raster before the image, pbi is the value of the raster after the image. The denominator is the sum of all multiplied values across the raster, scaled by the area of each cell, da.Weight and average the rasters For this method, we weighted the raster by the squared value of the probabilities of that scene. This has the effect of giving the concentrated raster a higher weight, thus weighting higher the raster that is closer in time to the image:$${w}_{a}=sum_{k=0}^{N} {p}_{ak}^{2}cdot da$$
(3)
and the weighted average at location i is:$${p}_{i}=frac{{p}_{ai}cdot {w}_{a}+{p}_{bi}cdot {w}_{b}}{{w}_{a}+{w}_{b}}$$
(4)
where wa is the weight for raster a, wb the weight for raster b (calculation analogous to wa’s in Eq. 3), pi is the probability in vessel density per km2 at location i.To determine whether we should multiply (Eq. 2) or average (Eq. 4) the probabilities, we compared the performance of these two metrics against a direct inspection of the detections. We found that at short intervals, multiplying the rasters and renormalizing often made probability values extremely small ( {d}_{d}cdot {p}_{d} + {p}_{f}$$
(5)
where ({p}_{v}) is the probability density of the vessel presence at the location of the SAR detection (the score listed above), ({p}_{d}) is the probability that the vessel is detected by SAR, ({d}_{d}) is the density of non-broadcasting vessels in the region, and ({p}_{f}) is the density of false detections in the scene. The greater ({p}_{d}), the more dark vessels there are in a scene, and the more likely it is that any given detection is a dark vessel instead of a vessel broadcasting AIS. The right-hand side of the equation ({d}_{d}cdot {p}_{d} + {p}_{f}) should roughly equal the number of detections per unit area that do not match to AIS in the region. In other words, the probability of the vessel with AIS being at that specific location and detected by SAR (left side of the equation) should be greater than the probability of a dark vessel or a false detection at that location (right side of the equation).The total number of unmatched vessels in each studied region normalized by total area covered gives a density of non-broadcasting vessels of 2.6–2.8 × 10–5 vessels km-2 (Indian Ocean) and 6.8–7.2 × 10–6 vessels km−2 (Pacific Ocean), similar to the thresholds estimated by analysts. For the most likely number of matched vessels, we use a threshold that is halfway between the higher and lower bound of the analyst (5 × 10–5 to 1 × 10–4), 2.5 × 10–5 which is also roughly equal to the theoretical estimate of the Indian Ocean.This threshold approach performed significantly better than a metric based on the distance between the SAR detection and the most likely location of the vessel, where the likely location is based on extrapolating speed and course of the position closest in time to the image (Fig. S4).Determining whether a vessel with AIS was within a sceneVessel positions from AIS are usually available before and/or after the SAR images, and sometimes it is unclear if a vessel should have been within the scene footprint at the time of the image.To estimate the probability that a vessel (with AIS) was within a scene, we used the multiplied probability raster, summing the values inside the scene boundaries. This provides an estimate of the likelihood that the vessel was within the scene footprint at the time of the image. We applied this to every vessel that had at least one AIS position within 12 h and 200 nautical miles of the scene footprint. The vast majority of vessels were either very likely inside or outside the scene footprints, with 516 vessels having a probability of  > 95% and only 16 having a probability between 5 and 95%. We filtered out all vessels that were definitely outside of the image footprint before matching.Estimating the likelihood of detecting a vessel with SARThe AIS data show that not all vessels broadcasting AIS were captured by the RADARSAT-2 images (Fig. 3a). Using the known lengths of detected vessels with AIS, we estimated the likelihood of detecting a vessel with SAR as a function of vessel length (Fig. 3a). For vessels shorter than 60 m, we approximated the detection rate as a linear function. Treating each vessel as an individual detection, we fitted the 50th percentile using quantile regression to approximate the detection rate. For vessels above 60 m, we assumed a constant detection rate as very few vessels above this length did now show up in the SAR images. Of the 46 unique vessels larger than 62 m, 42 were detected, implying a detection rate of ~ 91%. Given that it is highly likely that large vessels will be captured by medium-resolution SAR imagery, we manually reviewed these cases to confirm that they were (almost surely) inside the scene footprints at the time the images were taken.We should note that the probability of detecting a vessel in SAR also depends on the sea state, incidence angle, polarization, material of the vessel, and orientation of the vessel. We are unable, however, to measure these effects directly so we cannot explicitly model these effects.With sufficient scenes, these effects should be randomly distributed across our scenes, so they likely account for some of the variability in detectability and the inaccuracy in our length estimates from SAR.Estimating the number and length of non-broadcasting vesselsBecause SAR does not detect all vessels, and because the length as estimated by SAR can be incorrect, there are many possible distributions of actual non-broadcasting vessels that could have produced the distribution of unmatched SAR detections that we found in the scenes. To estimate the most likely such distribution, we built a model to combine the two key relationships—between vessel length and likelihood of detection, and between vessel length and the length as estimated by SAR. This model allowed us to estimate, based on the number and distribution of SAR vessels, the likely number and distribution of actual vessels present (Fig. 3c,d).We binned the likelihood of vessel detection as a function of length into 1 m intervals, yielding a vector (alpha) of length 400. We also binned into 1 m intervals the population of lengths of all detected vessels ((ell_{D})) as reported by AIS (i.e. number of vessels at each length bin), the population of expected SAR lengths ((ell_{E})), and the population of lengths of all vessels ((ell_{A}), the quantity we wish to estimate). Thus, (ell_{D}) can be expressed as the product of (alpha) and (ell_{A}):$$ell_{D} = {upalpha } odot ell_{{text{A}}}$$
(6)
where (odot) is the element-wise product. We then estimated a matrix (L_{{}}) that relates (ell_{D}) to (ell_{E}).$$ell_{E} = Lell_{D}$$
(7)
where each element (L_{ij}) represents the probability that a vessel with length in bin j would be estimated by SAR to be of length in bin i. We calculated these probabilities as lognormal probability density functions, with one distribution per column. To estimate the scale and shape parameters of these distributions, we first fitted a quantile regression using the (non-binned) lengths from AIS of detected vessels as the predictor for the lengths reported by SAR. Assuming that the predicted 1/3 and 2/3 quantiles (as shown in Fig. 3a) represent the quantiles of a lognormal distribution, allow us to calculate the shape and scale parameters. We chose a lognormal distribution because: 1) the variable of interest, length, was always greater than zero, 2) the population of lengths was skewed towards larger values, and 3) there is an explicit and relatively simple relationship between the lognormal quantiles and the shape and scale parameters that simplified the calculations.Combining Eqs. (6) and (7) provides a relation between (ell_{A}) and (ell_{E}):$$ell_{E} = {text{L}}left( {alpha odot ell_{A} } right)$$
(8)
To estimate ({mathcal{l}}_{A}) we minimized an objective function (O({mathcal{l}}_{E},{mathcal{l}}_{o})) between the vector of expected counts binned by length (({mathcal{l}}_{E})) and the vector of counts observed in SAR binned by length (({mathcal{l}}_{o})). For this objective function, we chose the sum of the Kolmogorov –Smirnov distance between length distributions and the squared difference of the total numbers of detections. The first term controls the shape of the resulting distribution while the second one controls the magnitude. Specifically:$$Oleft( {ell_{E} ,ell_{o} } right) = max left( {left| {C_{E} – C_{O} } right|} right) + left( {T_{E} – T_{O} } right)^{2}$$
(9)
where:$$T_{x} = mathop sum limits_{ } ell_{x}$$$$D_{x} = ell_{x} /T_{x}$$$$C_{x} = cumsumleft( {D_{x} } right)$$Assessing the uncertainty in the estimationTo test how accurately our approach predicts the correct number of vessels, we performed a bootstrap simulation. We computed the vector (alpha) and the matrix L from a random subset of vessels with AIS that had a high confidence ( > 95%) of appearing within the scenes. We then used our method on the SAR detections that matched the remaining vessels to predict the number of vessels they corresponded to ((ell_{text{A}})). By running 10,000 experiments we found a mean absolute percent error of + − 9% (Figs. S5 and S6). This provides a rough estimate of the uncertainty in our prediction due to the estimation process itself. We used the distribution of these samples to estimate the 90% confidence interval that we report with our estimates. We note that this uncertainty refers to the parametrization of the model and there may be other sources of error, such as the possibility that vessels without AIS have different radar properties (e.g. made out of materials with different reflectiveness), that we did not account for in our model.Catch and effort data in the overlapping area between WCPFC and IATTCWe downloaded gridded effort and catch data from the WCPFC and IATTC websites, and compared the reported number of hooks and catch from September to December of 2019 for the area between − 140 to − 150 longitude and − 5 to − 15 latitude, a bounding box that contains our study region in the Pacific and which is entirely within both the WCPFC and IATTC convention zones. We found that the reported number of hooks for Korea is three times higher for the IATTC as it is for the WCPFC (Fig. S7), and the numbers of hooks also disagree by more than 10% for most other flag states. Catch is also 2.5 times higher for IATTC than for WCPFC for Korea as well, with catch also differing by more than 10% for most other flag states. This finding suggests that the different RFMOs may not be accounting for the same vessels in the overlap region between the two RFMOs. More