# Avian extra-pair paternity in the last European primeval forest

Joanna Sudyka, Irene Di Lecce, Grzegorz Hebda, Patryk Rowiński, Robert Rutkowski, Anna Szczuka, Charles Perrier, Marta Cholewa, Marta Maziarz, Tomasz Wesolowski, et al.

► To cite this version:

Joanna Sudyka, Irene Di Lecce, Grzegorz Hebda, Patryk Rowiński, Robert Rutkowski, et al. Avian extra-pair paternity in the last European primeval forest. Functional Ecology, inPress, 10.1111/1365-2435.70223. hal-05484411

HAL Id: hal-05484411

https://hal.inrae.fr/hal-05484411v1

Submitted on 29 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Avian extra-pair paternity in the last European primeval forest

Joanna Sudyka<fnref n="1" />,² ID | Irene Di Lecce³ | Grzegorz Hebda⁴ | Patryk Rowiński⁵ |
Robert Rutkowski⁶ | Anna Szczuka⁷ | Charles Perrier⁸ |
Marta Cholewa⁹ | Marta Maziarz⁶ | Tomasz Wesołowski⁹,† | Marta Szulkin³

**Correspondence**

Joanna Sudyka

Email: joanna.sudyka@uj.edu.pl

**Funding information**

Ministerstwo Edukacji i Nauki, Grant/

Award Number: 2PO4C 075 27;

Narodowe Centrum Nauki, Grant/Award

Number: 2016/21/B/NZ8/03082

Handling Editor: Caroline Isaksson

**Abstract**

1. Primeval forests offer a reference baseline to understand the origins and evolution of mating systems, as their relatively undisturbed environment provides a glimpse into how ecological interactions and natural selection play out in their original context.

2. We established rates of extra-pair paternity (EPP) for two bird species referential to evolutionary ecology: blue tits and great tits, breeding in natural cavities in Białowieża National Park, the sole remnant of European lowland primeval forests.

3. Genotype-by-sequencing of 889 blue tits and 1256 great tits revealed that 48% of blue tit broods and 39% great tit broods were of mixed paternity, with 15% blue tit and 14% great tit nestlings not sired by their social father. These referential levels align with median values of EPP reported in nestbox studies in secondary and managed forests, suggesting that certain reproductive strategies are advantageous across environments. Observed EPP levels did not affect reproductive success or parental investment, indicating no or limited adaptive value for EPP in the context of a primeval forest ecosystem.

4. Our study challenges the assumption that EPP is maintained by direct or indirect fitness benefits. Rather than indicating recent adaptive responses, the observed patterns may be consistent with long-standing and resilient mating strategies that have persisted under increasing global ecological pressures.

**KEYWORDS**

avian mating systems, Białowieża, evolutionary ecology, extra-pair paternity, natural selection, primeval ecosystem, promiscuity

<footnote n="1">¹Institute of Environmental Sciences, Jagiellonian University, Kraków, Poland; ²Groningen Institute for Evolutionary Life Sciences (GELIFES), Groningen, Netherlands; ³Institute of Evolutionary Biology, Faculty of Biology, University of Warsaw, Warsaw, Poland; ⁴Institute of Biology, Opole University, Opole, Poland; ⁵Department of Forest Zoology and Wildlife Management, Institute of Forest Science, Warsaw University of Life Sciences SGGW, Warsaw, Poland; ⁶Museum and Institute of Zoology, Polish Academy of Sciences, Warsaw, Poland; ⁷Laboratory of Ethology, Nencki Institute of Experimental Biology, Polish Academy of Sciences, Warsaw, Poland; ⁸CBGP, INRAe, CIRAD, IRD, Montpellier SupAgro, University of Montpellier, Montpellier, France and ⁹Laboratory of Forest Biology, Wrocław University, Wrocław, Poland</footnote>

<footnote n="1">¹Deceased</footnote>
# 1 | INTRODUCTION

Primeval forests, characterized by long-term ecological continuity free from major direct human disturbance, are a rarity across the globe. In Europe, these forests are even scarcer and contribute only 0.4% to the global area of primary forests (FAO, 2020). Białowieża Forest is the sole iconic remnant of a mixed forest that once stretched across Central Europe and is a UNESCO World Heritage Site (UNESCO World Heritage Centre, n.d.). The area's unique history allowed its protection through the ages, long before large-scale anthropopopulation (Marris, 2008). The role of such ecosystems is irreplaceable, showing resilience to non-native species invasion (Mo et al., 2023) and contributing to global climate and biodiversity targets (Mo et al., 2023). Ecological networks inherent to Białowieża Forest are our best reference baseline to a specific biome of origin for many native European plant and animal species. The fundamental importance of intact landscapes is much discussed in biodiversity and community ecology (Watson et al., 2018). In contrast, mating systems and sexual selection, the driving forces of adaptation and evolution, remain understudied in these ecosystems. However, they hold a critical role in establishing evolutionary baselines for understanding how such processes may shift under global environmental change. Much of what we know about mating strategies comes from studies carried out in more accessible or impacted ecosystems, such as managed forests and fragmented habitats, or in controlled environments such as enclosures or laboratories. These strategies are deeply interwined with ecological interactions such as predator-prey dynamics and resource/breeding site competition. Because these ecological interactions can be disrupted in modified habitats, fragmentation, and anthropogenic disturbance can indeed alter mating behaviour, a phenomenon well documented in avian species (Di Lecce et al., 2024; Kempenaers et al., 2010).

Promiscuous mating behaviour results in extra-pair paternity (EPP), defined as cases where offspring are sired by males other than the female's social mate. EPP can benefit males, who increase their contribution to the gene pool, and females, who secure optimal genes for their offspring (Westneat & Stewart, 2003). These functional insights are mainly derived from avian populations breeding in secondary and managed forests, thereby experiencing altered breeding ecology compared with natural populations that evolved and optimized their reproduction in natural cavities of primeval forests. It remains unresolved if EPP rates among socially monogamous birds reported to date, where on average 11% of offspring have an extra-pair father (Griffith et al., 2002), represent the natural benchmark for adaptation and evolution of mating strategies. Moreover, nestbox provisioning in itself is another human impact on ecosystems, as it may increase breeding densities, alter predation pressures, and modify reproductive success (Di Lecce et al., 2023; Norris et al., 2018; Purcell et al., 1997; Sudyka et al., 2022; Wesołowski, 2011). EPP may increase with the density of breeding pairs (Mayer & Pasinelli, 2013), and the opportunity to engage in extra-pair copulations is greater in

anthropogenically altered habitats (Di Lecce et al., 2024), leading us to expect lower EPP levels in unaltered primeval habitats.

In this study, we examined EPP of two model species for avian mating behaviour (Brouwer & Griffith, 2019) in natural cavities of Białowieża primeval forest (Figure S1). Specifically, we genotyped 889 blue tits (*Cyanistes caeruleus*) at 14,019 SNP markers and 1256 great tits (*Parus major*) at 13,754 SNP markers with the aim of establishing reference baseline levels for EPP and assessing the applicability of existing nestbox-based ecological models that quantify EPP levels in natural populations. We further empirically tested how phenological and environmental parameters impacted EPP levels derived from natural cavities to uncover the ecological and evolutionary drivers of mating strategies under primeval conditions. Finally, we assessed whether EPP covaried with reproductive success and parental investment to establish its broad adaptive value and understand natural selection in this complex ecosystem. We used the number of nestlings that successfully fledged as a proxy for reproductive success (Weatherhead & Dufour, 2000) and nestling fledging age as a proxy for parental investment. Juveniles that stay longer in the nest require more parental resources, which enhances survival through increased parental care (Martin et al., 2018) but also increases predation risk (Keedwell, 2003), especially in primeval forests with abundant predators (Wesołowski et al., 2015).

# 2 | MATERIALS AND METHODS

## 2.1 | Ethics approval

All methods and sampling detailed in this manuscript were performed under the approval of the Local Ethical Committee under Warsaw University of Life Sciences, SGGW, no. 50/2002, and the Directorate of Białowieża National Park dated 19/04/2005.

## 2.2 | Study site, fieldwork, and sampling

Białowieża Forest is located in the heart of the European plain along the Polish-Belarusian border (Białowieża village: 52°41′N, 23°52′E). It is the last remaining fragment of the vast lowland forests that once covered much of temperate Europe, distinguished by its remarkable size and exceptional state of preservation (Wesołowski et al., 2002, 2010, 2015). Most of its 613 km² Polish section (about 45% of the total area) is managed, yet, a 48 km² (4800 ha) block of primeval old-growth forest within Białowieża National Park (BNP) remains strictly protected. This pristine forest is a multi-story, mixed-species ecosystem with uneven-aged stands featuring numerous veteran trees, including the Norway spruce (*Picea abies*), up to 52 m high (Wesołowski & Rowiński, 2014). Unmanaged deadwood accounts for 20%-25% of the total wood volume, fostering exceptional biodiversity (Wesołowski et al., 2015). Fieldwork was conducted in four large plots (33-54 ha each) within BNP, totalling 185 ha and spaced
1-2 km apart (Figure S1). Three plots (C, M, W) were in oak-lime-hornbeam stands dominated by hornbeam (Carpinus betulus), small-leaved lime (Tilia cordata), pedunculate oak (Quercus robur), spruce (Picea abies), and continental maple (Acer platanoides) (Maziarz et al., 2016; Wesołowski & Rowiński, 2012). The average breeding pair density per 10 ha in the oak-hornbeam-lime forest patches in BNP was 4.0 pairs for blue tits and 4.9 pairs for great tits (Wesołowski et al., 2010). The fourth plot (K), with similar densities, was in a swampy riverine forest of alder (Alnus glutinosa), ash (Fraxinus excelsior), and spruce, and birds breeding in plot K had access to nearby drier hornbeam-covered areas (Wesołowski et al., 2002). Natural cavities were abundant in all plots, and tits exclusively used these for breeding (Wesołowski et al., 2015). Intensive nest searches were conducted to locate all natural cavities in two of the oak-lime-hornbeam areas: plot C (48 ha) and plot M (54 ha) (Maziarz et al., 2016; Wesołowski & Rowiński, 2012). Additional nest searches were conducted less intensively during general mapping censuses in plots W (not to interfere with a concurrent flycatcher research; area of 50 ha) and K (because of its different habitat; 33 ha; Wesołowski & Rowiński, 2014).

Blue tit DNA samples were collected between 2005 and 2008 (Wesołowski & Rowiński, 2012), and great tit DNA samples were collected between 2007 and 2011 (Maziarz et al., 2016; Table S1). Even before breeding began, birds were tracked and caught from winter (November, December) till early spring (beginning of April). Birds were uniquely colour-ringed, aged as first-year breeders or older based on wing moult patterns, and their movements were recorded on field maps to estimate the number and distribution of breeding pairs, simplifying mapping and nest discovery. Natural cavities used for breeding were found by following birds mostly during nest-building and were regularly monitored to establish brood productivity and nest fate (failed or successful), with visit frequency adjusted based on the nesting stage (Maziarz et al., 2016; Wesołowski & Rowiński, 2012). Lower cavities (up to 6 m) were inspected from the ground or with a ladder, while higher ones by climbing with special tree-climbing spikes. This allowed us to acquire data on lay date, clutch size, hatch date, number of hatchlings, fledging date, number of successfully fledged nestlings, and fledging age of nestlings (i.e., duration of the nestling period, from hatching till fledging). An array of environmental variables that precisely described each cavity was also collected, such as cavity height, entrance and nest chamber dimensions, and cavity exposure (see Supporting Information: 'Phenological and environmental variables' and Figures S2 and S3 for a complete list of variables and their correlations). The effect of these phenological and environmental variables was further tested against EPP rates. The number of fledged nestlings and fledging age were treated as the best proxies for reproductive success and parental investment obtained in our study system. Both variables were tested for covariation with brood EPP levels. Twelve days after hatching, two tail feathers were collected and immediately stored in 96% ethanol from (i) nestlings safely removed from nests using a custom-made loop tool and (ii) from parents captured using mist nets or traps near nest entrances (if they had not been captured before, i.e., they

were unringed). All nestlings were ringed with a uniquely numbered ornithological ring.

## 2.3 | Genetic analysis and parentage assignment

Blue tit feather DNA was extracted with QIAamp DNA Investigator Kit (adults; optimized for keratinized cells) and QuickGene DNA tissue kit (nestlings; containing blood and tissue particles in the calamus; Qiagen, Hilden, Germany). Great tit samples were extracted with Quick-DNA™ Miniprep Plus Kit (Zymo Research, optimal for both adults and nestlings) according to the manufacturer's protocol. A section of the feather interior (shaft's superior umbilicus), where the blood clot is located, was cut into small pieces (2 mm each) for overnight lysis at 38°C. DNA sequencing was performed by Diversity Arrays Technology Pty, Ltd. (Canberra, AU) using DArTseqLD, a high-throughput genotyping by sequencing method that employs genomic complexity reduction using restriction enzyme pairs (Kilian et al., 2012). In blue tits, we obtained 908 samples genotyped at 14,019 SNP markers and in great tits, 1286 samples genotyped at 13,754 SNP markers. Details on the DArT sequencing technology can be found in Di Lecce et al. (2023) and Di Lecce et al. (2024). All subsequent analyses were run separately for each species in R (v. 4.3.2) (R Core Team, 2023). Using dartR (v. 2.9.7) (Gruber et al., 2018), we filtered out individuals with call rates below 20% and loci with call rates below 40% in blue tits, and both individuals and loci with call rates below 40% in great tits. These thresholds were selected to retain as much information as possible, given the low call rate in adults (Figure S4). After filtering, 12,248 SNP markers and 884 individuals were kept in blue tits (889 individuals when including five recruits, whose genotypes were used twice: first to determine their EPO status within their brood of rearing, and later as parental genotypes in subsequent breeding seasons), and 12,191 SNP markers and 1256 individuals were kept in great tits. We estimated genetic relationships among individuals using the function snpgdsGRM with the method GCTA (Yang et al., 2011) implemented in SNPRelate (version 1.36.0) (Zheng et al., 2012). snpgdsGRM was run with a minor allele frequency threshold of 0.005 and a missing rate of 0.05 in both species. A total of 2136 loci were used in blue tits and a total of 3711 loci were used in great tits. We confirmed that the effect of using these low filtering thresholds on genetic relatedness was negligible, by testing it on a subset of individuals with high call rates: relatedness was calculated using 40% of the SNP data and compared with relatedness calculated with 100% of the SNP data (Figure S5). We compared the resulting genomewide relatedness matrix (GRM) with a social pedigree of all individuals ringed in the field created using *ggroups* (v. 2.1.0) (Nilforooshan & Saavedra-Jiménez, 2020). Following Perrier et al. (2018) and Di Lecce et al. (2024), discrepancies between the GRM and the social pedigree were identified separately among (i) mothers with their offspring, (ii) social fathers with their offspring and (iii) offspring within their own nest, in nests where the social father was not sampled. Father-offspring pairs (social relatedness=0.5) showing GRM relatedness estimates below
0.3 (based on relationships between mothers and their offspring; Figure S6a) were classified as instances of EPP. In nests where the social father or both social parents were not sampled (e.g., when catching attempts failed or when a parent abandoned the nest—37% of nests in blue tits and 7% in great tits), we did not assign within-pair or extra-pair status to individual nestlings but only recorded whether the brood was of mixed paternity, that is, whether EPP occurred (these nests were included in all analyses on EPP occurrence, but not on the proportion of extra-pair offspring—EPO). Therefore, if there were pairs of siblings within a given nest (social relatedness=0.5) with GRM estimates between 0.15 and 0.35, the nest was classified as containing half-siblings; if all pairs of siblings within a given nest had GRM estimates above 0.35, the nest was classified as containing full-siblings. Nestlings with GRM relatedness estimates below 0.15 with both the social parents and the social siblings (social relatedness=0.5) were classified as instances of brood parasitism (Figure S6). The SNP parentage assignment was additionally validated with a microsatellite parentage analysis in blue tits, and the results were highly consistent (Supporting Information: Microsatellite parentage analysis and comparison with SNP data). Overall, the occurrence of mixed paternity in a brood was highly consistent between the methods (McNemar's test comparing EPP detected by microsatellites vs. SNPs: $\chi^2=0, df=1, p=1$), and the relationship between the number of offspring assigned as extra-pair per brood was strong between the methods (Spearman's $\rho=0.816, p<0.0001$).

## 2.4 | Statistical analyses

To examine the role of environmental and phenological breeding parameters on EPP, we fitted EPP occurrence (0 if none of the nestlings was extra-pair or 1 if at least one nestling was extra-pair in a nest) and EPO proportion (the number of EPO over the number of genotyped nestlings per brood, the ratio coded as a continuous value between 0 and 1) as dependent variables in two separate generalized linear mixed models (GLMMs) for each species separately in R (v. 4.2.2) (R Core Team, 2022). As there was mortality recorded between hatching and the time of sampling, and parents removed dead nestlings early on during development, we included all broods in which we successfully sampled and genotyped at least 50% of the original number of hatched nestlings, to gain a representative sample size of the original brood in all analyses involving EPP occurrence and EPO proportion (8 nests in blue tits and 6 in great tits were removed). Each model was run with a binomial distribution and logit link function, and the models on EPO proportion were additionally weighed for the number of genotyped nestlings per brood. To avoid multicollinearity, reduce model overfitting, and improve the reliability of model selection, we grouped phenological and environmental predictors into six sets based on high correlations within sets (Pearson's $r>0.7$). Each set represented a distinct environmental or phenological gradient: timing of breeding, reproductive investment, tree characteristics, cavity entrance characteristics, cavity depth dimensions and nest space dimensions (Supporting Information:

Phenological and environmental variables; Figures S2 and S3). From each set, we retained the single variable recorded in the highest number of nests to maximize sample size. Specifically, we selected: lay date, initial brood size, height, entrance area, depth, and nest bottom area. These six retained variables were then used as continuous predictors in model selection procedures, to which we also included nest exposure (reflecting thermal properties of the nest and categorized as: coolest, intermediate, warmest) as a fixed predictor. The random effects structure included year and study plot (C, M, which had the most nests located, and 'Other'—encompassing study plots W and K with fewer nests present, to balance out the sample sizes). We then performed AIC-based model selection with the *dredge()* function of the MuMln package (Bartoń, 2022), as we did not have specific assumptions to test which of the variables would contribute most to variance in EPP.

To test the relationship between EPP and available proxies of reproductive success and parental investment, we fitted separate LMMs per species, with responses defined as the number of nestlings that successfully fledged and nestling fledging age. In models on the number of fledged nestlings in great tits, more nests failed than in blue tits—thus, great tit fledgling number was more often 0. To model great tit reproductive success, we therefore used a GLMM with Poisson distribution as it provided a better fit, and lower AIC compared with fitting a separate zero-inflation process. We again fitted the same set of seven predictors, that is, lay date, initial brood size, height, entrance area, depth, nest bottom area and exposure, alongside plot and year as random effects. Here, we fitted separate models, testing two EPP-related predictors to assess the qualitative (EPP occurrence—fixed) and quantitative (EPO proportion—continuous) contributors to reproductive success. We used the *fixed* condition in the *dredge()* function of the MuMln package (Bartoń, 2022) to ensure these were always retained as variables of interest.

To identify the most parsimonious model, we performed model selection using AICc-based model averaging (MuMln::dredge). Predictors appearing in models with $\Delta$AICc <1.5 (rather than 2 as usually applied; Burnham & Anderson, 2002 to reduce model uncertainty and to include fewer meaningful predictors from more robust models) were retained and refitted in a single linear mixed-effects model to obtain parameter estimates, random effect variance, multicollinearity diagnostics (VIF), and marginal/conditional $R^2$. Due to problems with a singular fit at some models' selection stage, we removed random effects for which the variance was zero or close to zero, before running *dredge()*. We calculated $R^2$ as defined in Nakagawa and Schielzeth (2013) for each model with the MuMln package (Bartoń, 2022). We tested all final models for multicollinearity (VIF scores never exceeded 2). We checked each model fit (KS test), dispersion and identified potential outliers by simulating residuals and plotting residual diagnostics using the DHARMa package (Hartig, 2024). Outliers were identified and tested by re-running the models on adjusted datasets. In all cases, removing outliers did not alter the factor estimates either qualitatively or quantitatively, so all outliers were retained. The only exception was the model examining the proportion of EPO in great tits, where outliers did affect model
<figure><img src="image_3.png" /></figure>

FIGURE 1 (a) Primeval forest EPP occurrence and proportion of EPO in the brood in blue tits (left) and great tits (right), alongside their environmental, life-history and phenological contributors. There was no relationship between EPP occurrence or proportion of EPO in the brood and reproductive success (number of fledged offspring) or parental investment (age of nestlings at fledging) in Białowieża primeval forest. Green upside arrows denote positive correlations, red downward arrows—negative correlations, and blue NS signs – no correlation between EPP/EPO and the environmental/life-history variable of interest (see Tables 1 and 2 and Tables S2 and S3). Extra-pair broods with recorded EPP occurrence (as %) (b); and EPO proportion (as %) (c) for blue tits (upper panels) and great tits (lower panels) reported across studies in (Brouwer & Griffith, 2019) (black dots). Values from primeval forest natural cavities are reported as large green dots. Values from the only other known natural cavity study (Di Lecce et al., 2023) are reported as small green dots. The magenta dashed line represents average values, while the box represents the median with IQR and whiskers representing 1.5 IQR. No y-axis variable is plotted; the vertical position of points (jitter) is used to visually separate EPP/EPO data points to enhance graphical clarity (top=blue tit, bottom=great tit).

estimates. However, these values fell within the expected biological
range and were not highly influential. Therefore, we chose to retain
them to preserve ecological validity.

3 | RESULTS AND DISCUSSION

This first report on EPP in two model avian species breeding in natu-
ral cavities within a primeval forest provides baseline EPP levels that

are instrumental for a better understanding of avian mating strat-
egies. We report that EPP in blue tits was recorded in 48% of all
broods (46 out of 96 broods were of mixed paternity), with 15% of
offspring being extra-pair (56 nestlings out of 362). In great tits, EPP
was recorded in 39% of broods (47 out of 122 broods were of mixed
paternity), with 14% of EPO (111 nestlings out of 821) (Figure 1a).
These levels are in line with values reported in more commonly stud-
ied environments—nestbox populations of blue tits and great tits
breeding in secondary forests (Figure 1b).
<table><caption>TABLE 1 Generalized linear models examining the occurrence of extra-pair offspring in blue tit and great tit broods in Białowieża primeval forest.</caption><caption>Occurrence of extra-pair offspring in the brood, GLM binomial</caption><thead><tr><th></th><th colspan="4" style="text-align:center;">(n=80 nests)</th><th colspan="4" style="text-align:center;">(n=99 nests)</th></tr><tr><th></th><th>Estimate</th><th>SE</th><th>z</th><th>p</th><th>Estimate</th><th>SE</th><th>z</th><th>p</th></tr></thead><tbody><tr><td>(Intercept)</td><td>-0.94</td><td>1.67</td><td>-0.56</td><td>0.576</td><td>0.41</td><td>0.80</td><td>0.51</td><td>0.611</td></tr><tr><td>Depth</td><td>0.09</td><td>0.04</td><td>2.09</td><td><strong>0.037</strong></td><td>-0.02</td><td>0.02</td><td>-0.99</td><td>0.323</td></tr><tr><td>Height</td><td>0.12</td><td>0.05</td><td>2.56</td><td><strong>0.010</strong></td><td>-0.05</td><td>0.07</td><td>-0.66</td><td>0.507</td></tr><tr><td>Entrance area</td><td>-0.03</td><td>0.03</td><td>-0.91</td><td>0.361</td><td>-0.01</td><td>0.01</td><td>-0.82</td><td>0.413</td></tr><tr><td>Initial brood size</td><td>-0.12</td><td>0.13</td><td>-0.95</td><td>0.344</td><td></td><td></td><td></td><td></td></tr><tr><td>Exposure (intermediate)</td><td></td><td></td><td></td><td></td><td>1.55</td><td>0.59</td><td>2.63</td><td><strong>0.009</strong></td></tr><tr><td>Exposure (coolest)</td><td></td><td></td><td></td><td></td><td>0.92</td><td>0.52</td><td>-1.78</td><td>0.075</td></tr><tr><td>Nest bottom area</td><td></td><td></td><td></td><td></td><td>-0.003</td><td>0.002</td><td>-1.28</td><td>0.199</td></tr><tr><th>Random effects</th><th>Variance</th><th>Variance</th><th>Variance</th><th>Variance</th></tr><tr><td>Year</td><td>0.000 (removed)</td><td>0.000 (removed)</td><td>0.000 (removed)</td><td>0.000 (removed)</td></tr><tr><td>Plot</td><td>0.000 (removed)</td><td>0.000 (removed)</td><td>0.000 (removed)</td><td>0.000 (removed)</td></tr><tr><td>R<sup>2</sup>m/R<sup>2</sup>c</td><td>0.253/0.253</td><td></td><td></td><td>0.143/0.143</td></tr></tbody></table>

Note: The response variable was modelled as absence (0) or presence (1). Predictor variables were selected based on model averaging across all models with ΔAICc <1.5. Exposure 'warmest' was the reference category for the fixed predictor. Significant differences (p < 0.05) are in bold. Marginal (R²m) and conditional (R²c) R-squared are shown.

To contextualize the values observed in Białowieża, we derived population-specific values of EPP in blue tits and great tits from a recent review by Brouwer and Griffith (2019), in which both species were by far the most studied avian model systems in terms of EPP (9.5k blue tit and 14k great tit individuals assessed). We also added EPP values from a recent nestbox-natural cavity comparison to reference the only other study on blue tits and great tits carried out in natural cavities so far (Di Lecce et al., 2023; Figure 1b). Since EPP values for both species reported from Białowieża fall within the middle of distributions for extra-pair broods and EPO ranges of the otherwise unexplained variation in EPP levels (Brouwer & Griffith, 2019), they likely represent a natural reference point that is advantageous across a range of environmental conditions, pointing to the evolutionary stability of mating systems. This finding validates and supports foundational evolutionary biology principles, suggesting that certain reproductive strategies are advantageous across diverse environments with a degree of behavioural plasticity. It affirms the universality of specific selection pressures (e.g. sexual selection, mate choice, or territorial behaviour), extending the generalizability of ecological models. Future studies can thus use this primeval forest baseline to estimate the magnitude of shifts in EPP relative to natural levels.

Remarkably, observed variation in EPP levels in a primeval forest did not affect reproductive success or parental investment in either

species (Figure 1a, Tables S2 and S3). This suggests weak or neutral selection on promiscuity under the ecological conditions of this primeval forest ecosystem. In Białowieża, this might stem from high pressures by various predators [in particular forest dormouse, great spotted woodpecker, yellow-necked mouse, pine marten and weasel (Wesołowski & Rowiński, 2012)], which are largely diminished in secondary forests or in nestbox populations (Wesołowski et al., 2015). Under such conditions, EPP levels might be at equilibrium in environments where predation is the major cause of reproductive failure—it accounts for 80% of nest failures in blue tits (Wesołowski & Rowiński, 2012) and 69% in great tits (Maziarz et al., 2016). Thus our data suggest that there is no fitness advantage to altering EPP rates in a primeval forest. Despite demonstrated genetic benefits of EPP (Westneat & Stewart, 2003), promiscuous copulations, competition for partners and mate guarding might attract predators to the nest and increase energy expenditure or time spent away from more defensive or survival-focused activities (Brouwer & Griffith, 2019). When predation dominates as a source of reproductive failure, the costs of activities associated with promiscuity may balance any potential genetic benefits, resulting in selection that neither favours nor disfavors current EPP levels. Since in Białowieża primeval forest, EPP was not affecting reproductive success or parental investment directly (Figure 1a, Tables S2 and S3), observed EPP levels likely reflect this balance between competing selective pressures rather than an optimum driven by genetic benefits alone. The importance
<table><caption>TABLE 2 Generalized linear (mixed) models examining the proportion of extra-pair offspring per brood in blue tit and great tit broods in Białowieża primeval forest.</caption><caption>Proportion of extra-pair offspring per brood, GLMM binomial</caption><thead><tr><th colspan="5">(n = 37 nests)</th><th colspan="5">(n = 83 nests)</th></tr><tr><th></th><th>Estimate</th><th>SE</th><th>z</th><th>p</th><th>Estimate</th><th>SE</th><th>z</th><th>p</th></tr></thead><tbody><tr><td>(Intercept)</td><td>-5.52</td><td>1.13</td><td>-4.88</td><td>&lt;0.0001</td><td>2.11</td><td>0.79</td><td>2.68</td><td>0.007</td></tr><tr><td>Depth</td><td>0.08</td><td>0.02</td><td>3.18</td><td>0.001</td><td>-0.03</td><td>0.01</td><td>-2.61</td><td>0.009</td></tr><tr><td>Height</td><td>0.09</td><td>0.03</td><td>3.13</td><td>0.002</td><td>-0.02</td><td>0.03</td><td>-0.66</td><td>0.509</td></tr><tr><td>Lay date</td><td>0.07</td><td>0.04</td><td>1.74</td><td>0.082</td><td></td><td></td><td></td><td></td></tr><tr><td>Entrance area</td><td></td><td></td><td></td><td></td><td>-0.02</td><td>0.01</td><td>-2.41</td><td>0.016</td></tr><tr><td>Initial brood size</td><td></td><td></td><td></td><td></td><td>-0.26</td><td>0.07</td><td>-3.76</td><td>0.0002</td></tr><tr><td>Exposure (intermediate)</td><td></td><td></td><td></td><td></td><td>0.46</td><td>0.27</td><td>1.67</td><td>0.095</td></tr><tr><td>Exposure (coolest)</td><td></td><td></td><td></td><td></td><td>-0.10</td><td>0.27</td><td>-0.36</td><td>0.714</td></tr><tr><td>Nest bottom area</td><td></td><td></td><td></td><td></td><td>-0.001</td><td>0.001</td><td>-1.11</td><td>0.267</td></tr><tr><th>Random effects</th><th>Variance</th><th>Variance</th><th>Variance</th><th>Variance</th></tr><tr><td>Year</td><td>0.114</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>Plot</td><td>0.000</td><td>0.000</td><td>0.000</td><td>0.000</td></tr><tr><td>R<sup>2</sup>m/R<sup>2</sup>c</td><td>0.165/0.193</td><td></td><td></td><td>0.136/0.136</td></tr></tbody></table>

Note: Predictor variables were selected based on model averaging across all models with ΔAICc <1.5. Exposure 'warmest' was the reference category for the fixed predictor. Significant differences (p < 0.05) are indicated in bold. Marginal (R²m) and conditional (R²c) R-squared are shown.

of minimizing predation was also confirmed by lower numbers of fledged blue tits from cavities that had larger nest entrance areas, which were easier to depredate (Table S3A).

In general, we show that EPP levels vary with nest characteristics associated with quality and safety, but these relationships are species-specific and sometimes opposite in direction between species (Figure 1a). In blue tits, EPP increased both in terms of occurrence and EPO proportion in higher-located nests (Tables 1 and 2). Nests higher in the tree may be harder for the social father to guard, but more importantly, increasing height has been previously demonstrated to be the largest negative contributor to nesting success in blue tits, as high nests were exposed to higher predation risk (Wesołowski & Rowiński, 2012). Surprisingly, in deeper blue tit nests, which also tend to be safer (Wesołowski & Rowiński, 2012), there was more EPP (Tables 1 and 2). In great tits, we observed less EPP in nests with multiple indicators of higher quality: in deeper nests (Table 2), which were previously shown to be more successful and harder to predate (Maziarz et al., 2016), and in nests with warmer exposures (Table 1), which provide a favourable microclimate that allows reduction of costs of incubation and nestling thermoregulation (Sudyka et al., 2023). Moreover, we observed a lower proportion of EPO in great tit nests with larger broods and larger entrances (Table 2). Entrance area was also shown to contribute positively to breeding success in Białowieża great tits (Maziarz et al., 2016). This reduction of EPP in good-quality nests can be readily explained by

prioritization of current reproductive success over promiscuity-related risks and has been reported previously (Raj Pant et al., 2022; Vedder et al., 2011; Whittingham & Dunn, 2005). The divergent results with respect to nest depth between species are difficult to explain, but nest depth is a much stronger factor in nesting success for great tits, whose nests are generally deeper than those of blue tits (Maziarz et al., 2016; Wesołowski & Rowiński, 2012). Different predator pressures may therefore influence how depth contributes to nest quality perception in each species.

Ultimately, EPP levels observed in this primeval forest provide a reference point for better understanding promiscuity in natural populations, aligning with levels observed in human-modified landscapes. The absence of detectable fitness consequences for EPP variation suggests that selective pressures on promiscuity are context-dependent and possibly influenced by ecological factors such as predation regimes. This work offers a baseline for future studies seeking to disentangle the ecological and evolutionary drivers of mating systems, particularly across gradients of human disturbance.

## AUTHOR CONTRIBUTIONS

Joanna Sudyka, Grzegorz Hebda, Patryk Rowiński, Tomasz Wesołowski, and Marta Szulkin conceived the ideas and designed the methodology; Grzegorz Hebda, Patryk Rowiński, Marta Cholewa, Marta Maziarz, and Tomasz Wesołowski collected the
data; Joanna Sudyka, Robert Rutkowski, and Anna Szczuka performed lab analyses; Irene Di Lecce and Charles Perrier performed parentage assignment, and Joanna Sudyka performed the statistical data analysis; Joanna Sudyka led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

## ACKNOWLEDGEMENTS

We thank members of the Białowieża fieldwork team and numerous field assistants for their contribution to data collection. We also thank Hannah Dugdale for her support in data analysis.

## FUNDING INFORMATION

The study was funded by the Polish National Science Centre grant no. OPUS 2016/21/B/NZ8/03082 awarded to Marta Szulkin. Fieldwork in 2005–2007 was funded by a grant from the Polish Ministry of Science and Higher Education no. 2PO4C 075 27 awarded to Patryk Rowiński.

## CONFLICT OF INTEREST STATEMENT

The authors declare no competing interests.

## DATA AVAILABILITY STATEMENT

The R code used in this study and raw data can be found at https://zenodo.org/records/15487049 (Sudyka et al., 2025).

## STATEMENT ON INCLUSION

Our study brings together authors from a number of different countries, including scientists based in the country where the study was carried out. All authors were engaged early on with the research and study design to ensure that the diverse sets of perspectives they represent were considered from the onset. Whenever relevant, literature published by scientists from the region was cited; relevant work published in the local language was considered.

## ORCID

Joanna Sudyka https://orcid.org/0000-0003-3006-3806

## REFERENCES

* Bartón, K. (2022). MuMln: Multi-Model Inference. R package version 1.47.1 (1.47.1). https://cran.r-project.org/package=MuMln

* Brouwer, L., & Griffith, S. C. (2019). Extra-pair paternity in birds. *Molecular Ecology*, 28(22), 4864–4882. https://doi.org/10.1111/mec.15259

* Burnham, K. P., & Anderson, D. R. (2002). *Model selection and multi-model inference a practical information-theoretic approach*. Springer. https://doi.org/10.1017/CBO9781107415324.004

* Di Lecce, I., Perrier, C., Szulkin, M., & Sudyka, J. (2023). Extra-pair paternity, breeding density, and synchrony in natural cavities versus nestboxes in two passerine birds. *Ecology and Evolution*, 13(6), 1–16. https://doi.org/10.1002/ece3.10163

* Di Lecce, I., Sudyka, J., Perrier, C., & Szulkin, M. (2024). Extra-pair paternity in two passerine birds breeding in a gradient of urbanisation. *Molecular Ecology*, 33(17), e17481. https://doi.org/10.1111/mec.17481

* FAO. (2020). *Global forest resources assessment 2020*. FAO. https://doi.org/10.4060/ca9825en

* Griffith, S. C., Owens, I. P. F., & Thuman, K. A. (2002). Extra-pair paternity in birds: A review of interspecific variation and adaptive function. *Molecular Ecology*, 11(11), 2195–2212. https://doi.org/10.1046/j.1365-294X.2002.01613.x

* Gruber, B., Unmack, P. J., Berry, O. F., & Georges, A. (2018). dartR: An R package to facilitate analysis of SNP data generated from reduced representation genome sequencing. *Molecular Ecology Resources*, 18(3), 691–699. https://doi.org/10.1111/1755-0998.12745

* Hartig, F. (2024). DHARMa: Residual diagnostics for hierarchical (multi-level/mixed) regression Models. R package version 0.4.7.

* Keedwell, R. J. (2003). Does fledging equal success? Post-fledging mortality in the black-fronted tern. *Journal of Field Ornithology*, 74(3), 217–221. https://doi.org/10.1648/0273-8570-74.3.217

* Kempenaers, B., Borgström, P., Loës, P., Schlicht, E., & Valcu, M. (2010). Artificial night lighting affects dawn song, extra-pair siring success, and lay date in songbirds. *Current Biology*, 20(19), 1735–1739. https://doi.org/10.1016/j.cub.2010.08.028

* Kilian, A., Wenzl, P., Huttner, E., Carling, J., Xia, L., Blois, H., Caig, V., Heller-Uszynska, K., Jaccoud, D., Hopper, C., Aschenbrenner-Kilian, M., Evers, M., Peng, K., Cayla, C., Hok, P., & Uszynski, G. (2012). Diversity arrays technology: A generic genome profiling technology on open platforms. *Life Sciences*, 531(1), 67–89. https://doi.org/10.1007/978-1-61779-870-2_5

* Marris, E. (2008). Ecology: The heart of the wood. *Nature*, 455(7211), 277–280. https://doi.org/10.1038/455277a

* Martin, T. E., Tobalske, B., Riordan, M. M., Case, S. B., & Dial, K. P. (2018). Age and performance at fledging are a cause and consequence of juvenile mortality between life stages. *Science Advances*, 4(6), 1–9. https://doi.org/10.1126/sciadv.aar1988

* Mayer, C., & Pasinelli, G. (2013). New support for an old hypothesis: Density affects extra-pair paternity. *Ecology and Evolution*, 3(3), 694–705. https://doi.org/10.1002/ece3.489

* Maziarz, M., Wesołowski, T., Hebda, G., Cholewa, M., & Broughton, R. K. (2016). Breeding success of the Great Tit Parus major in relation to attributes of natural nest cavities in a primeval forest. *Journal of Ornithology*, 157(1), 343–354. https://doi.org/10.1007/s10336-015-1294-2

* Mo, L., Zohner, C. M., Reich, P. B., Liang, J., de Miguel, S., Nabuurs, G.-J., Renner, S. S., van den Hoogen, J., Araza, A., Herold, M., Mirzagholi, L., Ma, H., Averill, C., Phillips, O. L., Gamarra, J. G. P., Hordijk, I., Routh, D., Abegg, M., Adou Yao, Y. C., & Crowther, T. W. (2023). Integrated global assessment of the natural forest carbon potential. *Nature*, 624, 92–101. https://doi.org/10.1038/s41586-023-06723-z

* Nakagawa, S., & Schielzeth, H. (2013). A general and simple method for obtaining R2 from generalized linear mixed-effects models. *Methods in Ecology and Evolution*, 4(2), 133–142. https://doi.org/10.1111/j.2041-210x.2012.00261.x

* Nilforooshan, M. A., & Saavedra-Jiménez, L. A. (2020). Ggroups: An R package for pedigree and genetic groups data. *Hereditas*, 157(1), 17. https://doi.org/10.1186/s41065-020-00124-2

* Norris, A. R., Aitken, K. E. H., Martin, K., & Pokorny, S. (2018). Nest boxes increase reproductive output for Tree Swallows in a forest grassland matrix in central British Columbia. *PLoS One*, 13(10), 1–12. https://doi.org/10.1371/journal.pone.0204226

* Perrier, C., Delahaie, B., & Charmantier, A. (2018). Heritability estimates from genomewide relatedness matrices in wild populations: Application to a passerine, using a small sample size. *Molecular Ecology Resources*, 18(4), 838–853. https://doi.org/10.1111/1755-0998.12886

* Purcell, K. L., Verner, J., & Oring, L. W. (1997). A comparison of the breeding ecology of birds nesting in boxes and tree cavities. *Auk*, 114(4), 646–656. https://doi.org/10.2307/4089284

* R Core Team. (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing. (4.2.2).
R Core Team. (2023). *R: A language and environment for statistical computing*. R Foundation for Statistical Computing. (4.3.2). https://www.r-project.org/

Raj Pant, S., Versteegh, M. A., Hammers, M., Burke, T., Dugdale, H. L., Richardson, D. S., & Komdeur, J. (2022). The contribution of extra-pair paternity to the variation in lifetime and age-specific male reproductive success in a socially monogamous species. *Evolution*, 76(5), 915-930. https://doi.org/10.1111/evo.14473

Sudyka, J., Di Lecce, I., Hebda, G., Rowiński, P., Rutkowski, R., Szczuka, A., Perrier, C., Cholewa, M., Maziarz, M., Wesołowski, T., & Szulkin, M. (2025). Avian promiscuity in the last European primeval forest [Data set]. Zenodo. https://doi.org/10.5281/zenodo.15487049

Sudyka, J., Di Lecce, I., & Szulkin, M. (2023). Microclimate shifts in nest-boxes and natural cavities throughout reproduction. *Journal of Avian Biology*, 2023(1-2), 1-17. https://doi.org/10.1111/jav.03000

Sudyka, J., Di Lecce, I., Wojas, L., Rowiński, P., & Szulkin, M. (2022). Nest-boxes alter the reproductive ecology of urban cavity-nesters in a species-dependent way. *Journal of Avian Biology*, 2022(11-12), 1-15. https://doi.org/10.1111/jav.03051

UNESCO World Heritage Centre. (n.d.). *Białowieża Forest*. https://whc.unesco.org/en/list/33/

Vedder, O., Komdeur, J., & Van Der Velde, M. (2011). Polygyny and extra-pair paternity enhance the opportunity for sexual selection in blue tits. *Behavioral Ecology and Sociobiology*, 65, 741-752. https://doi.org/10.1007/s00265-010-1078-x

Watson, J. E. M., Evans, T., Venter, O., Williams, B., Tulloch, A., Stewart, C., Thompson, I., Ray, J. C., Murray, K., Salazar, A., McAlpine, C., Potapov, P., Walston, J., Robinson, J. G., Painter, M., Wilkie, D., Filardi, C., Laurance, W. F., Houghton, R. A., ... Lindenmayer, D. (2018). The exceptional value of intact forest ecosystems. *Nature Ecology & Evolution*, 2(4), 599-610. https://doi.org/10.1038/s41559-018-0490-x

Weatherhead, P. J., & Dufour, K. W. (2000). Fledging success as an index of recruitment in red-winged blackbirds. *Auk*, 117(3), 627-633. https://doi.org/10.2307/4089588

Wesołowski, T. (2011). Reports from nestbox studies: A review of inadequacies. *Acta Ornithologica*, 46(1), 13-17. https://doi.org/10.3161/000164511X589866

Wesołowski, T., Czeszczewik, D., Hebda, G., Maziarz, M., Mitrus, C., & Rowiński, P. (2015). 40 years of breeding bird community dynamics in a primeval temperate forest (Białowieża National Park, Poland). *Acta Ornithologica*, 50(1), 95-120. https://doi.org/10.3161/00016454A02015.50.1.010

Wesołowski, T., Mitrus, C., Czeszczewik, D., & Rowiński, P. (2010). Breeding bird dynamics in a primeval temperate forest over thirty-five years: Variation and stability in the changing world. *Acta Ornithologica*, 45(2), 209-232. https://doi.org/10.3161/000164510X551354

Wesołowski, T., & Rowiński, P. (2012). The breeding performance of Blue Tits *Cyanistes caeruleus* in relation to the attributes of natural holes in a primeval forest. *Bird Study*, 59(4), 437-448. https://doi.org/10.1080/00063657.2012.722189

Wesołowski, T., & Rowiński, P. (2014). Do blue tits Cyanistes caeruleus synchronize reproduction with caterpillar peaks in a primeval forest? *Bird Study*, 61, 37-41. https://doi.org/10.1080/00063657.2014.899307

Wesołowski, T., Tomiałojć, L., Mitrus, C., Rowiński, P., & Czeszczewik, D. (2002). The breeding bird community of a primeval temperate forest (Białowieża National Park, Poland) at the end of the 20th century. Acta Ornithologica, 37(1), 27-44. https://doi.org/10.3161/068.037.0105

Westneat, D. F., & Stewart, I. R. K. (2003). Extra-pair paternity in birds: Causes, correlates, and conflict. *Annual Review of Ecology, Evolution, and Systematics*, 34, 365-396. https://doi.org/10.1146/annurev.ecolsys.34.011802.132439

Whittingham, L. A., & Dunn, P. O. (2005). Effects of extra-pair and within-pair reproductive success on the opportunity for selection in birds. *Behavioral Ecology*, 16(1), 138-144. https://doi.org/10.1093/beheco/arh140

Yang, J., Lee, S. H., Goddard, M. E., & Visscher, P. M. (2011). GCTA: A tool for genome-wide complex trait analysis. *The American Journal of Human Genetics*, 88(1), 76-82. https://doi.org/10.1016/j.ajhg.2010.11.011

Zheng, X., Levine, D., Shen, J., Gogarten, S. M., Laurie, C., & Weir, B. S. (2012). A high-performance computing toolset for relatedness and principal component analysis of SNP data. *Bioinformatics*, 28(24), 3326-3328. https://doi.org/10.1093/bioinformatics/bts606

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**Figure S1.** Map of study sites within Białowieża National Park.

**Figure S2.** Correlation heatmap of the full set of phenological and environmental variables gathered in Białowieża primeval forest blue tit nests.

**Figure S3.** Correlation heatmap of the full set of phenological and environmental variables gathered in Białowieża primeval forest great tit nests.

**Figure S4.** Call rate by individual (a) and by locus (b) in adult blue tits and great tits.

**Figure S5.** Comparison of genetic relatedness calculated with 40% and 100% of the SNP data in blue tits and great tits.

**Figure S6.** GRM relatedness values.

**Table S1.** Sample sizes.

**Table S2.** Generalized and general linear (mixed) models examining variation in (a) number of fledged young and (b) fledging age with EPP occurrence.

**Table S3.** Generalized and general linear (mixed) models examining variation in (a) number of fledged young and (b) fledging age with EPO proportion per nest.

**How to cite this article:** Sudyka, J., Di Lecce, I., Hebda, G., Rowiński, P., Rutkowski, R., Szczuka, A., Perrier, C., Cholewa, M., Maziarz, M., Wesołowski, T., & Szulkin, M. (2025). Avian extra-pair paternity in the last European primeval forest. *Functional Ecology*, 00, 1-9. https://doi.org/10.1111/1365-2435.70223