Methods

Ethics. Our study was purely observational and conducted in accordance with the national regulations of the countries where it took place. All the experimental procedures were approved by the Lower Saxony State Office for Consumer Protection and Food Safety (LAVES: Niedersächsisches Landesamt für Verbraucherschutz und Lebensmittelsicherheit) in Germany. Moreover, the procedures were approved by the research coordinators at Osnabrück Zoo and WildtierPark Edersee (hereafter WPE) in Kellerwald, both in Germany, and by the directors of Riserva Naturale Regionale Lago di Penne (RNRLP), Parco Faunistico del Monte Amiata (PFMA) in Arcidosso and Centro di Recupero Animanatura Wild Sanctuary (CRAWS) in Semproniano, all in Italy. Risk assessments were carried out collaboratively with the keepers of each study group, with approval granted only if the procedures were deemed to pose no risk to the animals. As all facilities allowed visitors and our observations followed standard visitor practices (i.e., observing the animals without any experimental intervention), the study was considered to pose no adverse effects to the animals.

Study subjects. We studied two groups of non-admixed wolves (N = 20) and three groups of admixed wolves (N = 13). All individuals in the admixed groups exhibited various combinations of anomalous morphological traits (e.g., dark or pale coats, dewclaws, depigmented nails; Ciucci et al. 2003, Galaverni et al. 2017). Genetic analysis using biparental and/or uniparental markers confirmed they were introgressed individuals from second- or later-generation backcrosses to wolves (Caniglia et al. 2020). At the time the observations, all individuals had reached sexual maturity and were considered adults (see Table 1, for the sex and age composition of the five study groups). All admixed wolves had been removed from the wild for conservation management and relocated to captivity, following implementation of management responses aimed at preventing hybridization in wild wolf populations (Bocci et al. 2015). All the study subjects lived in captive conditions, in relatively large enclosures with their conspecifics. In all groups, animals were exposed to visitors and keepers, who usually fed them every one to three days by placing dead chickens, rabbits and other potential prey in their enclosures.

Group | Location | Subject | Sex | Age | Rank | Centrality
Non-admixed wolves | Osnabruck | Angel | Female | 5 | 0.000 | 0.627
Ano | Male | 5 | 0.033 | 0.713
Eren | Male | 7 | 0.224 | 1.000
Grey | Female | 4 | 0.015 | 0.982
Jackal | Female | 4 | 0.311 | 0.805
Kaz | Male | 7 | 0.488 | 0.850
Lagerta | Male | 5 | 0.549 | 0.737
Noname | Female | 5 | 0.83 | 0.794
One | Female | 5 | 0.205 | 0.682
Ragnar | Male | 4 | 1.000 | 0.920
Snow | Male | 7 | 0.123 | 0.882
Witcher | Female | 4 | 0.138 | 0.410
Kellerwald | Amik | Male | 5 | 0.000 | 0.219
Anouk | Male | 5 | 0.207 | 0.727
Bagira | Male | 4 | 0.188 | 0.020
Cara | Female | 2 | 0.213 | 0.431
Charly | Male | 2 | 0.787 | 0.656
Cowalski | Male | 2 | 0.515 | 0.637
Noki | Female | 8 | 0.816 | 1.000
Nuk | Male | 9 | 1.000 | 0.845
Dog-admixed wolves | Penne | Baloo | Male | 7 | 0.048 | 0.625
Balto | Male | 7 | 0.657 | 1.000
Licaone | Male | 7 | 0.000 | 0.390
Lucio | Male | 7 | 0.454 | 0.812
Neo | Male | 7 | 1.000 | 0.825
Syd | Male | 7 | 0.554 | 0.707
Semproniano | Bo | Male | 11 | 0.481 | 1.000
Nero | Male | 11 | 1.000 | 0.965
Vincent | Male | 11 | 0.000 | 0.978
Arcidosso | Koda | Female | ≥ 8 | 0.204 | 0.684
Luna | Female | 10 | 0.988 | 0.949
Ora | Female | ≥ 8 | 1.000 | 1.000
Rossa | Female | 10 | 0.000 | 0.000

The Osnabrück group included 12 Hudson wolves (Canis lupus hudsonicus): three males from the same litter who were born in the Amsterdam Zoo (Netherlands), and nine individuals (i.e., six females and tree males) from two different litters that were born in the Osnabrück Zoo. Three of the six males had been vasectomized at the time of the study (i.e., testes were intact and gonadal endocrine function was expected to be maintained). Their 0.45 ha enclosure included open woods and a walkway along the edges, which allowed a clear view on the group. The WPE group included 8 European wolves (Canis lupus lupus): a breeding couple and part of their offspring from three previous litters. The 1.04 ha enclosure included grassland and shrubs, and a steep forested hillside, where wolves were not always visible. All the individuals in the group were still intact at the time of the study.

The RNRLP group included 6 male admixed wolves, all siblings, captured as pups at about 1 month of age in 2017 in Gran Sasso–Monti della Laga National Park from a dog-admixed wolf pack. They were introgressed at nuclear loci and/or at the K-locus and all shared a dog Y-haplotype (Caniglia et al. 2020), most likely representing second- or later generation backcrosses to wolves. They were initially raised in kennel boxes but relocated, at 2.5 months of age, to the RNRLP enclosure, which measured approximately 0.50 ha and included a pond and a grove. All were vasectomized at six months of age and therefore were hormonally intact. The group could be observed from an elevated visitor area that provided visual access to most of the enclosure. The CRAWS group included 3 male admixed wolves, all siblings from a litter of 10 that who was captured in May 2014 near Siena, when the pups had less than one week of age (Zingaro et al. 2014). Similar to the RNRLP group, they were introgressed at nuclear loci and/or at the K-locus, and all shared a dog Y-haplotype (Caniglia et al. 2020), most likely representing third- or later generation backcrosses to wolves. After capture, they were initially hand raised and bottle-fed by caregivers for about one month, and six (three males and three females) were then transferred to the CRAWS, where the three males had always remained together. All three males were neutered through orchiectomy during their first year of life. Their enclosure measured around 0.30 ha and included a grass clearing, a hill and, beyond the hill, another green area with trees that was not visible to the observer, who conducted the observations in front of the clearing. The PFMA group included 4 admixed wolves, all females, who were captured in different occasions in Tuscany (G. Romeo, Tuscany Region, pers. comm.). Two sisters, along with a third that died prior to the onset of this study, were from the wild admixed litter captured in 2014 near Siena and hosted in the CRAWS, from where they were transferred at 8 years of age. Two additional females were live captured as adults in 2012, one near Firenze and one near Grosseto, and were housed in separate enclosures until joining the two previous sisters in 2016 and 2022. In addition to anomalous morphological traits, all females were introgressed at the nuclear loci, likely representing third- or later generation backcrosses to wolves (Caniglia et al. 2020). All females had been sterilized, the two sisters via ovariohysterectomy, and the other two individuals via tubal ligation (Amici et al. 2024). The PFMA enclosure measured approximately 9.00 ha and included a grass clearing atop a hill surrounded by woods. Observations were conducted from a visitor turret, which provided sight of only a portion of the enclosure.

Data collection and coding. We conducted behavioural observations in Osnabrück from July to December 2023, at the WPE in Kellerwald from January to May 2024, at the RNRLP in Penne from April to August 2024, at the PFPMA in Arcidosso from September to October 2024, and at the CRAWS in Semproniano from May to September 2025. In each group, we collected three main types of data. First, we used all-occurrence sampling to assess group hierarchy and determine individual ranks. We recorded all visible instances of aggression (i.e., bites, lunges, chases), dominance behaviours (i.e., standing over a recumbent individual in a T-position or placing a paw on it while standing for more than two seconds), and submissive behaviours (e.g., licking a partner while maintaining a hunched posture and/or a tucked tail; Frézard and Pape 2003; Packard 2003; Pifarré et al. 2012). For each interaction, we noted the identities of the winner and loser, as well as the date. Second, we conducted 15-min individual focal observations to assess dyadic bond strength and individual social prominence, based on the proportion of time spent in proximity to other group members. Each study subject was observed in 22 focal sessions distributed across different days in a pseudorandomized order and at varying times of day to minimize temporal biases.

During focal sessions, we recorded the exact duration that the focal individual spent within 5 m of other group members, noting their identities. The 5 m threshold was selected based on previous observations of the study groups suggesting that this distance reliably captured affiliative social proximity between individuals while excluding more incidental spatial associations (Amici et al. 2024). Third, whenever we detected a yawning event (Fig. 1), we recorded the identities of all visible group members and coded (i) whether they observed the yawn (i.e., had their head oriented toward the yawner during the event) or not, and (ii) whether they yawned or not within the subsequent two minutes. If it was not possible to understand with reasonable certainty whether another group member had observed the yawn or not, we did not include the observation in the dataset (Fig. 1). We selected a two-minute response window based on previous findings indicating that most contagious yawning in wolves occurs shortly after the triggering event (Romero et al. 2014), and because longer observation windows would have reduced reliability when monitoring multiple individuals simultaneously. Most yawning events were recorded during focal observations, although additional events were collected ad libitum to increase sample size and robustness.

Fig. 1: Yawn event in the non-admixed wolf group at the Osnabrück Zoo, Germany. The second wolf from the left was the yawner; the first and fourth wolves from the left were considered as having observed the yawn (i.e., their head was oriented toward the yawner during the yawn event); the third wolf from the left, in the background, was not included in the dataset, as it was not possible to understand with reasonable certainty whether the wolf had observed the yawn event or not; and the last wolf from the left was considered as not having observed the yawn

Statistical analyses. We conducted all analyses in R (version 4.4.2; R Core Team 2024). First, we determined the dominance hierarchy and extracted individual ranks (Table 1) using the Elo-rating method (EloRating package, version 0.43), based on all dyadic agonistic interactions with a clear winner-loser outcome Neumann and Kulik (2025). Elo scores were standardized to range from 0 to 1, with 0 corresponding to low-ranking individuals and 1 to high-ranking individuals. Second, we assessed dyadic bond strength using the Composite Sociality Index (CSI; Silk et al. 2009). CSI values were calculated by dividing the proportion of time two individuals spent in proximity by the mean proportion of time all dyads within the same group spent in proximity. CSI scores below 1 indicated lower-than-average relationship strength and values above 1 indicated higher-than-average strength. Third, we quantified individuals’ social network centrality by constructing adjacency matrices based on the time dyads spent in proximity. We used the packages asnipe (version 1.1.17; Farine 2013) and igraph (version 2.0.3; Csardi et al. 2006) to calculate individual centrality measures from these matrices (Table 1). Centrality scores ranged from 0 to 1 and reflected an individual’s level of social integration within the group: lower values indicated weaker integration, whereas higher values indicated stronger or more numerous social connections (Farine 2017).

Finally, we used generalized linear mixed models (GLMMs) to test our predictions. We constructed the dataset by entering one data point for each yawning event and each visible group member that could potentially observe it (N = 775). For each data point, we recorded: (i) the observation day (as multiple yawning events could occur on the same day); (ii) the event identity (as multiple partners could be associated with a single event); (iii) the identity of the yawner; (iv) the identity of the partner; (v) whether the partner had observed the event (0/1); and (vi) whether the partner yawned within the subsequent two minutes (0/1). We additionally included (vii) the CSI for each yawner-partner dyad, as well as the yawner’s (viii) rank and (ix) centrality as measures of social relevance within the group. Including both observers and non-observers allowed us to account for the baseline probability of yawning occurring independently of any contagion effect.

We fitted GLMMs using the glmmTMB package (Magnusson et al. 2021) with a binomial error distribution to model the log-odds of a partner yawning within two minutes of the initial event. As test predictors, we included three separate three-way interactions between introgression (admixed vs. non-admixed wolf, observation of the event (0/1), and (i) CSI, (ii) yawner rank and (iii) yawner centrality. As control predictors, we included yawner and partner sex (given evidence that female wolves may be more likely to exhibit behavioural contagion; Romero et al. 2014). Random intercepts included yawner identity, partner identity and event identity nested within observation day. Group identity was not included as a random effect because the number of groups was too small to reliably estimate variance components. We compared each full model to a corresponding null model that excluded the test predictors, using likelihood ratio tests (Dobson and Barnett 2008), available as R function ‘‘anova’‘. When the full–null comparison was significant, we used the drop1 function to assess the contribution of individual test predictors. If the three-way interaction was not significant, we refitted the model excluding that interaction and retaining the corresponding two-way interactions. We conducted post hoc pairwise comparisons for significant categorical interactions using the emmeans package with Tukey adjustments (Lenth 2020), with estimates back-transformed to the response scale for easier interpretability. We performed model diagnostics using the DHARMa (Hartig 2022) and performance (Lüdecke et al. 2021) packages. We detected no issues regarding residual distribution, model convergence, overdispersion or multicollinearity (maximum VIF across models = 1.68; Miles 2005).