Social network community structure and the contact-mediated sharing of commensal E. coli among captive rhesus macaques (Macaca mulatta)

Study location and subjects

The study was conducted at the California National Primate Research Center (CNPRC) and the School of Veterinary Medicine (SVM), University of California at Davis. Data were collected on 299 adult rhesus macaques (90 males, 209 females) between 3–29 years of age (mean = 7.7 years), across three social groups (Table 1). The groups were housed in separate, 0.2 ha outdoor enclosures. Animals were fed a standard diet of monkey chow twice per day at approximately 0,700 h and between 1,430 and 1,530 h. They were provided fresh fruit or vegetables once a week, with seed-mixture being provided daily. Water was available ad libitum, sporadically as natural puddles but mostly via artificial sources such as taps. For more information regarding the study groups, and intergroup differences in sociodemographic characteristics, see Table 1 or Balasubramaniam et al. (2016). The protocols used for this research were approved by the UC Davis Institutional Animal Care and Use Committee (IACUC; Protocol #: 18525; Office of Laboratory Animal Welfare (OLAW) Assurance Number: A3433-01), and were in accordance with the legal requirements of the jurisdictions in which the research was conducted.

Behavioral data collection

Behavioral and biological data were collected during a 6-week sampling period per group, with two groups being observed in the spring (Group I: March–April 2013; Group III: March–April 2014) and one being was observed in the fall (Group II: September–October 2014). For each group, three observers collected data for 6 h on 4 days per week from 0,900–1,200 h and 1,300–1,600 h. Observers used an Event Sampling design to record both mild and severe aggressive interactions, and Scan Sampling to record affiliative grooming and huddling interactions (Altmann, 1974). Further details regarding the precise definitions of behaviors and sub-categories (for aggression) may be found in Balasubramaniam et al. (2016).

The event sampling approach has been previously proposed as being useful to optimize reliable data collection in large social groups to improve statistical power, and navigate non-independence issues that may affect the computation of social network measures (Balasubramaniam et al., 2016; Farine & Whitehead, 2015; McCowan et al., 2011; Vandeleest et al., 2016). Further, our frequency of scan sampling of affiliative interactions—once every 20 min during a six-hour duration of sampling per day—was also intense, generating approximately 432 scans in total. This scan sampling regime was sufficient to generate biologically meaningful social networks; a recent study on wild Japanese macaques revealed that for a given duration, frequent instantaneous scan sampling generates identical amounts of grooming data to focal sampling (Romano et al., 2016). Nonetheless, to verify that our sampling effort was adequate, we first computed the values of three social network measures—specifically Newman’s eigenvector-based modularity, mean degree, and network density (see Table S1 for definitions and R packages used)—computed from each social network used in the analyses. We then compared these observed measures to the corresponding measures calculated from 1,000 permuted networks generated by bootstrapping increasingly smaller subsets (100%–10%, with decrements of 10%) of edges from the original network (Croft et al., 2011; Farine & Whitehead, 2015; Lusseau, Whitehead & Gero, 2008). Plots of network measures vs. percentage sampling effort revealed asymptotic trends, indicating that sufficient sampling effort had been reached (e.g., Newman’s modularity: Fig. 1A–1C; Table S2).

Bacterial isolations and DNA fingerprinting

All individuals within a particular group were sampled on the same day; this was critical to ensure comparability, since there could be significant turnover of E. coli genotypes in the mammalian gut (Anderson, Whitlock & Harwood, 2006). Further, sampling was conducted on a pre-selected day on the final week of the behavioral observation phase, in order to facilitate the detection of contact-based sharing of E. coli attributable to the animals’ recent history of interactions. Prior to fecal collection, each animal was immobilized (10 mg/kg of ketamine) and given standard physical examinations by veterinary staff (e.g., checked for injuries, weighed). Two fresh fecal swabs were collected from every macaque at the end of the behavioral observation period, following previously published methods (Good, May & Kawatomari, 1969). A sterile cotton-tip swab was inserted into the rectum of each individual, rotated gently to collect fecal material, and immediately immersed into a 15 ml test-tube (labeled with the animal ID containing a sterile growth medium (Group I: Phosphate Buffer Saline (PBS); Groups II and III: Tryptic Soy Broth (TSB; BD, Franklin Lakes, NJ, USA))); a duplicate sample was taken using another sterile swab and placed into a second tube. The samples were incubated within 4 h of collection, with orbital rotation of 100 rpm at (1) 25 °C for 2 h, (2) 42 °C for 8 h, and (3) held static at 6°C overnight. Commensal E. coli was isolated from the TSB enrichment. First, 10 uL of the enrichment was streaked for isolation onto MacConkey agar plates (BD, Franklin Lakes, NJ, USA) and incubated at 37 °C for 18–24 h. From these, suspect colonies were streaked on to Eosine Methylene Blue agar (EMB) (BD, Franklin Lakes, NJ, USA) and incubated under similar conditions for a 24-hour period. Following an additional cycle of isolation, streaking and incubation on MacConkey plates, all suspect E. coli isolates were biochemically confirmed using Triple Sugar Iron (TSI) (Remel, Lenexa, KS, USA), Citrate (Remel, Lenexa, KS, USA), and Urea (BD, Franklin Lakes, NJ, USA), Methyl Red-Voges-Proskauer (MR-VP) (BD, Franklin Lakes, NJ, USA), and Indole (BD, Franklin Lakes, NJ, USA). These tests confirmed commensal E. coli from the majority of individuals sampled within each group (Table 1). Remaining individuals were deemed ‘untypable’, and were excluded from further analyses. Confirmed isolates were then banked and frozen within a −80 °C freezer for subsequent bacterial fingerprinting and phylogenetic reconstruction.

We used PFGE (PulseNet Pulsed Field Gel Electrophoresis, using the CDC protocol) to generate DNA fingerprint profiles for symbiotic E. coli. This technique is a well-established and valid method for a surface comparison across isolates (Cesaris et al., 2007; Kilonzo et al., 2011; Kondo et al., 2010; Ribot et al., 2006). While more novel techniques, such as whole genome sequencing, can be used to reconstruct the phylogenetic relationships of organisms using nucleotide datasets, we preferred PFGE owing to its sufficiency in assessing genus-typical bacterial diversity. Further, given (a) the large number of samples (299 macaques), and (b) our focus on a specific, non-pathogenic inhabitant of the gut microbiome (commensal E. coli), implementing metagenomic processing to develop complete profiles of gut microbiota was beyond the scope of this study. Finally, PFGE has been shown to perform well in previous research that links microbial sharing with the spatial and social contact networks of African ungulates (VanderWaal et al., 2013a; VanderWaal et al., 2014a).

From a single banked isolate from each individual macaque, we streaked bacterial colonies onto Tryptic Soy Agar (TSA), which in turn facilitated the production of agarose plugs containing the lysed bacterial colonies. Following this, bacterial DNA was digested with Xba-1 restriction enzyme, loaded and run through an agarose gel, and stained with ethidium bromide to visualize DNA banding patterns. DNA fingerprints were then grouped and standardized using the Bionumerics software (version 6.6, Applied Maths, Inc). Following this, we used a bandmatching procedure to extract a binary, bipartite matrix of bacterial genotypes (Table S3). Rows represented individual monkeys, and the columns each of 66 band positions identified and optimized using automated analytical parameters to fit the dataset. Cells in this matrix indicated either a presence (1) or absence (0) of a band for each individual in each band position. We also reconstructed three phylogenetic trees of bacterial similarity, one for each of the three study groups (e.g., Group I: Fig. 1). For this, we used the UPGMA (Unweighted Paired Group Method with Arithmetic-mean) procedure, which provides reliable topologies via bandmatching of fingerprint data. From each group-specific tree, we extracted a “similarity” matrix of pairwise cophenetic coefficients, i.e., the correlation between the similarity in densitometric curves and the phylogenetic branch length distance between each pair of macaques.

Bacterial isolations from environmentally deposited feces

Our study system and design were particularly well-suited for examining evidence for microbial sharing via socially-mediated fecal-oral contact-routes. First, there was little variation in climatic conditions at CNPRC within and across the observation periods. The monkeys experienced generally hot and dry conditions, with extended periods of sunlight. Although some strains of E. coli may generally persist in the environment for between eight-25 weeks, hot and dry environmental conditions are considered highly unfavorable for the survival of gut E. coli deposited into the environment (Habteselassie et al., 2008; Sinton, Hall & Braithwaite, 2007; Van Elsas et al., 2011). Further, the water supply from the artificial taps tested negative for commensal E. coli. For these reasons, the sharing of E. coli on account of environmental space-use overlap, and/or from using the same artificial water-sources was highly unlikely. Nonetheless, to rule out E. coli sharing on account of independent acquisition from environmental feces, we compared, for one of our study groups, the E. coli strains isolated and fingerprinted from macaque rectal swabs to isolates from environmentally deposited feces. Specifically, we collected 18 environmental fecal samples, three from within each of six, equally divided sections of the enclosure, on the same day of macaque rectal swab collection. To collect each sample, we used a sterile sponge held with a pair of forceps to wipe a roughly 0.5 m ×0.5 m area on a man-made enrichment surface (e.g., perching frames, rails) that contained freshly deposited macaque feces. The sponges were then immediately immersed in sterile, pre-labeled bags containing TSB, which were then incubated and processed simultaneously with the macaque rectal swabs for E. coli confirmation and fingerprinting, using the procedures described above.

Population genetics, social network, and statistical analyses

To examine whether the overall population genetic diversity of E. coli was influenced by group membership controlling for other individual attributes, we performed a series of Analysis of Molecular Variance (AMOVA: Excoffier, Smouse & Quattro, 1992) tests with 1,000 random permutations of the data, using the pegas package in R (Paradis, 2010). In each test, the “outcome” was a square matrix of Euclidean distances of E. coli genetic similarity between each pair of individuals, calculated from the bipartite matrix of band presence-absence extracted using the bandmatching procedure (Table S3). On these matrices, we ran three hierarchical or nested AMOVAs. With group membership (I, II, and III) used as a level-two variable, we ran models with (T1) species-typical age-class (old (>13 years), prime (between 4 and 13 years), (T2) sex-class (males, females), and (T3) matriline membership (individuals within the same matriline were grouped together) each nested within groups. This last analysis was run for just groups I and II, since the matriline structure for Group III was highly fragmented on account of this group being primarily composed of younger individuals introduced from various other groups (Table 1; Balasubramaniam et al., 2016). To examine whether the overall E. coli genetic similarity among individuals was different across the three groups, we first calculated the row-wise mean % similarity of each individual’s E. coli genotype with its group members, from the similarity matrices extracted from the group-specific phylogenetic trees. We then ran a one-way ANOVA with a Tukey HSD posthoc test for multiple comparisons, with individuals’ mean similarity coefficient set as the continuous variable, and “group ID” set as the fixed factor (Levene’s homogeneity of variance test: F = 2.12, df = 2, p = 0.12).

To examine whether the frequencies of dyadic social behavioral interactions influenced the % similarity in E. coli genotypes, we ran multiple, univariate Multiple Regression Quadratic Assignment Procedure (or MR-QAP) models with double dekker semi-partialling and 1,000 permutations (Dekker, Krackhardt & Snijders, 2007; Hubert, 1987; Krackhardt, 1987). MR-QAP accounts for the non-independence of dyadic datasets (Hanneman & Riddle, 2005) by coercing matrices into vectors. After performing a standard linear regression across the corresponding cells of a dependent matrix and one or more independent or co-variate matrices, the procedure uses a Monte Carlo method to randomly permute the rows and columns of the dependent matrix. It thus re-computes regression coefficients 1,000 times to generate a distribution of coefficients against which the observed coefficients may be compared. The semi-partialing approach has been shown to be fairly robust to the distribution (normal, gamma, negative binomial) of values in the outcome matrices (Dekker, Krackhardt & Snijders, 2007). Thus, although the distribution of E. coli % similarity deviated significantly from normality for all three groups (e.g., Group I: Shapiro–Wilcoxon test: w = 0.96, p < 0.01), we still ran and interpreted linear MR-QAP matrix regressions. For each macaque group, we ran four univariate models. In each model, the dependent network was the “similarity” matrix of cophenetic correlation coefficient matrix of % similarity in E. coli. We used the netlm linear function in the SNA R package (Butts, 2008). Independent networks included frequencies of (m1) grooming, (m2) huddling, and (m3) aggression. We included both mild and severe aggression in the aggression networks; although mild aggression does not involve direct contact, it was included since it may be linked to the likelihood of occurrence of moderate and severe aggression in captively housed macaques that interact frequently (Balasubramaniam et al., 2016). Since our social networks were not collinear (range of Pearson’s (r) row-wise matrix correlation coefficients for all pairs of networks and all groups: 0.02 < r < 0.68), we also ran a single multivariate model per group that included all three types of networks as predictors of E. coli similarity (Table S4). Previous studies have suggested that a major component of nonhuman primate gut microbiota may be evolutionarily conserved, or inherited (Ley et al., 2008; McCord et al., 2014). Further, social interactions among macaques may be more likely to occur among closely related maternal kin (Berman, 2011; Chapais, 2006). For these reasons, we also regressed a (m4) binary matrix of kinship (1: close kin dyads with a relationship coefficient (r) of ≥0.125; 0: distant kin or unrelated individuals (r < 0.125)), on E. coli genotypic similarity. We ran the kinship analyses only for Groups I and II since the matriline structure for Group III was highly fragmented (see above).

Despite its utility in handling interdependencies in the data, the MR-QAP method has several restrictions. In contrast to the ordinary least squares method, it is not possible to calculate degrees of freedom, statistical power, or effect sizes in MR-QAP regression (Ferrin, Dirks & Shah, 2006). R-squared values also tend to have little meaning (Gibbons, 2004; Zagenczyk et al., 2013). Further, there is controversy in the use of goodness-of-fit statistics (e.g., AIC) for MR-QAP. So rather than likelihood-based model selection criteria,we interpreted all results using the β coefficients, and the p values computed based on permutation tests which is the primary statistic of interest in MR-QAP analyses (Zagenczyk et al., 2013).

To determine whether E. coli sharing is more readily detectable across sets of closely interacting or spatially associated individuals, we used the Data Cloud Geometry (or DCG) method to reconstruct behavioral community structures (Fushing & McAssey, 2010; Fushing et al., 2013; McCowan et al., 2016). DCG identifies network community structure at multiple levels by performing a random walk through an empirical network guided by the data. Cumulatively, these random walks produce a similarity matrix describing the pairwise similarity in social connections, from which a hierarchical tree of clustering is generated (Fushing & McAssey, 2010; Fushing et al., 2013). From a biological perspective, a DCG cluster may be therefore defined as a subset of group members whose social ties are both closer to, and stronger among each other, than they are to other group members (VanderWaal et al., 2014b). Such “closer” individuals, whether close kin (Berman, 2011), non-kin allies (Seil et al., 2017), and/or strong social bond investors (Silk, Alberts & Altmann, 2003), tend to cluster together at a lower level of this tree than individuals with fewer connections and/or similarities in connections. We used DCG because it offers specific advantages over other commonly used methods to identify cluster structure, such as Hierarchical Clustering (Corpet, 1988; Johnson, 1967) (summarized in Fushing et al., 2013; VanderWaal et al., 2014b). First, it does not require dyadic relationships to be binary and instead, utilizes the strength of relationships (e.g., frequencies of behavioral interactions). Second, in comparison with Hierarchical Clustering trees, DCG trees are more robust, less sensitive to measurement errors, and provide information on the intrinsic scales embedded within the data cloud. This is because of the implementation of stricter rules for assigning nodes to the same cluster (ultrametric or strong triangle inequality rule), making DCG more accurate in identifying cluster structure (Fushing et al., 2013; VanderWaal et al., 2014b).

We constructed three DCG trees per macaque group, one each from dyadic grooming, huddling, and aggression. We used Monte Carlo tests to determine whether the observed cluster membership of macaques was significant. At each hierarchical level for each tree, we generated 1,000 random clustering configurations. We then compared the observed mean behavioral frequency within a cluster, to a distribution of mean frequencies of behaviors from the randomly generated cluster configurations. We thus considered a particular level to be significant if this mean fell within the 95th percentile of the permuted distribution (as in VanderWaal et al., 2014b). To determine whether % similarity in E. coli was greater among individuals within the same, versus across different DCG communities, we used a series of Wilcoxon rank-sum tests. We favored a non-parametric test since the dyadic E. coli similarity coefficients were deviated significantly from normality (see above). Nevertheless, we also ran a series of randomization tests, which compared the observed mean within-cluster E. coli similarity coefficients to a distribution of coefficients generated from 1,000 permuted datasets in which individuals’ behavioral cluster membership was assigned randomly. Finally, to rule out E. coli sharing on account of macaques’ shared exposure to environmental feces, we used an additional Wilcoxon rank-sum test to compare the mean % similarity of macaque-macaque E. coli isolates to the % similarity of macaque-environmental fecal isolates for Group II. All statistical analyses were performed using R (ver 3.1.3), with the value of α being set after a Bonferroni correction for some relevant tests (AMOVA tests: α = 0.02; MR-QAP tests: α = 0.01; Wilcoxon rank-sum tests: α = 0.02.

Intergroup variation in bacterial genotypic diversity

From 299 individual macaques, we isolated, confirmed, and generated a fingerprint profile for 243 individuals, 79 in Group I, 78 in Group II, and 86 in Group III. Table 2 shows the results from the AMOVA tests. As predicted, group membership was responsible for the highest proportion of the observed genotypic variance, despite the nesting of age-category, sex-category, and matriline membership within groups (Table 2). In comparison to group membership, these nested variables had little effect on E. coli genetic variance. Further, permutation tests associated with the AMOVAs showed that the genetic variance across different groups was significantly greater than chance datasets in which group membership was assigned randomly to the isolates. In contrast, the variance across age-categories, sex-categories and matrilines nested within groups failed to reach significance.

We found significant differences in the mean % similarity of E. coli of individuals across groups (one-way ANOVA: F_2,243 = 52.8, p < 0.01; Fig. 2). A post-hoc Tukey test revealed that the mean similarity was significantly higher among the macaques of Group III (M = 50.31, SD = 4.33), in comparison to macaques in Groups I (M = 45.39, SD = 5.08), and II (M = 43.53, SD = 3.68) (Fig. 2).

Effects of dyadic social behavioral interactions and kinship on E. coli similarity

Contrary to our predictions, MR-QAP models showed no clear associations between pairwise E. coli % similarity coefficients and the dyadic frequencies of social behavioral interactions, and/or kinship (Table 3). Such a lack of association was largely consistent across all three types of contact—grooming, huddling, and aggression—frequencies. Although huddling and aggression frequencies both showed a positive impact on E. coli similarity in Group II, the effects were not significant after correcting for multiple comparisons (Table 3). The directions of coefficients were highly inconsistent across the co-variate matrices. For instance, although both huddling and aggression appeared to have positive associations with E. coli similarity for Group II, grooming showed a negative coefficient. Further, E. coli similarity was also not higher among close-kin dyads. Finally, such inconsistencies between dyadic social behaviors and E. coli % similarity were also a feature of the results from multivariate models that combined all three types of social networks (Table S4).

Bacterial sharing among macaque social network communities

Reconstructions of social network community structures using the DCG approach revealed multiple hierarchical levels of clustering, in which individuals were embedded in higher-order communities and sub-communities. The characteristics of the DCG tree for the different study groups and types of behavioral networks are summarized in Table 4. Figure 3 shows the trees for Group I. We identified between three and five hierarchical levels of clustering in all the trees. Permutation tests run on each tree and at each hierarchical level established that the observed cluster memberships at all these levels were significantly different compared to membership within 1,000 randomly-generated communities (p < 0.01). We hence resorted to assigning communities based on cluster membership at the intermediate hierarchical level of each tree (e.g., level-2 in a 3-level tree, level-3 in a 5-level tree), to ensure both optimum sizes and numbers of communities (Table 4).

Wilcoxon rank-sum tests established that the mean similarity in E. coli among individuals within the same behavioral communities were consistently and significantly greater than the similarity in E. coli among individuals between different communities (Figs. 4A–4C; Table 5). Group III showed the greatest differences, as indicated by the highest z coefficients for all three types of behavioral communities. In Groups I and II, the z coefficients were higher for certain types of communities, specifically huddling communities in Group I and grooming in Group II. Further, randomization tests strongly supported our findings from the Wilcoxon rank-sum tests. These showed the greatest, most consistent support for contact-mediated community-wide sharing among Group III, in comparison to groups I and II in which the extent of support was contingent on the type of behavioral community. For Group III, the degree of within-community E. coli similarity was significantly greater than expected by chance for all three types of communities. In comparison, the results were less consistent for Group I (only the huddling network reached significance) and Group II (grooming and aggression networks, but not huddling, reached significance) (summarized in Table 5).

Out of the 18 environmental fecal samples that we processed for Group II, we confirmed and generated an E. coli genotypic profile from 15 samples. A Wilcoxon rank-sum test revealed that isolates from fecal swabs of individual macaques were more similar to each other than they were to these environmental fecal isolates (Z = 98.5, p < 0.01). This confirmed that E. coli sharing among captive macaques is more likely to occur via social contact-mediated routes than via independent acquisition from environmental feces.