Spatiotemporal Variation in Hatching Success and Nestling Sex Ratios Track Rapid Movement of a Songbird Hybrid Zone

Hybridization often occurs at the parapatric range interface between closely related species, but fitness outcomes vary: hybrid offspring exhibit diverse rates of viability and reproduction compared with their parental species. The mobile hybrid zone between two chickadee congeners (Poecile atricapillus × Poecile carolinensis) has been well studied behaviorally and genetically, but the viability of hybrids and the underlying mechanisms contributing to hybrid fitness have remained unclear. To better characterize the fitness costs of hybridization in this system, we analyzed 21 years of data from four sites, including more than 1,400 breeding attempts by the two species, to show that rates of hatching success changed substantially as the zone of hybridization moved across the landscape. Admixture-associated declines in hatching success correlated with reduced proportions of heterogametic (female) offspring, as predicted by Haldane’s rule. Our data support an underlying mechanism implicating genetic admixture of the homogametic (male) parent as the primary determinant of offspring sex ratio, via incompatibilities on the hemizygous Z chromosome. Our long-term study is the first to directly measure changes in fitness costs as a vertebrate hybrid zone moves, and it shows that changes in these costs are a way to track the distribution of a hybrid zone across the landscape.


Introduction
Hybridization-the interbreeding of species that differ in one or more heritable characters-results in a variety of fitness consequences for both parents and offspring. The severity of survival or other fitness costs are significantly influenced by the degree of genetic divergence between the interbreeding taxa (Barton and Hewitt 1989;Harrison 1990;Mallet 2007;Abbott et al. 2013). Hybrid offspring experience a range of fitness costs relative to unadmixed parents: they may be fully viable or completely inviable. During periods of isolation and evolutionary divergence, two lineages may accumulate intrinsic reproductive barriers that only partially prevent gene flow between closely related species (Hewitt 1988;Harrison 1993;Barton 2001). Tension zones occur at range interfaces where incomplete reproductive isolation is maintained through immigration and interbreeding of parental taxa, balanced by selection against admixed offspring (Key 1968).
Environmental changes can shift species distributions and subsequently the geographic location of tension zones (Buggs 2007;Taylor et al. 2014b;Billerman et al. 2016). Tension zone movement is likely if hybridizing species differ in the strength of their response to environmental change, in their competitive ability, or both (McQuillan and Rice 2015; Taylor et al. 2015;Billerman et al. 2016). When examined using long-term data sets, tension zones are useful natural laboratories for understanding how species adapt and interact (Scriber 2011;De La Torre et al. 2014;Taylor et al. 2015). As a tension zone moves over time, a site will harbor a single-species population before hybrid encroachment, followed by genetic admixture in an actively hybridizing population with heterospecific breeding pairs; finally, there can be a shift toward an unmixed population of the second species as hybridization moves away (Buggs 2007;Rheindt and Edwards 2011;Taylor et al. 2015;Aguillon and Rohwer 2022). Research that spans this shift to and from hybridization at a site is of great value because the variation in phenotype and fitness of a population can be compared before, during, and after hybridization. This permits direct comparisons between conspecific and heterospecific pairs in a single location, eliminating geography as a variable and focusing on the role of genotypic differences to overall fitness. However, data sets for addressing such issues require large sample sizes over many years to accumulate enough variance in hybrid and conspecific pair composition to test the influence of pair genotype on fitness; such data sets are rare, and they often compare time points separated by decades (Arntzen and Wallis 1991;Engebretsen et al. 2016).
Some of the survival costs incurred by hybrid offspring occur during development (Price and Bouvier 2002). In birds this can manifest as reduced hatching success, which leads to fewer viable offspring even if clutch size is unaffected by admixture (Lijtmaer et al. 2003;Carling and Zuckerberg 2011). Haldane's rule states that the heterogametic sex is more likely to experience reduced fecundity or survival compared with the homogametic sex (Haldane 1922;Wu et al. 1996), which leads to differences in the cost of hybridization between the sexes. Evidence supporting Haldane's rule derives from several sex determination systems (Laurie 1997;Presgraves and Orr 1998;Brothers and Delph 2010), including the avian ZZ/ZW system: heterogametic female birds (ZW) may exhibit lower fitness than homogametic males (ZZ; Tegelström and Gelter 1990;Carling and Brumfield 2008). Reduced introgression of maternally inherited loci further suggests that hybrid female birds suffer greater fitness disadvantages (Carling and Zuckerberg 2011;Jacobsen and Omland 2012;Gowen et al. 2014).
Black-capped chickadees (Poecile atricapillus; hereafter, BCCH) and Carolina chickadees (Poecile carolinensis; hereafter, CACH) are common, widespread North American passerines (family Paridae) that are sister taxa (Harris et al. 2014) with an estimated divergence time of~2.5 million years (Gill et al. 2005). These species hybridize at the interface of their ranges, along a narrow contact zone in the eastern United States from Kansas to New Jersey: BCCH occur in the north, and CACH occur in the south (Tanner 1952;Brewer 1961). The hybrid zone between these species exhibits the classic characteristics of a tension zone (Bronson et al. 2003;Taylor et al. 2014a). Characterization of chickadee behavior, morphology, and hybridization has occurred in multiple transects across the hybrid zone in Ohio, Pennsylvania, and Virginia (Robbins et al. 1986;Bronson et al. 2005;Reudink et al. 2007;Sattler et al. 2007;Olson et al. 2010;Davidson et al. 2013;Taylor et al. 2014a).
Preliminary evidence suggests that hybridizing chickadees experience a reduction in fitness due to low hatching success of admixed offspring (Bronson et al. 2003(Bronson et al. , 2005Huynh and Rice 2019), but these studies are potentially confounded by small sample sizes, geographic effects of the particular microclimate where measurements took place, and stress on chickadees that were artificially placed in heterospecific pairs by researchers. Hybrid chickadees exhibit cognitive impairment in learning and memory Rice and McQuillan 2018), but the relationship between these deficiencies and fitness is not fully understood.
We used nesting records and genetic samples collected over 21 breeding seasons (1998-2018) across a transect in eastern Pennsylvania to investigate the spatiotemporal consistency of signatures of selection in the chickadee hybrid zone. We examined northward shifts in the geographic location of active hybridization and the relationship between hybridization frequency and changes in hatching success and nestling sex ratios. We hypothesized that increased admixture between the species in the region of active hybridization would negatively impact hatching success owing to intrinsic genetic incompatibilities between paired adults. We therefore examined (1) whether hatching success, assessed at the whole-population level, varied through time across our sampling sites and (2) whether variation in degree of admixture of parent birds affected hatching within a hybrid zone population. We further hypothesized that female hybrids should experience stronger intrinsic incompatibilities (Haldane's rule) such that nestling sex ratios should be increasingly male biased in populations experiencing elevated admixture.

Field Methods
We monitored breeding populations at four sites ( fig. 1; see also maps in Reudink et al. 2007;Taylor et al. 2014b). Previous genomic analysis (Taylor et al. 2014b) detected differences in population genetic structure among these sites between 2002 and 2012 as follows. Only CACH occupied the southernmost site, Great Marsh (GM). Nolde Forest (NF) hosted a mixed population in 2002, but CACH predominated there by 2012. Hawk Mountain (HM) comprised a primarily BCCH population when sampled in 2002, but by 2012 both species were present and hybridizing. As of 2012 the northernmost site, Tuscarora State Park (TU), supported an unadmixed BCCH population.
We obtained data from breeding chickadees using artificial nest snags (based on Grubb and Bronson 1995) and a few natural cavities. Data reported here spanned 1998-2018 at GM and NF, 2000-2018 at HM, and

Fitness Costs Track a Moving Hybrid Zone 265
2006-2018 at TU. Whenever possible, we recorded the date of the first egg laid in each nest (clutch initiation date), total clutch size, and the proportion of eggs that successfully hatched (hatching success). We included hatching success data only for clutches not subject to predation, snag destruction (from weather, treefall, bears), or parental abandonment that occurred before the end of the typical 12-13-day incubation period.
We captured adult birds using mist nets at active nests and sampled nestlings before fledging. We obtained blood (30-40 mL) by puncturing the brachial vein and then preserved the blood samples in 300 mL of lysis buffer (100 mM Tris, pH 8, 100 mM NaCl, 0.5% sodium dodecyl sulfate buffer, double-distilled H 2 O) and stored them at 47C.

Molecular Methods
We extracted DNA using a DNeasy blood and tissue kit or Puregene core kit A (Qiagen, Germantown, MD). We determined each nestling's sex via DNA amplification of the chromohelicase DNA-binding gene (CHD) on the sex chromosomes using the forward primer 1237L and the reverse primer P2 (Griffiths et al. 1998;Kahn et al. 1998). To characterize species identity of parent chickadees, we performed restriction fragment analysis on eight speciesdiagnostic single-nucleotide polymorphism

Data Analyses
To determine sex ratio within each brood, we calculated the proportion of offspring surviving to an age permitting blood sampling (≥8 days after hatch) that were male. A small percentage of eggs or corresponding nestlings (!3% of all eggs) disappeared before we could sex the brood; our analyses omitted these. We also lacked data about the sex of unhatched eggs. Visual inspection of the contents of a subset of unhatched eggs during the 2014-2016 breeding seasons at HM (n 1 100) yielded no evidence of developing embryonic tissue from which we could extract DNA for sexing. We scored hybrid identity by determining the ancestral genotype of each species-diagnostic locus. To determine the proportion of alleles from each parental species, we divided the number of CACH alleles by the total number of alleles assayed for each bird. We sought to compare hybrid individuals (F 1 or equivalent; 0.50 CACH) with unadmixed (i.e., parental) birds (either 1.0 CACH or 0.0 CACH). The absolute difference of a bird's proportion of CACH alleles subtracted from the value of an F 1 hybrid (0.50 CACH) yielded a scale from the most hybrid individuals (0.0) to the least admixed species individuals (0.50; Bronson et al. 2003).
For breeding pairs with genotype scores available for both parents, we calculated a pair genotypic index (GI) using the following formula: pair GI p proportion CACH alleles male 1 proportion CACH alleles female 2 2 0:5 : We calculated a compatibility index (CI; Bronson et al. 2003) for each breeding pair by measuring the proportion of heterozygous alleles at all eight loci. We used the formula where k represents the total number of loci (in our case, eight) and d represents the measure of homozygosity between parents, with d p 1 for identical homozygotic parents at a given locus, 0.5 for heterozygotic pairs, and 0 for opposite homozygotes. Genotyping focused on breeding adults at HM because this population experienced more change in hybridization frequency and genetic admixture throughout the study than did our other populations (Taylor et al. 2014a). Genotyping of breeding adults from GM, NF, and TU did not provide new information compared with previous assessments of genetic ancestry at these sites (Taylor et al. 2014a; see above).

Statistical Analysis
We generated all statistical analyses and figures in R (ver. 3.5.1; R Core Team 2019) using the packages lme4 (Bates et al. 2015), blme (Chung et al. 2013), ggplot2 (Wickham 2009), and optimx (Nash and Varadhan 2011). We used generalized linear mixed effects models that assumed binomial errors and that specified the identity of the individual females and males as crossed random effects to account for repeated pairings between the same pairs as well as new pairings over time. For analyses of GI and CI, we weighted indexes by the number of loci from each sex that were included to adjust for any discrepancy in the information used to calculate each index. Data for GI and CI were limited to only HM, and we did not have sufficient resolution to reliably estimate nonzero random effects for both male and female breeder identities. Therefore, for analyses of end points, we used Bayesian generalized mixed models, using the package blme (Chung et al. 2013), that specify a weak Wishart prior distribution for random effects, allowing more stable nonzero estimates of random effects' standard deviations and associated confidence intervals. We based inferences for all of the analyses on likelihood ratio tests (and confirmed using corrected Akaike information criterion to ensure consistency) comparing our most inclusive models (i.e., a model containing all focal fixed and random effects and their interactions) with simplified models that excluded the focal parameter. We interrogated model fits by examining profile zeta and density plots (Bates et al. 2015), and where appropriate (i.e., strong evidence of well-behaved quadratic likelihood profiles), we report Wald z statistics and modelgenerated averages with Wald-approximated 95% confidence intervals for comparisons of fixed effects (e.g., location) unless stated otherwise. We used the prediction function from Bolker (2015) to generate and plot model predictions and confidence intervals for continuous covariates in our models (e.g., time and GI). For comparison between sites and time frames, we used Welch's two-sample t-tests.

Sample Sizes
The total number of nests used in our analyses varied among the four study sites and variables measured (table 1). Analyses that follow focus especially on HM because of its central position relative to the hybrid zone (Taylor et al. 2014a).
Hatching success totals depended on years of monitoring and number of active nests. We had the highest number of broods with complete sex ratios at HM that experienced hybridization during part of the 19-year sampling period.

Spatiotemporal Variation in Hatching Success
Hatching success varied over time differently among the four sites in ways that reflected regional movement of the hybrid zone ( Combining data from all sites and years, hatching success was independent of clutch size (fig. S1; x 2 p 0:004, P p :95, n p 1,189 clutches). Clutch initiation date also did not influence hatching success ( fig. S2; x 2 p 0:82, P p :37, n p 895 clutches).

Ancestry Patterns in the Mixed Hawk Mountain Population
We

Variation in Nestling Sex Ratios
Across all years and all four sites, nestling sex ratios were slightly male biased (0.53% male, n p 373 broods, 95% confidence interval p 0:51-0:56). Across three sites for which we collected data on both hatching success and nestling sex ratio (GM, NF, and HM), nestling sex ratios did not show a relationship with decreased hatching success ( fig. S5; hatching success#nestling sex ratio interaction: x 2 p 0:02, P p :90, n p 174 nests). At HM, increasingly hybrid GIs of pairs were associated with increasingly male-biased nestling sex ratios ( fig. 4a; x 2 p 6:05, P p :01, n p 139 pairs). In contrast, the pair CI of breeding pairs was not clearly associated with variation in nestling sex ratio (fig. S6; x 2 p 2:30, P p :13, n p 130). The bias in nestling sex ratio toward males increased with GI of the male parent, such that the most hybrid males tended to produce proportionally the most sons ( fig. 4b; x 2 p 9:00, P p :003, n p 157). GI of the female parent alone did not affect nestling sex ratio ( fig. 4c; x 2 p 0:19, P p :67, n p 152).

Discussion
Within a transect of a chickadee hybrid zone in eastern Pennsylvania, we show how hybrid fitness, as measured by decreased hatching success, mirrored changes in hybrid distribution as hatching success continually declined at one study site (HM) and as high hatching success was regained at another site (NF) across 19 breeding seasons. We measured how increasing hybridization over time in one population (HM) was associated with decreased hatching success and increasingly male-biased nestling sex ratios. These findings constitute evidence that hybrid fitness costs, including reduced hatching success and male-biased nestling sex ratios, increased with breeding pair admixture. Our patterns of male-biased sex ratios are also consistent with expectations of Haldane's rule, which predicts that genetic incompatibilities in the heterogametic sex (in this case, females) will cause sex-biased survival in hybrid offspring.
Hatching success decreased at HM as the hybrid indexes of parental breeding birds increased. Our study presents the first evidence of long-term changes in the distribution of hybrid fitness in a natural avian hybrid zone, providing the first temporal component to previously reported associations between a decline in hatching success and increasingly admixed parental genotype using short-term or single-site observations (Lijtmaer et al. 2003;Carling and Zuckerberg 2011;Huynh and Rice 2019). The HM population had a wide range of average hatching success across different breeding seasons that was correlated with genetic admixture in the population. This suggests that the reduced hatching success was a result of genetic incompatibilities between chickadee species, independent of local geography. The genetically regulated, as opposed to environmentally mediated, reduction in hatching success provides a potential mechanism to explain the consistent genomic signatures of reproductive isolation reported by Wagner et al. (2020). That study identified genes with reduced introgression across the same Pennsylvanian transect of the chickadee hybrid zone and found enrichment for gene ontology (GO) categories that may contribute to dysfunction in the hybrids. Our results provide a specific developmental timing of selection for some of these gene functions. One identified GO category, calcium ion release, is required for skeletal development, and in the embryonic development of birds, calcium is additionally absorbed through the eggshell (Simkiss 1967). Transfer RNA methylases, another GO category identified by Wagner et al. (2020), are upregulated in the embryonic tissues of mammals during differentiation (Kerr 1971 Figure 4: Variation in the proportion of male nestlings in relation to the genotypic index of parent(s) at Hawk Mountain. The proportion of male offspring decreased with increasing unadmixture of hybrid breeding pairs (a) and increasingly unadmixed male parents (b), conforming to Haldane's rule, but did not vary with species identity of the female parent (c). Lines denote model predictions with 95% confidence intervals, data points represent individual pairs or breeders, and larger points denote multiple coincident values. d, Heterogametic (female) offspring inherit only their hemizygous Z chromosome from the male parent. Incompatibilities conferred from the single copy of Z chromosome alleles are therefore inherited only from the homogametic (male) parent; the comparatively small W chromosome may not contain as many potentially detrimental loci. CACH p Carolina chickadees; BCCH p black-capped chickadees. Original artwork by Kristen Orr. the developmental timing of substantial selective pressures in chickadee hybridization and provide a potential life stage at which the identified GO categories impact fitness. Reduced hatching success, combined with the reduced spatial memory and problem-solving of hybrid chickadees (Mc-Quillan et al. 2018) and potential difficulty of hybrids surviving cold winter temperatures (Olson et al. 2010;Taylor et al. 2014b), presents a more complete view of the reduced fitness in hybrids that occurs at all life stages.
Using nestling sex ratio as a measurement, we present the most comprehensive counts of offspring sex supporting Haldane's rule for any wild bird system. The chickadee hybrid zone follows Haldane's rule, with increasingly male hybrid offspring associated with parents with an increasingly hybrid pair GI. Haldane's rule is typically illustrated by genetic evidence such as reduced cline widths or increased differentiation of sex chromosome and mitochondrial loci (Jacobsen and Omland 2012;Gowen et al. 2014). Haldane's rule specifically states, "when in the F 1 offspring of two different animal races one sex is absent, rare, or sterile, that sex is the heterozygous sex" (Haldane 1922, p. 7). Our counts of male and female offspring are a way to measure Haldane's rule that is not inferred by genetic clines or differentiation. Our data are consistent with the previously reported decrease in Z chromosome cline widths across the same transect (Taylor et al. 2014b;Wagner et al. 2020).
Our results are consistent with a hypothesized mechanism underlying Haldane's rule. In the homogametic sex, if genetic incompatibility occurs between a locus on one sex chromosome and a locus on an autosome, the allele on the homologous sex chromosome deriving from the second parental species is potentially available to compensate. However, a second copy does not exist in the heterogametic sex. The heterogametic sex, with one Z chromosome (X in mammals) and one W chromosome (Y in mammals), uniformly inherits its Z chromosome from the homogametic parent ( fig. 4d). Therefore, the heterogametic sex does not pass its Z chromosome on to heterogametic offspring. Admixture of the gametes in the homogametic parent largely mediates the increased susceptibility to genetic incompatibility in heterogametic offspring. We show that increase in male sex bias of nestlings is not correlated with the hybrid index of the female parent but correlates instead with the hybrid index of the male parent ( fig. 4). This provides evidence that genetic incompatibility through autosome and sex chromosome interactions could be a contributing factor in this manifestation of Haldane's rule. Importantly, we were unable to confirm whether unhatched eggs were fertilized. Therefore, we cannot rule out the additional possibility of prezygotic barriers that prevent sperm with admixed Z chromosomes from fertilizing eggs. However, in other nonhybridizing parids, the vast majority of eggs that fail to hatch do show evidence of fertilization: only 2.2% and 1.3% of unhatched eggs were unfertilized in blue tits (Cyanistes caeruleus; n p 419) and great tits (Parus major; n p 375), respectively, when evaluated using methods permitting detection of embryos that died before a visible zygote was discernible (Hemmings and Evans 2020). Our current work using the same methods aims to determine whether fertilization failure is higher in hybrid zone chickadee populations and to reveal the sex ratio among embryos that die before hatching.
Although we did not account for differences in physiological or social condition of adults or for quality of the local environment, it is possible that the conditions of hybrid adults are not the same as those of unadmixed individuals and that these conditions influence hormonal levels in the female parent. The stress hormone corticosterone appears to play a role in proximate control of sex ratio manipulation (Pike andPetrie 2005, 2006). Corticosterone may increase in response to stressors that include low territory quality or low mate quality (Pike andPetrie 2005, 2006), which can subsequently cause a bias toward more female offspring. This bias may be the result of a regulatory role of corticosterone in chromosomal segregation during the first meiotic division (Pike and Petrie 2003). In parids, as in many passerines, local resource competition may be greater among philopatric sons than among dispersing daughters, potentially making it adaptive for mothers to bias nestling sex ratios toward females (Andreu and Barba 2006;Song et al. 2016). Because our results show a sex bias in favor of males, and not females, it is likely that sex ratios among nestlings of hybrid fathers stem from genetic incompatibilities following Haldane's rule (Pike andPetrie 2005, 2006) and not local resource competition.
Our study compares parental genotype, brood hatching success, and sex ratio, but we did not fully consider complicating effects of extrapair paternity. Extrapair copulations, which may maximize fitness by allowing birds to produce offspring with breeding partners other than their social mate, are widespread in birds (Ferretti et al. 2016;Bolton et al. 2017;Brouwer and Griffith 2019). Breeders in unadmixed BCCH populations engage in extrapair copulations (Otter et al. 1994;Mennill et al. 2003), and extrapair paternity accounted for 26% of nestlings from 2000 to 2003 in the hybridizing NF population (Reudink et al. 2006). If female chickadees perceive heterospecific or admixed social mates as suboptimal, they could respond by engaging in extrapair copulations with conspecific males (Veen et al. 2001). However, in our chickadee populations, the prevalence of extrapair offspring was not tied to genetic distance between mated parents, and females did not preferentially produce extrapair offspring with conspecific males (Reudink et al. 2006). Additionally, extrapair Fitness Costs Track a Moving Hybrid Zone 271 copulations in our system would obscure the relationships that we detected between social parent genotype, hatching success, and nestling sex ratio. Therefore, correlations between these measurements would likely have been even stronger had we analyzed parentage exhaustively: in our analyses, some extrapair nestlings in broods associated with hybrid/mixed pairs may have hatched because they were not hybrids, whereas some fertilized eggs in clutches associated with conspecific social pairs may have died because they were hybrid products of extrapair mating. We still detect strong correlations of admixed genotype and decreased hatching success despite the possible impact of extrapair copulations.

Conclusions
We used a data set spanning 21 breeding seasons and four localities to document how genotypic variation has impacted the fitness of hybrid offspring between sister chickadee species. We showed that hatching success decreased as the genetic admixture of parents in the pair increased, for the first time comprehensively characterizing the fitness variability that dovetails into well-studied behavioral differences in this hybrid zone. We suggest that genetic incompatibilities during early embryonic development, or possibly during fertilization, cause the decreased hatching success. Nestling sex ratio data were consistent with the expectations of Haldane's rule as a mechanism for hybrid sex ratio biases but was related to the admixture proportion of the male and not the female parent. This result implies that the hybrid fraction of the paternally inherited Z chromosome could be responsible for observed incompatibilities and is a compelling example of parent-of-origin effects on genetic incompatibilities in natural hybrid zones. We showed that fitness costs can be used to track hybrid movement at both the trailing and the leading edges of a hybrid zone. As opposed to genotypic cline widths that indirectly approximate historical movements, our long-term data set uniquely used direct measurements of fitness costs (hatching success) to measure the movement of the hybrid zone across the landscape. These fitness costs represent a new way to track hybrid zone movement in real time and could be implemented for monitoring hybrid distribution and dynamics in future studies of natural hybrid zones. and Pickering Creeks Conservation Trust, the Nature Conservancy, and the Pennsylvania Department of Conservation and Natural Resources for site access. Keith Bildstein, Laurie Goodrich, Todd Baumann, and HMS staff provided generous logistical assistance. We thank many Villanova students-especially Sue Guers, Matt Reudink, Lindsay Zemba Leigh, Steve Van Pelt, Breanna Bennett, and Sarah Polekoff-for contributions to data collection and laboratory work; Keith Danielson for laboratory guidance and editorial input; Scott Taylor for detailed manuscript feedback; and Sue McRae, Sara DeLeon, and Peri Bolton for additional editorial assistance. We thank Kristen Orr for providing the chickadee illustration, Lu Xu for Chinese abstract translation, and three anonymous reviewers and editorial staff for their substantial input on the manuscript. For funding support, we thank Villanova University and the many societies providing student grants (to R.J.D., E.S.B., and other participants), including the Wilson Ornithological Society, the American Ornithological Society, the Association of Field Ornithologists, and Sigma Xi. This is Hawk Mountain conservation science contribution no. 358.

Data and Code Availability
All data sets and R code required to replicate analyses have been deposited and are available in the Dryad Digital Repository (https://doi.org/10.5061/dryad.c866t1g8d; Driver et al. 2022).