Mass gain is an important fitness correlate for survival in highly seasonal species. While many physiological, genetic, life history, and environmental factors can influence mass gain, more recent work suggests the specific nature of an individual’s own social relationships also influences mass gain. However, less is known about consequences of social structure for individuals. We studied the association between social structure, quantified via social network analysis, and annual mass gain in yellow-bellied marmots (Marmota flaviventer). Social networks were constructed from 31,738 social interactions between 671 individuals in 125 social groups from 2002 to 2018. Using a refined dataset of 1,022 observations across 587 individuals in 81 social groups, we fitted linear mixed models to analyze the relationship between attributes of social structure and individual mass gain. We found that individuals residing in more connected and unbreakable social groups tended to gain proportionally less mass. However, these results were largely age dependent. Adults, who form the core of marmot social groups, residing in more spread apart networks had greater mass gain than those in tighter networks. Yearlings, involved in a majority of social interactions, and those who resided in socially homogeneous and stable groups had greater mass gain. These results show how the structure of the social group an individual resides in may have consequences for a key fitness correlate. But, importantly, this relationship was age-dependent.

Data collection

We examined the mass gain-social structure interface in the population of yellow-bellied marmots located in and around the Rocky Mountain Biological Laboratory (RMBL) in the Upper East River Valley, Gothic, Colorado (38°57’N, 106°59’W; ca. 2900 m elevation), which have been studied continuously since 1962. Yellow-bellied marmots are a facultatively social, harem-polygynous species of ground-dwelling squirrel that live in matrilineal colonies with one or two territorial males (Frase and Hoffman 1980; Armitage 1991). Throughout their five-month active season from early May to mid-September, marmots must accumulate enough fat reserves to survive hibernation. Thus, they allocate a large amount of time to foraging and digesting (Armitage et al. 1996). 

From 2002 to 2018, marmots were observed and repeatedly live trapped during their active season. Individuals were trapped using baited Tomahawk-live traps near burrow entrances, and immediately transferred to cloth handling bags to record body mass, sex, and age category (pups [< 1 year], yearlings [= 1 year], and adults [≥ 2 years]). Only adults are reproductively mature. All marmots were given two uniquely numbered permanent metal ear tags (Monel self-piercing fish tags #3, National Band and Tag, Newport, KY) and marked on their dorsal pelage with a nontoxic Nyanzol fur dye (Greenville Colorants, Jersey City, NJ, U.S.A.) to be identified from a distance. Virtually all marmots in our study population are trapped and marked annually, permitting us to accurately identify interacting individuals. For the few individuals that might lose their marks and not be recaptured after they moult, the fact that most other marmots at their respective colony site were marked, permitted us to identify them. 

Marmots were studied annually at the same sites in the Upper East River Valley. Colony sites can be grouped into a higher classification of up- and down-valley (5 are up-valley, 7 are down-valley). Up-valley is at a higher elevation and experiences harsher weather conditions than down-valley (Van Vuren and Armitage 1991; Blumstein et al. 2006; Maldonado-Chaparro et al. 2015b). Upon sexual maturity, nearly half of females and most males disperse (Armitage 1991) with most dispersal resulting in movement out of the valley. Social matrices were constructed for each group of connected individuals that appeared naturally within a valley location (up- or down-valley) and that did not interact with any other individuals.

Detailed social interactions in this population have been recorded since 2002. Behavioural observations were made during hours of peak activity (0700-1000 h and 1600-1900 h; Armitage 1962) using binoculars and spotting scopes from distances that did not disrupt normal social behaviour (20-150 m; Blumstein et al. 2009). For each interaction we then classified the behaviour as either affiliative (e.g., greeting, allogrooming, play) or agonistic (e.g., fighting, chasing, biting; detailed ethogram in Blumstein et al. 2009). We also recorded the initiator and recipient, time, and location of each interaction. Most interactions (81.69%) occurred between identified individuals. The initiator and/or recipient of 18.31% interactions could not be identified because the marmot’s dorsal fur mark was not visible because of either the marmot’s posture or visual obstructions. Excluding these interactions between unidentified individuals should not significantly influence social structure (Silk et al. 2015). Our social matrices only consisted of yearlings and adults because these cohorts were present early in the season, when social interactions are most common. We excluded pups from our matrices because most pups emerge in July and were therefore present only a fraction of the year. We also filtered out individuals with fewer than five interactions in a given location to eliminate those dispersing and that were not actually part of the social group (Wey and Blumstein 2012; Fuong et al. 2015; Yang et al. 2017; Blumstein et al. 2018).

Mass gain model

To predict body masses during the growing season, we calculated best linear unbiased predictions (BLUPs) by fitting linear mixed effects models from the repeated body mass recordings taken for yearlings and adults from 2002-2018 to predict 1 June and 15 August body mass. Models were fitted in R (R Development Core Team 2020; version 3.6.3) using the “lmer” package (Bates et al. 2015; version 1.1-23). Data used in our BLUPs consisted of 7,164 observations across 4,077 individuals and 56 years. There was a mean of 3.4 observations per individual (range: 1.0 – 24.0; Median = 2.0). Martin and Pelletier (2011) showed that BLUPs can make accurate body mass predictions when there are on average greater than three measurements of body mass per individual, a criterion our data meet. The repeatability of body mass in our models varies between 0.35 and 0.47 depending on the age-sex specific model. We acknowledge that using BLUPs in follow-up analysis (such as our linear models discussed below) can lead to higher rates of Type 1 error (Hadfield et al. 2010; Houslay and Wilson 2017). However, our large dataset used to produce the BLUPs helps to mitigate this error (Dingemanse et al. 2019). Additionally, we do not use mean and mode of the posterior distribution of each BLUP after fitting in lme4 (as proposed by Hadfield et al. 2010 and Houslay and Wilson 2017) because the estimates from lme4 are equivalent. We included individual identity, year, and site as random effects in the models, producing individual- and year-specific intercept predictions (Maldonado-Chaparro et al. 2015b; Kroeger et al. 2018; Heissenberger et al. 2020). Therefore, to calculate individual yearling and adult proportional mass gain, we divided individual body mass on 15 August by the body mass on 1 June.

Social network measures

Using social observation data collected from 2002 to 2018 and the R package “igraph” (Csardi and Nepusz 2006; version 1.2.5), we constructed weighted and directed social interaction matrices based on observed affiliative interactions between individuals for each year. These affiliative networks consisted of 31,738 social interactions between 671 individuals in 125 social groups. 13,668 of these interactions occurred down valley and 18,068 occurred up valley. From these matrices we calculated seven social network measures to quantify social structure (described in Table 1). Our observations of marmot social groups across their entire active season and low rate of unknown individuals involved in social interactions facilitates the reliability of the seven social network measures (mean across years per individual = 33.1, range = 8.9 to 91.3; Supplementary Table 1; Silk et al. 2015; Davis et al. 2018; Sánchez‐Tójar et al. 2018). Because some of the network measures could not be calculated for certain group sizes or group configurations (e.g., transitivity for a group of two or a linear group) and because mass gain rates could not be calculated for some individuals, we systematically removed all N/A’s for network measures across all models and removed individuals without mass gain rate values. This can be attributed to some individuals only being weighed once in a year, only observed a few times a year, or due to their membership in a small group (e.g., a group of two). This final 17-year dataset used in our analysis consisted of 1,022 annual observations of mass gain and group metrics, from 587 unique individuals, that lived in 81 different social groups (Supplementary Table 1).

Data Analysis

To test the relationships between social structure (quantified via seven network-level measures [Table 1]) and proportional mass gain, we fitted linear mixed models in R (R Development Core Team 2020; version 3.6.3) using the “lmer” package (Bates et al. 2015; version 1.1-23). Each model had a different network measure as the primary predictor variable; mass gain was the response variable across models. All models included group size (number of individuals in a social group), number of mass recordings, sex, age, and valley location as fixed effects. We included year, individual ID, and group ID as random effects (random effects are crossed as an individual ID may be seen in multiple years, and thus in multiple group IDs across years). Categorical variables sex, age, and valley location were mean-centered following Schielzeth (2010). As such, we coded females, yearlings, and down-valley individuals as ‘+1’ while males, adults, and up-valley individuals were coded as ‘-1’..

Group size was included as a fixed effect to account for network measures that may differ as a function of group size. The number of mass observations was included as a fixed effect to account for variation in the certainty of BLUP estimates. We included individual identity as a random effect to account for individuals that had observations over multiple years and colony ID to account for multiple members of the same group that shared a network measure within a given year. We also included interactions between the social network measure and sex, age, and valley location because these three variables are well-known correlates with mass gain; there are significant differences between the mass gain of females and males (Armitage 1998), yearlings and adults (Armitage et al. 1976), and individuals down-valley versus individuals up-valley (Van Vuren and Armitage 1991; Blumstein et al. 2006; Maldonado-Chaparro et al. 2015b). 

Proportional mass gain rate (response variable) and the number of mass observations (predictor variable) were log₁₀transformed and all variables then were standardized (mean-centered and divided by one SD using the base “scale” function in R; Becker et al. 1988). We checked for collinearity by calculating correlation coefficients between continuous predictors. No models had a correlation coefficient > 0.8 between the network measure and a fixed effect (Franke 2010; Shrestha 2020). After fitting each model, we calculated the marginal and conditional part R² values to estimate the variance explained by each of our fixed and random effects, using the “partR2” package in R (Stoffel et al. 2021; Nakagawa and Schielzeth 2013; version 0.9.1). The marginal part R² gives an estimate of the variance explained by each fixed effect whereas the conditional part R² gives an estimate of the variance explained by each fixed effect plus the variance explained by all the random effects. We estimated 95% confidence intervals for our part R² values using 100 parametric bootstrap iterations. Then we used the “check_model” function from the “performance” package in R (Lüdecke et al. 2020; version 0.6.1) to ensure each model met the assumptions of linear mixed models. Graphs were generated using “ggplot2” package in R (Wickham 2016; version 3.3.3).