Decrease in social cohesion in a colonial seabird under a perturbation regime

Social interactions, through influence on behavioural processes, can play an important role in populations’ resilience (i.e. ability to cope with perturbations). However little is known about the effects of perturbations on the strength of social cohesion in wild populations. Long-term associations between individuals may reflect the existence of social cohesion for seizing the evolutionary advantages of social living. We explore the existence of social cohesion and its dynamics under perturbations by analysing long-term social associations, in a colonial seabird, the Audouin’s gull Larus audouinii, living in a site experiencing a shift to a perturbed regime. Our goals were namely (1) to uncover the occurrence of long-term social ties (i.e. associations) between individuals and (2) to examine whether the perturbation regime affected this form of social cohesion. We analysed a dataset of more than 3500 individuals from 25 years of monitoring by means of contingency tables and within the Social Network Analysis framework. We showed that associations between individuals are not only due to philopatry or random gregariousness but that there are social ties between individuals over the years. Furthermore, social cohesion decreased under the perturbation regime. We sustain that perturbations may lead not only to changes in individuals’ behaviour and fitness but also to a change in populations’ social cohesion. The consequences of decreasing social cohesion are still not well understood, but they can be critical for the population dynamics of social species.

Scientific Reports | (2020) 10:18720 | https://doi.org/10.1038/s41598-020-75259-3 www.nature.com/scientificreports/ However, many colonial species are philopatric, thus annual association in breeding places between individuals may not necessarily reflect individual social ties but a shared tendency to breed in the same birthplace (i.e. natal philopatry) or breeding place (i.e. breeding philopatry). Additionally, the decision to join or leave a particular group may be based on general preferences for being with others (gregariousness) or preferences for particular individuals within the subgroup 28 (i.e. social ties). Thus the challenge lies in disentangling whether an annual association between individuals is only due to philopatry (natal or breeding), to random gregariousness, or also due to the existence of social ties within groups of individuals over time 31 . The existence of social ties between neighbouring pairs in breeding colonies are rarely considered in behavioural and ecological studies 32 and, if true, such associations may suggest the evolution of social cohesion for exploiting the evolutionary advantages of social living (including social information sharing) for individual fitness prospects. Many colonial seabirds evolved in very variable environments, with very patchy and unpredictable resources (both food and safe breeding habitats), which means that they are forced to make decisions constantly, and social coping and social information sharing may be more relevant than previously thought 4,24 . Social network theory, which originated in sociology to study human relationships and social organization 11,33,34 now provides both a conceptual framework and the analytical tools to explore social cohesion and social processes in animal populations 20,[35][36][37][38] . Network theory is now being simultaneously developed in a number of fields, including statistical physics, sociology, molecular biology, and computer science. As a result, the field is changing at a rapid pace. While not all developments can or should be applied toward the study of animal societies 39 , this rush of novel ideas from other disciplines is enriching behavioural ecology 40 .
To assess the existence of social cohesion and its dynamics under perturbations, we studied an ecological system involving a colonial, social vertebrate (the Audouin's gull Larus audouinii) living in a site experiencing a shift to a perturbed regime. Interestingly from a social point of view, the species breeds aggregated in spatially discrete patches within large colonies. Each breeding season, some patches go extinct and some are colonized 41,42 , with a patch occupation expectancy over time of 1.5 years 24 , forcing individuals to breed in patches different from the ones they were born in or they bred in the previous year. These colonization-extinction processes dismiss the possibility of aggregation only due to natal or breeding philopatry. Additionally, all patches have very similar habitat features, dismissing also the possibility of aggregation due to similar habitat choices between individuals, an important issue when analysing social aggregation in many study systems. Thus, our particular study system may allow us to disentangle whether annual breeding aggregation among individuals is an annual random association, or it rather results from social cohesion among individuals over the years (i.e. there are social ties between individuals).
An extensive long-term individual monitoring program has been carried out since 1988 at the Ebro Delta, including the main breeding site for the species during the last two decades, the Punta de la Banya (40° 34′ 48.0" N 0° 41′ 24.0" E) 43 . At this breeding site, population dynamics have undergone different phases: an initial growing phase after site colonization, a stable phase of dynamic equilibrium, and a final transition phase to population collapse 24,44 (Fig. 1). This collapse was due to the arrival of terrestrial predators, which led this colony to hold from 70% to only 3% of the total world population in only a decade (32% mean annual decline) 41,43 (Fig. 1).
The perturbation regime caused changes in the spatial distribution of patches at the site, changes in age structure, decrease in fecundity, and a progressive increase of dispersal to other sites 7,41 . The response of this population to predators has not been immediate probably due to strong philopatry, high site-suitability inertia, and social behavioural processes, such as conspecific attraction 24 . This raise in dispersal was probably caused by social processes, as social copying 24 , however it remains to assess how social cohesion among individuals, if it occurred, was affected by dispersal processes. One possibility is that dispersal would break social cohesion by individual heterogeneity in the willingness to disperse 45 . In contrast, social cohesion can be maintained over time when dispersal between breeding patches occurs collectively at the scale of social groups to the same sites 31,46 . Taking advantage of the long-term monitoring of this species, the knowledge of its population dynamics, and the use of tools recently developed in the Social Network Analysis (SNA) framework, we specifically addressed the following questions: (1) are there any long-term social ties between individuals breeding in the same patch?; and (2) have perturbations, in this case a perturbation regime, affected social cohesion? We finally discuss the role and consequences of social cohesion in population dynamics and resilience in social species.

Results
We analysed a total of 1,610,922 dyadic interactions during the first period (2002-2011) and 368,142 during the second period (2012-2017) (Supplementary Table S1).
Occurrence of social ties that persist over time. When assessing the social ties with the contingency table approach during the period of stability, the assumption that breeding aggregations in Audouin's gull were at random was rejected, and those individuals that bred in the same patch during the sub-period 2002-2006 had a higher probability of breeding together during the sub-period 2007-2011 ( χ 2 1 = 64.685, P < 0.0001) (Supplementary Table S1). When randomly reducing the sample size of the data set, results were still statistically significant in more than 95% of the cases (1000 randomizations).
Accordingly, when assessing the social ties with dependent regression terms in the AME function, we showed that the probability of breeding together during the second sub-period (2007-2011) depended on whether they have bred previously together in the first sub-period (2002)(2003)(2004)(2005)(2006), with a statistically significant coefficient of regression parameter (Table 1; Fig. 2). When analysing this data set with permutated data, we concluded that there was indeed a non-random association of individuals within the patches, with our statistic being among the 5% extreme values.

Scientific Reports
| (2020) 10:18720 | https://doi.org/10.1038/s41598-020-75259-3 www.nature.com/scientificreports/ Impact of perturbations on social cohesion. When we analysed the social ties during the transition to collapse phase, we observed that the probability of breeding together during the period 2012-2017 did not depend on whether they have bred in the same patch the five previous years (χ 1 2 = 1.8143, P value = 0.178) (Supplementary Table S1) and we could not reject the hypothesis of a random association between individuals during that period. Besides, the SNA approach showed that breeding aggregations in Audouin's gull during the transition phase period did not depend on whether they have bred in the same patch the five previous years (Table 1; Fig. 2).

Discussion
By studying a particular ecological system of a colonial long-lived species that experiences a perturbation regime, we showed that individuals do not associate at random but there are social ties among individuals that persist over the years when conditions are unperturbed. Additionally, we observed that perturbations may decrease social cohesion in animal populations.
The characteristic breeding behaviour of the study species that aggregate in patches that change annually, allowed us to show that there is some group stability, with individuals establishing social ties that persist over time. Even if we cannot rule out that habitat preference is behind the observed individual non-random aggregations, we consider this hypothesis improbable in our study system as there is no obvious habitat differences between the patches. Our study system resembles what was recorded for Slender-billed gulls (Chroicocephalus genei), a colonial breeder with weak inter-annual breeding-site fidelity: some individuals bred in the same patch across breeding seasons and some social groups showed tenacity despite the colony often moving every year 46 .  44 . Red arrow indicates the arrival of predators (mainly foxes) to the colony. The shadow area shows the temporal window for data analysis for comparing the periods before and over the collapse (separated by a red dashed line). Table 1. Results of the AME regression function to test if there were social ties between individuals while breeding during the stable period and during the transition phase to the collapse period. The null hypothesis is that individuals aggregate annually at random for breeding. ".dyad": coefficient of the dependent regression term considering the previous dyadic relationship between individuals; "pmean": posterior mean estimate; "psd": posterior standard deviation; "z-stat ": nominal z-score. www.nature.com/scientificreports/ Group stability can emerge as a product of network self-organization, although may then provide the necessary conditions for the evolution of other social processes 47,48 . Our results would support the idea that social aggregation during breeding would provide other advantages than the mere defence against predators 49,50 , such as social information sharing (e.g. [51][52][53], especially in colonial seabirds that evolved in changing and unpredictable environments 42 . Social information sharing is also crucial for decision-making in risky behaviours, such as dispersal, and previous studies showed that the perturbed regime in our study site caused dispersal to other sites, including the colonization of new habitats 41,54 . In our case study, sociality may have played a major role in driving dispersal and thus population dynamics, both during the exponential growth after colonization and the collapse after the perturbation regime 24,42 . This idea is also reinforced with a mechanistic dynamical model that shows that population dynamics of Audouin's gulls at the study site can only be explained by dispersal runaway caused by social copying 24 . The importance of social information compared to private information increases under perturbations, even when the quality of social information does not increase compared to a non-perturbed regime 55,56 . For instance, Maldonado et al. 57 show that experimental disturbances applied to a social bird may impact its foraging efficiency by social instability caused by the split of social groups. In colonial birds, breeding failure, which is a proxy of environmental stress, may trigger splitting of the social groups (e.g. 46 ). At the demographic level, the alteration of social network structure may affect the behaviour of populations. For instance, under stress conditions, sociality may operate through feedback loops such as social copying for dispersal, causing non-linear population dynamics and playing a critical role in the resilience of populations (e.g. 24 ). We showed here that after a perturbation, not only the number of individuals in the population may decrease (by increased mortality or dispersal) but also its social cohesion, likely reducing but also altering the information transfer within the social network composed by those individuals that remain in the site where perturbation occurs. Among other demographic processes, dispersal may alter social connections of both individuals remaining and those dispersing, with consequences for social network structure 45 . The perturbation regime suffered by the study population has likely triggered a social transition 58 in collective behaviour from philopatric to dispersal and with the fast diffusion of innovations such as the colonization of harbours by a large number of individuals, a habitat safe from predators never occupied before 54 . Previous studies have shown that responses of populations to perturbations may also depend on individual personalities in the population 22,59-61 . For example, dispersers are different from non-dispersing individuals for a suite of phenotypic traits, including their behavioural profile [62][63][64] . Heterogeneities in personalities for dispersal decision-making may have also played a role in our studied population, with most individuals dispersing to other sites after a period of disturbance, while some individuals remaining philopatric. This change may have also further consequences for social network stability, as performance in social groups may improve with heterogeneity in individual personalities 64,65 .
Our study opens new research avenues about the resilience of social organisms under perturbations; if perturbations affect social cohesion and heterogeneity in personalities in the population, we may wonder whether this population would be equally resilient to future perturbations. Additionally, in our study population, sociality seemed to operate not right after the first perturbation episode but after cumulative maintained perturbation 54 ;

Material and methods
Study species and study area. The Audouin's gull is a long-lived seabird with more than 80% of the global population breeding in the western Mediterranean (https ://www.iucnr edlis t.org/detai ls/22694 313/0 43 . The species was critically endangered until the early 1980s when it colonized a new site, the Punta de la Banya in the Ebro Delta (Fig. 2). Here, the large availability of both suitable breeding habitat and food resulted in rapid and exponential growth, ending with the site holding more than 70% of the total world population in 2006. The global population dynamics were mainly driven by this colony and after the exponential growth, the species was downgraded to a conservation category of "least concern" 67 . However, the Punta de la Banya colony is now collapsing and even if the species is colonizing new sites, the global population is decreasing at a 5% annual rate 43 .
In 1997, first carnivores (mainly foxes, but also badgers, beech martens, and least weasels) arrived at Punta de la Banya, likely due to their increasing densities in recent decades 68 and the attractiveness of the site in terms of food availability and lack of competition. Since then the site has been perturbed by the presence of carnivores. Annual censuses of breeding pairs at every patch within colonies at the Ebro Delta area have been carried out since colonization in 1981 to 2017 (Fig. 1, 3; Supplementary Table S2). In the Ebro Delta there are three colonies: Punta de la Banya colonized in 1981 and occupied throughout the study period, Sant Carles de la Ràpita harbour, occupied from 2011 to 2015, and Sant Antoni, occupied from 2013 up to now ( Fig. 3; Supplementary Figure S1). Within colonies, individuals are distributed in patches 69 . As patch location may change from 1 year to another, we annually geolocalized, mapped and defined the breeding patches ( Fig. 3; Supplementary Figure S2).
The species is monogamous, there is assortative mating by age, and from an evolutionary point of view is a bet-hedger, laying commonly 3 eggs, although few chicks survive, except in years with high food availability, when the strength of density-dependence is low 70 .
Individual data. During 1988-2017 a total of 30,290 individuals were captured and ringed as chicks at the Punta de la Banya 71,72 . From 2002 to 2017, resightings were made using spotting scopes from a distance all over the western Mediterranean with a total of 63,106 resights in the study area and 5593 different individuals resighted. Each year we recorded the breeding patch for each individual. To make sure that individuals were breeding and that they did so in a particular patch, we only selected those individuals seen during the breeding season in a particular patch showing unequivocal breeding behaviour: specifically, individuals making alarm calls, incubating eggs, or with chicks. Even if monitoring of the species started in 1988, resights made before 2002 could not be included in this analysis because we did not record breeding behaviour that assures us that the individual was breeding on that particular patch or just resting or prospecting there. After the selective filter of unequivocal breeders, our final database included 3548 individuals. www.nature.com/scientificreports/ We considered the breeding association between two individuals when individuals were recorded breeding in the same patch the same year, and not associated if they were breeding in different patches. If one of the individuals was missed or individuals do not occur in overlapping timeframes, were considered as missing data. SNA framework. Our social network, defined as the observed pattern of breeding association, was constructed taking individuals (N = 3548) as the nodes of the network and each edge dyad (i.e. two individuals) representing the fact that individuals breed in the same patch. We ended with a global sociomatrix, i.e. the matrix representation of the dyadic relationships among individuals, of 3548* 3548. Edges showed if two individuals bred in the same patch at least once in a certain period (see below). The network was not directional. Based on previous results on population dynamics, and the population size of this species and colony 43,54 , we divided our dataset into two main periods: a period defined as "stable phase" from 2002 to 2011, and a period of "transition phase to collapse", from 2012 to 2017 (Fig. 1).
We used the recently developed AME function from the AMEN package 73,74 , which can be applied to binary, ordinal, and continuous network data. This new approach is a random-effects regression model; uses an iterative Markov chain Monte Carlo (MCMC) algorithm that provides Bayesian inference of the parameters in the social relations regression model (SRM 75 ) using additive and multiplicative effects and combining the linear regression model with the covariance structure of the SRM 74 . The AME method, not currently used in research on animal social networks is also able to cope with missing and censored data, our data set complying with the assumption that individuals are missing at random. Through Markov chain Monte Carlo algorithm, the AME function approximates the joint posterior distribution of all unknown quantities in the system, including parameters and missing data. Specifically, as in almost all Bayesian modeling methods, at each iteration of the Markov chain, values for missing data are iteratively simulated from their conditional distribution, given the observed data values and the current values of the unknown model parameters. For that being accurate, "individuals should be missing at random", an assumption we feel confident with because we have no reason to suspect that those individuals that we fail to observe each year in the field are not a random sample of all the individuals. Coping with missing data is highly relevant when analysing sociality on wild populations, as detection rate for individuals is almost always imperfect, and properly controlling for missed observations is a very important step in social network analysis 76,77 . To create and visualize our networks we used the packages Amen 73 , Asnipe 78 , gdata 79 and igraph in R 80 .
Occurrence of social ties that persist over time. We investigated if individuals create social ties that persist over time longer than one breeding occasion by means of two approaches: (1) contingency tables and (2) the inclusion of time-dependent regression terms in the AME modelling framework 73 (see the previous section). We used both methods because this is the first application of the AME approach in an ecological context. We analysed data of the period of stability, from 2002 to 2011, dividing this period into two sub-periods of 5 years (2002-2006 and 2007-2011). In the contingency table approach, we tested if the probability of two individuals breeding in the same patch at least once during a 5-year period (2007-2011) was independent of having bred in the same patch at least once in the previous 5-year period (2002)(2003)(2004)(2005)(2006). We built a 3 × 3 table of frequencies, showing the frequencies of two individuals breeding or not together at least once during the second sub-period depending on whether they bred in the same patch or not at least once during the first sub-period, pulling apart those dyads with missing data. We then tested for deviation of random frequencies by χ 2 test.
With the SNA approach, we analysed the social ties between individuals using the AME function provided in the Amen R-package and including data from the first period (five previous years) as predictors of the association during the second period. We considered that this time window was not too large to include important death events, but large enough to account for the imperfect detection of individuals. To achieve convergence, we increased the number of iterations to 100,000 from the default value of 10,000 and lengthened the burn-in period to 500.

Impact of perturbations on social cohesion.
To assess if perturbations affected social structure in this species, we analysed as previously, with both the contingency table approach and the SNA approach, the social ties during the "transition phase to collapse" (2012-2017). To do so, we tested if the probability of breeding in the same patch in this phase (2012-2017) was independent of having bred in the same patch during five previous years (2007)(2008)(2009)(2010)(2011). We then compared these results from those previously observed during the "stability phase".
A potential concern was the reduced power during the collapse period because the number of individuals decreased from the stability period. In order to have similar power in both analyses, we performed the contingency table analysis during the stability period by drawing at random a number of observed associations equal to the number of observed associations during the collapse. We did this repeatedly (1000 times) and calculated which percentage of times the resulting chi-square was significant at the 5% level.
Regarding the SNA approach, it is advised 81 to do permutations on the raw data prior to the analysis and compare the result of some relevant statistic obtained with the original data to the distribution of the same statistic over the permutations. We chose the regression coefficient of the association of a dyad on the previous association of the same dyad as our statistic of social cohesion. This statistic is provided by the function ame of the package amen 73 . We calculated this statistic on the original data. Then, within each year, we rearranged randomly the individuals among the patches, keeping the same number of individuals within each patch. We did this 200 times and calculated each time the regression coefficient. Then, we situated the value of the regression coefficient from the original data among the distribution of regression coefficients from the permutated data and examined how extreme it was. If it was among the 5% extreme values, we concluded that there was indeed a non-random association of individuals within the patches.