Genetic Variation in the Social Environment
Contributes to Health and Disease

Amelie Baud, Megan K. Mulligan, Francesco Paolo Casale, Jesse F. Ingels, Casey J. Bohl, Jacques Callebert, Jean-Marie Launay, Jon Krohn, Andres Legarra, Robert W. Williams, et al.

► To cite this version:

Amelie Baud, Megan K. Mulligan, Francesco Paolo Casale, Jesse F. Ingels, Casey J. Bohl, et al.  
Genetic Variation in the Social Environment Contributes to Health and Disease. PLoS Genetics,  
2017, 13 (1), Non paginé. 10.1371/journal.pgen.1006498 . hal-01603345

HAL Id: hal-01603345

https://hal.science/hal-01603345

Submitted on 26 May 2020

**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>
RESEARCH ARTICLE

# Genetic Variation in the Social Environment
Contributes to Health and Disease

Amelie Baud¹*, Megan K. Mulligan²*, Francesco Paolo Casale¹*, Jesse F. Ingels², Casey J. Bohl², Jacques Callebert³, Jean-Marie Launay³, Jon Krohn⁴, Andres Legarra⁵, Robert W. Williams², Oliver Stegle¹*

1 European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, United Kingdom, 2 Department of Genetics, Genomics and Informatics, University of Tennessee Health Science Center, Memphis, Tennessee, United States of America, 3 AP-HP, Hôpital Lariboisière, Department of Biochemistry, INSERM U942, Paris, France, 4 Wellcome Trust Centre for Human Genetics, Oxford, United Kingdom, 5 INRA, UMR 1388 GenPhySE, Castanet Tolosan, France

© These authors contributed equally to this work.

* abaud@ebi.ac.uk (AB); oliver.stegle@ebi.ac.uk (OS)

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

## Abstract

Assessing the impact of the social environment on health and disease is challenging. As social effects are in part determined by the genetic makeup of social partners, they can be studied from associations between genotypes of one individual and phenotype of another (social genetic effects, SGE, also called indirect genetic effects). For the first time we quantified the contribution of SGE to more than 100 organismal phenotypes and genome-wide gene expression measured in laboratory mice. We find that genetic variation in cage mates (i.e. SGE) contributes to variation in organismal and molecular measures related to anxiety, wound healing, immune function, and body weight. Social genetic effects explained up to 29% of phenotypic variance, and for several traits their contribution exceeded that of direct genetic effects (effects of an individual's genotypes on its own phenotype). Importantly, we show that ignoring SGE can severely bias estimates of direct genetic effects (heritability). Thus SGE may be an important source of "missing heritability" in studies of complex traits in human populations. In summary, our study uncovers an important contribution of the social environment to phenotypic variation, sets the basis for using SGE to dissect social effects, and identifies an opportunity to improve studies of direct genetic effects.

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

OPEN ACCESS

**Citation:** Baud A, Mulligan MK, Casale FP, Ingels JF, Bohl CJ, Callebert J, et al. (2017) Genetic Variation in the Social Environment Contributes to Health and Disease. PLoS Genet 13(1): e1006498. doi:10.1371/journal.pgen.1006498

**Editor:** Marcus W. Feldman, Stanford University, UNITED STATES

**Received:** March 2, 2016

**Accepted:** November 21, 2016

**Published:** January 25, 2017

**Copyright:** © 2017 Baud et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** Gene expression data from the experiment with inbred strains are available from ArrayExpress E-MTAB-5276. Phenotype data for the same experiment are provided as S7 Table.

**Funding:** The High-Throughput Genomics Group at the Wellcome Trust Centre for Human Genetics is funded by Wellcome Trust grant reference 090532/ Z/09/Z and MRC Hub grant G0900747 91070. AB was supported by fellowships from the EMBL Interdisciplinary Postdoc Programme under Marie Curie COFUND Actions and the Wellcome Trust

## Author Summary

Daily interactions between individuals can influence their health both in positive and negative ways. Often the mechanisms mediating social effects are unknown, so current approaches to study social effects are limited to a few phenotypes for which the mediating mechanisms are known *a priori* or suspected. Here we propose to leverage the fact that most traits are genetically controlled to investigate the influence of the social environment. To do so, we study associations between genotypes of one individual and phenotype of another individual (social genetic effects, SGE, also called indirect genetic effects). Importantly, SGE can be studied even when the traits that mediate the influence of the
(105941/Z/14/Z). MKM, JFI, CJB, and RWW are supported in part by NIAAA grants (U01 AA016662, U01 AA013499, U01 AA014425) and the UTHSC Center for Integrative and Translational Genomics. The funders played no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing Interests:** The authors have declared that no competing interests exist.

social environment are not known. For the first time we quantified the contribution of SGE to more than 100 organismal phenotypes and genome-wide gene expression measured in laboratory mice. We find that genetic variation in cage mates (i.e. SGE) explains up to 29% of the variation in anxiety, wound healing, immune function, and body weight. Hence our study uncovers an unexpectedly large influence of the social environment. Additionally, we show that ignoring SGE can severely bias estimates of direct genetic effects (effects of an individual's genotypes on its own phenotype), which has important implications for the study of the genetic basis of complex traits.

## Introduction

Social interactions contribute to health and disease (e.g. peer smoking increases one's risk of taking up smoking). So far, quantifying social effects has required a clear hypothesis about the mechanisms mediating the influence of the social environment (in the example above peer smoking is the trait that mediates the social influence). For many phenotypes however, such hypotheses do not exist. Therefore we propose an alternative strategy to study social effects: we investigate effects on an individual's phenotype that arise from genotypes of social partners (social genetic effects, SGE, also called indirect genetic effects[1, 2]). SGE constitute the genetic basis of social effects and can be detected without prior knowledge of the phenotypes through which the social influence is exerted. SGE have been reported for interactions between mothers and offspring (maternal genotypes indirectly affect offspring phenotypes)[3-8] and more recently for interactions between adult individuals, in livestock and wild animals[2, 9-17] For example, growth rate in farm pigs has been found to be in part determined by the genetic makeup of the other pigs in the pen[2]. However, the extent to which SGE explain variation in biomedical traits is largely unknown.

If SGE do contribute to such traits, they are a promising approach to quantify effects of the social environment. Additionally, they provide an anchor to investigate causal paths and dissect the mechanisms underlying social effects. Finally, in studies of direct genetic effects carried out by the broad community, SGE may be used to account for social environmental effects.

Our study aimed at quantifying the contribution of SGE to multiple biomedical traits. We uncover unexpectedly large social genetic effects on multiple organismal and molecular phenotypes.

## Results

To investigate whether SGE explain variation in biomedical traits, we considered two experiments in laboratory mice involving complementary genetic designs, and assessed both organismal and gene expression traits.

### Experiment with two inbred strains

We first carried out an experiment with two inbred strains. We chose C57BL/6J (B6) and DBA/2J (D2), the progenitor strains of the largest mouse reference population, the BxD recombinant inbred panel[18]. We co-housed 86 mice at weaning as B6/B6, B6/D2 or D2/D2 pairs. After a period of six weeks during which the mice interacted undisturbed in their home cages, we collected 50 organismal phenotypes relevant to unconditioned anxiety, helplessness (a measure of depressed mood), general locomotor activity, stress, social dominance, wound
a

<table><thead><tr><th>Experimental test</th><th>Number of measures</th><th>Measures of</th><th>Age (weeks)</th></tr></thead><tbody><tr><td>Elevated plus maze (EPM)</td><td>4</td><td>Anxiety and locomotor activity</td><td>9</td></tr><tr><td>Activity chamber (AC)</td><td>14</td><td>Anxiety and locomotor activity</td><td>9</td></tr><tr><td>Tail suspension test (TST)</td><td>5</td><td>Helplessness (depression)</td><td>10</td></tr><tr><td>Forced swim test (FST)</td><td>5</td><td>Helplessness (depression)</td><td>10</td></tr><tr><td>Body weight</td><td>9</td><td>Body weight</td><td>3-10</td></tr><tr><td>Adrenal gland weights</td><td>2</td><td>Stress</td><td>11</td></tr><tr><td>Gonad and preputial gland weights</td><td>2</td><td>Social dominance</td><td>11</td></tr><tr><td>Brain noradrenaline and dopamine</td><td>5</td><td>Stress</td><td>11</td></tr><tr><td>Brain and blood serotonin</td><td>3</td><td>Social dominance</td><td>11</td></tr><tr><td>Ear punch</td><td>1</td><td>Wound healing</td><td>11</td></tr><tr><td>Gene expression in the prefrontal cortex</td><td>12,898</td><td></td><td>11</td></tr></tbody></table>

b

<figure><img src="image_4.png" /><figcaption>**Fig 1. Experiment with two inbred strains.** (A) Phenotypes collected in this experiment. Organismal phenotypes were measured on each of the 86 mice; after quality control, gene expression data was available for 79 mice (Methods). (B) Example phenotype (Wound area, from an ear punch, a measure of the rate wound healing) showing that phenotypic variation arises from variation in the strain of the focal (i.e. phenotyped) mouse (direct genetic effects, DGE), variation in the strain of its cage mate (social genetic effects, SGE), and an interaction between the two whereby the effect of the strain of the cage mate depends on the strain of the focal mouse. B6: C57BL/6J, D2: DBA/2J.</figcaption></figure>

doi:10.1371/journal.pgen.1006498.g001

healing and body weight (Fig 1A, S1 Table). We also profiled genome-wide gene expression in the prefrontal cortex (PFC). The PFC was selected because it is involved in coordinating behavioural responses based on sensorimotor information, motivation and affect, all of which may be affected by the social environment; as a result individual differences in PFC expression levels may reflect behavioural responses to the social environment [19–21].

In this design with two inbred strains there are three potential genetic sources of phenotypic variation: differences in the strain of the focal (i.e. phenotyped) mice (DGE), differences in the strain of their cage mates (SGE) and an interaction between the two, whereby the effect of the strain of the cage mate depends on the strain of the focal mouse (Fig 1B, S1 Fig and Methods). Variance partitioning and model selection (see Methods) provided evidence that interactions between DGE and SGE are common across traits (S2 Fig), with SGE typically affecting a specific trait in one strain (B6 or D2) but not the other (e.g. Fig 1B, S3 Fig). Thus, for each trait, we modeled the measurements collected in B6 mice and D2 focal mice separately (S1 Fig).

For eleven out of 50 organismal phenotypes we found significant SGE in either B6 or D2 mice (P < 0.05, FDR < 43%—see Methods; Table 1, S2 and S3 Tables). Strain-specific differences existed as different measures were affected by SGE in B6 and D2 mice. Of these eleven phenotypes, two were measures of stress and six were measures of anxiety, providing evidence that variation in the genetic makeup of cage mates causes variation in stress-related phenotypes. The direction of effect was consistent across the six measures of anxiety (S3 Fig). We also detected strong SGE on the rate of wound healing (measured from an ear punch, P = 6.6 $10^{-3}$, Q = 0.36), showing that SGE are not limited to behaviors. (Table 1). Importantly, SGE
<table><caption>**Table 1. Organismal phenotypes significantly affected by SGE in the experiment with inbred strains (P < 0.05).** B6 and D2 focal mice were analysed separately. SGE Q value: see Methods. SGE variance (%): Percentage of phenotypic variance explained by SGE +/- standard error (see Methods). * the expression in the PFC and in B6 focal mice of genes associated with the regulation of dopamine metabolism is also affected by SGE (Table 2).</caption><thead><tr><th>Focal mice</th><th>Phenotype</th><th>Measure of</th><th>SGE P value</th><th>SGE Q value</th><th>SGE variance (%)</th></tr></thead><tbody><tr><td>B6</td><td>Wound area</td><td>Wound healing</td><td>6.6 10<sup>-3</sup></td><td>0.36</td><td>18 +/- 12</td></tr><tr><td>B6</td><td>Midbrain concentration of noradrenaline</td><td>Stress</td><td>7.2 10<sup>-3</sup></td><td>0.36</td><td>15 +/- 10</td></tr><tr><td>B6</td><td>Midbrain concentration of dopamine*</td><td>Stress</td><td>2.6 10<sup>-2</sup></td><td>0.38</td><td>10 +/- 9</td></tr><tr><td>D2</td><td>Time spent not moving in center first min (AC)</td><td>Anxiety</td><td>1.1 10<sup>-2</sup></td><td>0.38</td><td>15 +/- 11</td></tr><tr><td>D2</td><td>Total arm entries (EPM)</td><td>Locomotor activity</td><td>1.6 10<sup>-2</sup></td><td>0.38</td><td>12 +/- 10</td></tr><tr><td>D2</td><td>Ratio time spent not moving periphery/center min 2 to 5 (AC)</td><td>Anxiety</td><td>2.6 10<sup>-2</sup></td><td>0.38</td><td>14 +/- 11</td></tr><tr><td>D2</td><td>Ratio time spent not moving periphery/center min 6 to 30 (AC)</td><td>Anxiety</td><td>3.2 10<sup>-2</sup></td><td>0.38</td><td>11 +/- 9</td></tr><tr><td>D2</td><td>Number immobility bouts first 2 min (FST)</td><td>Helplessness</td><td>3.4 10<sup>-2</sup></td><td>0.38</td><td>8 +/- 8</td></tr><tr><td>D2</td><td>Ratio time spent periphery/center min 2 to 5 (AC)</td><td>Anxiety</td><td>3.5 10<sup>-2</sup></td><td>0.38</td><td>12 +/- 10</td></tr><tr><td>D2</td><td>Ratio ambulatory time periphery/center min 6 to 30 (AC)</td><td>Anxiety</td><td>3.8 10<sup>-2</sup></td><td>0.38</td><td>9 +/- 8</td></tr><tr><td>D2</td><td>Proportion of entries in open arms (EPM)</td><td>Anxiety</td><td>4.8 10<sup>-2</sup></td><td>0.43</td><td>10 +/- 9</td></tr></tbody></table>

doi:10.1371/journal.pgen.1006498.t001

explained a considerable proportion of phenotypic variance (up to 18%), showing that social effects of genetic origin are important contributors to phenotypic variation in these two strains.

We also found evidence that SGE affect gene expression in the PFC. A gene set enrichment analysis based on the contribution of SGE to gene expression levels (see Methods) revealed “integrin-mediated signaling pathway” (P = 2.5 10⁻⁶, Q = 0.024) and “regulation of dopamine metabolic process” (P = 7.0 10⁻⁴, Q = 1) as the most significant Gene Ontology (GO) terms in D2 and B6 mice respectively (Table 2). The latter enrichment, although statistically less robust, is strikingly consistent with our finding that SGE affect dopamine levels measured by HPLC in B6 mice (Table 1). Our results therefore converge to show that variation in the genetic makeup of cage mates causes variation in behavioral, biochemical, and gene expression traits relevant to stress. Although the large number of tests limits the statistical power to detect SGE on individual genes (12,898 genes tested), we identified three genes with significant effects (Srsf2, Phlpp2, and Ppid; FDR < 33%; S4 Table).

<table><caption>**Table 2. Gene set enrichment analysis based on the contribution of SGE to gene expression in the prefrontal cortex (PFC) in the experiment with inbred strains.** The analysis was carried out in B6 and D2 focal mice separately. Enrichment was calculated based on the contribution of SGE to gene expression variance (see Methods). The 5 most significant GO terms for each of strain are shown. * dopamine measured by HPLC in the midbrain is also strongly affected by SGE (Table 1).</caption><thead><tr><th>Focal mice</th><th>GO term ID</th><th>GO term</th><th># annotated genes</th><th>P value</th><th>Q value</th></tr></thead><tbody><tr><td>B6</td><td>GO:0042053</td><td>regulation of dopamine metabolic process*</td><td>12</td><td>7.0 10<sup>-4</sup></td><td>1</td></tr><tr><td>B6</td><td>GO:0045216</td><td>cell-cell junction organization</td><td>127</td><td>2.8 10<sup>-3</sup></td><td>1</td></tr><tr><td>B6</td><td>GO:0072376</td><td>protein activation cascade</td><td>20</td><td>2.8 10<sup>-3</sup></td><td>1</td></tr><tr><td>B6</td><td>GO:1902017</td><td>regulation of cilium assembly</td><td>15</td><td>3.3 10<sup>-3</sup></td><td>1</td></tr><tr><td>B6</td><td>GO:0042987</td><td>amyloid precursor protein catabolic process</td><td>19</td><td>3.4 10<sup>-3</sup></td><td>1</td></tr><tr><td>D2</td><td>GO:0007229</td><td>integrin-mediated signaling pathway</td><td>57</td><td>2.5 10<sup>-6</sup></td><td>0.024</td></tr><tr><td>D2</td><td>GO:0043269</td><td>regulation of ion transport</td><td>431</td><td>1.7 10<sup>-4</sup></td><td>0.81</td></tr><tr><td>D2</td><td>GO:0030001</td><td>metal ion transport</td><td>543</td><td>1.4 10<sup>-3</sup></td><td>1</td></tr><tr><td>D2</td><td>GO:1900006</td><td>positive regulation of dendrite development</td><td>59</td><td>1.5 10<sup>-3</sup></td><td>1</td></tr><tr><td>D2</td><td>GO:0008542</td><td>visual learning</td><td>43</td><td>1.9 10<sup>-3</sup></td><td>1</td></tr></tbody></table>

doi:10.1371/journal.pgen.1006498.t002
## Reanalysis of a large outbred mice dataset

We next investigated SGE in a dataset from outbred (Heterogeneous Stock) mice[22–24]. This design better represents traditional mouse housing conditions (groups of four, five and six mice mostly rather than pairs; S4 Fig) and genetic variation in natural populations (high genetic diversity and genetically unique individuals rather than two inbred strains). The dataset comprises more than 100 organismal phenotypes measured in 2,448 mice (S5 Table), and gene expression in hippocampus for a subset of 457 mice. To accommodate the genetic design of this study, we fitted random effects models with variance components for DGE, SGE and their covariance. Our models are inspired from models published in the literature[2, 25] yet differ in some important ways, which we explain and justify in S1 Note. Because genome-wide genotypes were not available for a subset of mice (526 out of 2,448), we used pedigree information to estimate pairwise genetic covariance for the mice with no genotype information (see Methods). We found that the results obtained using both pedigree and genotype data for estimation of genetic similarity were in agreement with those obtained using genotype data only from the subset of mice that were in cages where all mice had been genotyped (S5 Fig). All models were fitted using LIMIX[26, 27].

Simulations showed that DGE and SGE were unbiased (S6 Fig and S7 Fig). Of the 117 organismal phenotypes available in this dataset, 43 were significantly affected by SGE ($P < 0.05$, FDR $< 5.7\%$; Table 3, S5 Table). SGE explained up to 29% of phenotypic variation and an average of 8.9% across the 43 significantly affected traits. Importantly, the estimated contribution of SGE was greater than that of DGE for 8 of the 43 traits.

Among the organismal phenotypes most significantly and strongly affected by SGE in this dataset were measures of lymphocyte activation (in particular size and number of CD4+ T cells and B cells, collected by fluorescence-activated cell sorting and full blood count, Table 3). In contrast, measures related to other leucocytes (neutrophils, basophils and monocytes) and natural killer T cells. Altogether, these results indicate that genotypes of cage mates influence humoral immunity, although this is unlikely that this represents spread of a disease as the mice were kept in a clean mouse facility.

In addition, rate of wound healing, the measure most significantly and strongly associated with SGE in the experiment with two inbred strains, was also significantly affected by SGE in the outbred dataset ($P = 1.6 \times 10^{-3}$, $Q = 3.7 \times 10^{-3}$).

Three measures of body weight (collected on weeks 6 and 7) also figured among the traits significantly affected by SGE. Other measures of body weight, collected at weeks 9 and 10, showed no SGE (S5 Table). The two sets of measures likely reflect a different phenotype as the normal physiology of the mice was disrupted between weeks 7 and 9 by aggressive phenotyping (tests of conditioned anxiety involving foot shocks, airway sensitization by allergen, intraperitoneal injection of glucose). Another potential but unlikely explanation for the higher contribution of SGE on earlier body weight measures is a sharp decrease in social effects between weeks 7 and 9.

Finally, a subset of measures of anxiety, blood biochemistry, and lung function were also affected by the genotypes of cage mates.

In contrast, there was no statistical evidence for SGE on gene expression levels in the hippocampus (smallest nominal P value $4.1 \times 10^{-5}$ corresponding to a Q value of 64%, S6 Table). This is most likely the result of a much smaller sample size (457 or 5% of mice) for expression traits, and is consistent with published power analyses[28].

We next explored whether studies focused on DGE can safely ignore SGE, focusing on the estimation of the collective effect of additive DGE on phenotypic variation (narrow-sense heritability). In the outbred dataset, cage mates are genetically more similar to each other than
<table><caption>**Table 3. Organismal phenotypes significantly affected by SGE in the outbred mice dataset (P < 0.05).** Q value: see Methods. STE: standard error. Bold text highlights eight phenotypes for which the contribution of SGE is greater than that of DGE (contribution of SGE-STE > contribution of DGE + STE).</caption><thead><tr><th>Measure</th><th>Measure of</th><th>SGE P value</th><th>SGE Q value</th><th>SGE variance (%)</th><th>DGE variance (%)</th></tr></thead><tbody><tr><td>CD4 fluorescence intensity</td><td>Lymphocyte function</td><td>&lt;1.00 10<sup>-16</sup></td><td>&lt;1.00 10<sup>-16</sup></td><td>29 +/- 5</td><td>21 +/- 2</td></tr><tr><td>Glycemia before glucose injection</td><td>Basal glycemia</td><td>1.30 10<sup>-10</sup></td><td>3.20 10<sup>-09</sup></td><td>11 +/- 4</td><td>11 +/- 2</td></tr><tr><td>Size of CD4+ cells (forward scatter)</td><td>Lymphocyte function</td><td>2.50 10<sup>-08</sup></td><td>4.40 10<sup>-07</sup></td><td>27 +/- 6</td><td>5 +/- 1</td></tr><tr><td>Size of B cells (forward scatter)</td><td>Lymphocyte function</td><td>4.00 10<sup>-08</sup></td><td>5.20 10<sup>-07</sup></td><td>25 +/- 6</td><td>5 +/- 1</td></tr><tr><td>Size of CD8+ cells (forward scatter)</td><td>Lymphocyte function</td><td>8.00 10<sup>-08</sup></td><td>8.20 10<sup>-07</sup></td><td>27 +/- 6</td><td>4 +/- 1</td></tr><tr><td>Chloride</td><td>Blood biochemistry</td><td>7.70 10<sup>-07</sup></td><td>6.60 10<sup>-06</sup></td><td>17 +/- 5</td><td>4 +/- 2</td></tr><tr><td>Calcium</td><td>Blood biochemistry</td><td>9.90 10<sup>-07</sup></td><td>7.30 10<sup>-06</sup></td><td>18 +/- 5</td><td>5 +/- 2</td></tr><tr><td>Absolute number of lymphocytes</td><td>Lymphocyte function</td><td>3.60 10<sup>-06</sup></td><td>2.30 10<sup>-05</sup></td><td>7 +/- 4</td><td>9 +/- 3</td></tr><tr><td>Weight at 6 weeks of age</td><td>Body weight</td><td>5.30 10<sup>-06</sup></td><td>3.00 10<sup>-05</sup></td><td>7 +/- 4</td><td>11 +/- 2</td></tr><tr><td>Sodium</td><td>Blood biochemistry</td><td>5.90 10<sup>-06</sup></td><td>3.00 10<sup>-05</sup></td><td>17 +/- 5</td><td>2 +/- 2</td></tr><tr><td>Alanine transaminase</td><td>Liver function</td><td>1.00 10<sup>-05</sup></td><td>4.70 10<sup>-05</sup></td><td>11 +/- 4</td><td>2 +/- 2</td></tr><tr><td>White blood cells</td><td>Lymphocyte function</td><td>3.50 10<sup>-05</sup></td><td>0.00014</td><td>6 +/- 4</td><td>7 +/- 3</td></tr><tr><td>Mean startle (following bang only) before conditioning with foot shock</td><td>Unconditioned anxiety</td><td>0.00011</td><td>4.00 10<sup>-04</sup></td><td>15 +/- 4</td><td>24 +/- 4</td></tr><tr><td>Weight before food hyponeophagia test (7 weeks)</td><td>Body weight</td><td>0.00015</td><td>0.00052</td><td>12 +/- 4</td><td>6 +/- 2</td></tr><tr><td>Glycemia 15 minutes after glucose injection</td><td>Glucose tolerance</td><td>0.00035</td><td>0.0011</td><td>9 +/- 4</td><td>9 +/- 3</td></tr><tr><td>High density lipoprotein</td><td>Blood biochemistry</td><td>0.00046</td><td>0.0014</td><td>8 +/- 4</td><td>35 +/- 3</td></tr><tr><td>Proportion of CD3+ T-cells that are CD4+</td><td>Lymphocyte function</td><td>0.00049</td><td>0.0014</td><td>7 +/- 5</td><td>9 +/- 3</td></tr><tr><td>Red blood cell distribution width</td><td>Red blood cell function</td><td>0.00053</td><td>0.0014</td><td>6 +/- 4</td><td>28 +/- 4</td></tr><tr><td>Weight after food hyponeophagia test</td><td>Body weight</td><td>0.00093</td><td>0.0024</td><td>10 +/- 4</td><td>7 +/- 2</td></tr><tr><td>Percentage of B220+ cells</td><td>Lymphocyte function</td><td>0.0011</td><td>0.0027</td><td>10 +/- 5</td><td>32 +/- 4</td></tr><tr><td>Wound (ear hole) area</td><td>Wound healing</td><td>0.0016</td><td>0.0037</td><td>6 +/- 3</td><td>32 +/- 3</td></tr><tr><td>Expiratory time (after metacholine)</td><td>Lung function</td><td>0.0018</td><td>0.0039</td><td>3 +/- 3</td><td>11 +/- 3</td></tr><tr><td>Enhanced pause (before metacholine)</td><td>Lung function</td><td>0.0019</td><td>0.0041</td><td>10 +/- 4</td><td>9 +/- 3</td></tr><tr><td>Distance traveled in the center of the EPM</td><td>Unconditioned anxiety</td><td>0.0021</td><td>0.0043</td><td>4 +/- 3</td><td>27 +/- 4</td></tr><tr><td>Mean startle (following tone + bang) before conditioning with foot shock</td><td>Unconditioned anxiety</td><td>0.0032</td><td>0.0063</td><td>9 +/- 4</td><td>25 +/- 4</td></tr><tr><td>Percentage of CD8+ cells</td><td>Lymphocyte function</td><td>0.0039</td><td>0.0073</td><td>4 +/- 3</td><td>48 +/- 4</td></tr><tr><td>Distance traveled in the closed arms of the EPM</td><td>Locomotor activity</td><td>0.0042</td><td>0.0077</td><td>2 +/- 3</td><td>18 +/- 3</td></tr><tr><td>Distance traveled in the open arms of the EPM</td><td>Unconditioned anxiety</td><td>0.0057</td><td>0.01</td><td>2 +/- 3</td><td>19 +/- 3</td></tr><tr><td>Number of entries into the open arms of the EPM</td><td>Unconditioned anxiety</td><td>0.0063</td><td>0.011</td><td>2 +/- 3</td><td>19 +/- 3</td></tr><tr><td>Albumin</td><td>Blood biochemistry</td><td>0.0066</td><td>0.011</td><td>5 +/- 4</td><td>4 +/- 2</td></tr><tr><td>Number of K167+ cells</td><td>Adult neurogenesis</td><td>0.01</td><td>0.017</td><td>4 +/- 6</td><td>38 +/- 7</td></tr><tr><td>Measured mean cell hemoglobin concentration</td><td>Red blood cell function</td><td>0.012</td><td>0.019</td><td>7 +/- 4</td><td>15 +/- 3</td></tr><tr><td>Beam breaks from ambulation in home cage during last 5 minutes</td><td>Locomotor activity</td><td>0.013</td><td>0.02</td><td>2 +/- 3</td><td>11 +/- 3</td></tr><tr><td>Faecal steroids in males and females</td><td>Stress</td><td>0.015</td><td>0.022</td><td>6 +/- 5</td><td>8 +/- 4</td></tr><tr><td>Time spent in the open arms of the EPM</td><td>Unconditioned anxiety</td><td>0.015</td><td>0.022</td><td>1 +/- 3</td><td>17 +/- 3</td></tr><tr><td>Aspartate transaminase</td><td>Blood biochemistry</td><td>0.023</td><td>0.032</td><td>4 +/- 4</td><td>2 +/- 2</td></tr><tr><td>Alkaline phosphatase</td><td>Blood biochemistry</td><td>0.026</td><td>0.036</td><td>1 +/- 3</td><td>44 +/- 3</td></tr><tr><td>Creatinine</td><td>Blood biochemistry</td><td>0.027</td><td>0.036</td><td>5 +/- 5</td><td>2 +/- 3</td></tr><tr><td>Glucose</td><td>Basal glycemia</td><td>0.032</td><td>0.041</td><td>11 +/- 4</td><td>9 +/- 3</td></tr></tbody></table>

(Continued)
<table><caption>**Table 3.** (Continued)</caption><thead><tr><th>Measure</th><th>Measure of</th><th>SGE P value</th><th>SGE Q value</th><th>SGE variance (%)</th><th>DGE variance (%)</th></tr></thead><tbody><tr><td>Time freezing during cue</td><td>Conditioned anxiety</td><td>0.034</td><td>0.042</td><td>7 +/- 4</td><td>15 +/- 4</td></tr><tr><td>Percentage of CD3+ cells</td><td>Lymphocyte function</td><td>0.034</td><td>0.042</td><td>5 +/- 4</td><td>31 +/- 4</td></tr><tr><td>Respiratory rate (after metacholine)</td><td>Lung function</td><td>0.04</td><td>0.048</td><td>1 +/- 3</td><td>17 +/- 4</td></tr><tr><td>Expiratory time (before metacholine)</td><td>Lung function</td><td>0.049</td><td>0.057</td><td>2 +/- 3</td><td>15 +/- 3</td></tr></tbody></table>

doi:10.1371/journal.pgen.1006498.t003

average (S8 Fig). As a result, DGE and SGE are correlated (an example of gene-environment correlation). Thus, we hypothesized that failing to account for SGE would bias estimates of DGE. We also hypothesized that fitting cage effects might be sufficient to eliminate this bias. To investigate both hypotheses, we compared DGE estimates obtained using a linear mixed model for DGE that does not account for cage effects nor SGE, a model that accounts for DGE and cage effects, and one that accounts for DGE, SGE, corresponding social environmental effects, and cage effects (“full model”). Models that did not account for SGE led to substantially larger DGE estimates, and estimates from the model with DGE and cage effects were intermediate between those from the model with DGE only and those from the full model. In simulated traits based on the real genotypes and generated from DGE, SGE and cage effects (see Methods), we found that models that did not account for SGE yielded inflated DGE estimates, whereas joint modeling of DGE, SGE and cage effects resulted in unbiased estimates (Fig 2C and S6 Fig). Importantly, fitting cage effects but no SGE did not eliminate the bias. The simulation results strongly suggest that, in the real data, the estimates obtained from the full model are most accurate, and that models that ignore SGE overestimate heritability. This problem is particularly acute when direct and social random genetic effects are positively correlated (i.e. $\sigma_{A_{DS}} > 0$, see Methods; Fig 2A and 2B). The problem we highlight here is general and likely to affect other studies in which social partners are related, including twin and family studies used to estimate heritability in humans (see Discussion).

## Discussion

Using two complementary genetic designs—one using two mouse inbred strains and one using outbred mice—we estimated the contribution of social genetic effects to a variety of organisational phenotypes and gene expression traits. The experiment with two inbred strains was designed to investigate SGE and focused on behaviours (anxiety and helplessness), as there is strong evidence that behaviours are socially affected[29–33]. To test whether SGE can be detected in outbred populations and survey a broader range of phenotypes, we re-analysed a large dataset from outbred mice and quantified the contribution of SGE to more than 100 phenotypes. The design of our study raises important questions: are positive results (i.e. evidence of SGE) in the experiment with two inbred strains expected to replicate in the outbred dataset? Are some phenotypes expected to be affected by SGE and some not? We now discuss these points.

Some phenotypes were measured in both experiments with similar protocols. Wound area, the measure of wound healing, was collected in both experiments using the exact same protocol. It is significantly affected by SGE in both experiments. The protocols for measuring body and adrenal gland weight are fairly simple thus reducing technical variation between experiments, and body weight was measured at about the same time point (around 50 days of age) in both experiments. There was strong evidence for SGE on body weight in the outbred dataset but no evidence in the experiment with two inbred strains. No SGE on adrenal gland weight
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Ignoring SGE leads to biased estimates of DGE (heritability) in the outbred mice dataset. (A,B) DGE estimates for all organismal phenotypes. The colour of the dots indicates the sign of the covariance $\sigma_{ADs}$ between DGE and SGE (blue: negative; red: positive; see Methods); the intensity of the colour indicates the magnitude of SGE. (A) Comparison of DGE estimates from model with DGE only (x-axis) and full model with DGE, SGE, (y-axis). (B) Comparison of DGE estimates from model with DGE and cage effects (x-axis) and full model (y-axis). (C) Difference between estimated and simulated DGE. The contribution of SGE to simulated phenotypes varied and is shown on the x-axis. See Methods for values given to other parameters. Simulations were analyzed with three different models: DGE only (red), DGE + cage effects (green), DGE + SGE + cage effects + social environmental effects (blue).</figcaption></figure>

doi:10.1371/journal.pgen.1006498.g002

were detected in either dataset. Finally, a partially overlapping set of measures of unconditioned anxiety was significantly affected by SGE in both experiments. While reviewing results
from the two experiments in parallel is informative, positive results in one experiment are not strictly expected to replicate in the other. Indeed, although the variants that give rise to SGE in the experiment with two inbred strains also segregate in the outbred population (the two strains used in our experiment were among the eight founders of the outbred population), they have recombined with many additional variants from the six other founders. Moreover, the housing conditions were very different in the experiment with inbred strains and outbred experiment (group size of 2 vs. 2 to 7 respectively, and unfamiliar mice vs. familiar mice housed together). Therefore, one should not expect the overall contribution of SGE be the same in the two experiments. Rather, combining the two experiments provides a first hint at the generalizability of our results. Published studies of SGE provide additional information on this matter, and suggest that SGE may contribute to variation in body weight across species[2].

That social genetic effects contribute to variation in anxiety probably does not come as a surprise but their contribution to wound healing maybe more so. This result is however supported by significant p-values in both experiments of our study and large effect sizes (18 and 6%). When interpreting this result, it is important to bear in mind that social effects on wound healing can (and will, necessarily) be mediated by traits of cage mates that are different from wound healing. For example, social effects on wound healing could be mediated by social grooming, which could either mechanically disrupt the healing process or chemically enhance it[34]. Any traits of cage mates that may induce a systemic stress response in the focal animal could also mediate social effects on wound healing[35, 36]. Thus, social effects on wound healing are not unlikely, and, similarly, social effects may affect any phenotype (e.g. by through the induction of a systemic stress response,).

Because any phenotype may *a priori* be affected by social effects and the mechanisms at play are rarely known, SGE offer an attractive alternative to investigate social effects. First, as we have shown, they can be used to quantify social effects, effectively providing a lower bound estimate of social effects (as only the genetic component is captured). Second, SGE can be used to test whether a particular trait of social partners has an effect on a phenotype of interest. Establishing a causal relationship between two phenotypes is always difficult because of the risk of reverse causation and independent action of hidden confounders on both traits; SGE provide an anchor to test causality.

Independent of their relevance for studying social effects, we show that ignoring SGE can lead to biased estimates of heritability (i.e. the collective effect of DGE). In our study (outbred dataset), DGE and SGE are correlated by design (mice that share a cage are more genetically similar than average), and we show that this correlation leads to biased estimates of heritability if unaccounted for. Fitting cage effects, which has the primary goal of accounting for environmental effects shared by cage mates (e.g. noise levels), does not eliminate the bias. Our results are of interest to the broad genetics community as DGE and SGE are correlated in most if not all experimental designs traditionally used to estimate heritability in humans and model organisms, and SGE may thus have caused widespread bias. For example, in twin designs, MZ twins not only share 100% of their genotypes but they also share 100% of the genotypes of their sibling; DZ twins in comparison share both 50% of their genotypes and 50% of the genotypes of their sibling. Thus, SGE can contribute to increased concordance between MZ twins compared to DZ twins. If SGE are not modelled, heritability may be overestimated and appear “missing” when compared to genome-wide association results obtained from unrelated individuals[37]. Note that when the covariance $\sigma_{ADS}$ between direct and social random genetic effects is negative (competition effects), ignoring SGE may lead to underestimating heritability. SGE in humans were considered once before (“sibling effects” [38]) but were never, to the best of our knowledge, modelled in heritability studies. Because we found that fitting cage effects
was not sufficient to eliminate the bias due to SGE, we suspect that accounting for a “common environment” shared by family members, as is commonly done in human studies[39–41], will not eliminate SGE-induced bias.

It is not the first time that unaccounted for gene-environment correlations are put forward as potential causes of bias (e.g. Conley *et al.* investigated the correlation between genetics and urban setting [42]). However, the impact of the correlation between DGE and SGE is likely to be particularly severe as we have shown that SGE affect a wide range of phenotypes and DGE and SGE are correlated in most experimental designs used to estimate heritability.

Our study sheds light on an important component of the genetic architecture of complex traits, one that lies outside the individual, in social partners. Social genetic effects have already been shown to play an important role in artificial selection of livestock[43] and have important evolutionary consequences[44, 45]. Our results provide evidence that SGE are also an important component of health and disease.

# Methods

## Experiment with inbred strains

**Animals and experimental design.** Forty-five C57BL/6J (B6) and 45 DBA/2J (D2) about 3-weeks-old (±3 days) were shipped from The Jackson Laboratory. Mice were weaned on the day they were shipped and littermates were shipped together so as to minimise novel social interactions prior to experiments. Three D2 mice died during shipment and to maintain balance one B6 mouse was excluded. On arrival at UTHSC, mice were transferred to standard (75 square inches) mouse cages and allowed to recover overnight with their original littermates. Ears were punched for identification, and mice were transferred to new cages. Mice were co-housed in pairs to make up sixteen B6/D2 (BD) cages, 14 B6/B6 (BB) cages, and 13 D2/D2 (DD) cages assembled using mice from 15 B6 and 15 D2 litters (3 mice per litter originally). Littermates were never paired in BB and DD cages so cage mates were always unfamiliar with each other at the beginning of the experiment. Littermates were distributed across the three groups (BB, BD, and DD), and the cages ordered as successive BB, BD, and DD triplets on the shelves. Mice were housed in standard mouse cages with corncob bedding and a cardboard shelter on a fixed 12:12 hour light: dark cycle with *ad libitum* access to food and water. Phenotypes were collected after six weeks of co-housing.

**Phenotyping.** All protocols were approved by the UTHSC Institutional Animal Care and Use Committee. Mice were phenotyped between 7am and 2pm, in a random order different for each test. The order in which tests were performed, as well as the measures collected in each test, are summarised in S1 Table.

**Elevated plus maze:** We used an elevated plus maze (Columbus Instruments) with a surface ~30 cm above the floor that contained four 51-cm long, 11.5 cm-wide arms arranged at right angles. Closed arms had opaque walls 30 cm high, extending the length of the arm. Each mouse was placed in the center of the maze facing an open arm and allowed 5 minutes to explore the maze. Each cage mate was tested in turn and precaution was taken so that cage mates did not interact just after the first cage mate had finished the test and before the second one was taken to the maze.

**Activity chamber:** We used six 17" X 17" Plexiglas activity chambers (Med Associates Inc., ENV-510) with photobeam around the edge to measure movement. The beams were interfaced to a computer to record movement. Each activity chamber was located into a closed section of a cupboard, which had white light on. Each mouse was placed in the centre of the chamber and given 30 minutes to explore the activity chamber. Cage mates were tested at the same time.
**Tail suspension test:** Adhesive tape was wrapped around each tail three quarters of the distance from the base, and used to attach a suspension hook. Animals were suspended for six minutes and video recorded from the side. Cage mates were tested at the same time.

**Forced swim test:** Mice were placed in six cylinders (height: 20 cm, diameter: 15 cm) with water at about 25°C for 6 minutes and their behaviour recorded from the side. Upon completion mice were dried. Cage mates were tested at the same time.

**Body weight:** Mice were weighed on the day after arrival at UTHSC and every week thereafter.

**Sacrifice, tissue collection, and gland weights:** Mice were euthanized by injection of 250 mg/kg Avertin (20 mg/ml; tribromoethanol) followed by cardiac puncture and exsanguination under deep anaesthesia. The blood was collected in lithium heparin tubes, and an aliquot was taken in another lithium heparin tube for serotonin quantification. The aliquot was immediately snap frozen in liquid nitrogen and then stored at -80C.

The brain was dissected and placed in an ice-cold coronal mouse brain matrix. A slab containing the prefrontal cortex (PFC) was removed, and the PFC taken, snap frozen and stored at -80°C. The part of the brain caudal to the amygdala was dissected, snap frozen and stored at -80°C. Then the midbrain was dissected by letting the brain soften on ice and removing the cerebellum. The midbrain was stored at -80°C and shipped on dry ice for analysis by high-pressure liquid chromatography.

Both preputial glands were dissected and weighed together, as were both gonads. Right and left adrenal glands were weighed separately.

**High-pressure liquid chromatography:** Brain samples were weighed frozen. They were then crushed in 1,5 ml of 10⁻³ M HCl containing sodium metabisulfite and ascorbic acid (antioxidants), and centrifuged at 22,000 g for 20min at 5°C. The supernatants were collected and filtered through a 10-kDa membrane (Nanosep, Pall) by centrifugation at 7,000 g. Then, a 20 µl aliquot of sample was analyzed for serotonin by fluorometric detection[46]. The amounts of catecholamines (dopamine and noradrenaline) and their metabolites (DOPAC, HVA, VMA) as well as the serotonin metabolite 5HIAA were measured by electrochemical detection on a serial array of coulometric flow-through graphite electrodes (CoulArray, ESA)[47]. The analysis, data reduction and peak identification were fully automated.

Ten µl of blood were used for serotonin quantification. Serotonin was extracted by adding 400 µl of 10⁻³ M HCl, followed by 2 centrifuge runs including one on Nanosep 10 kD. Quantification was carried out as for brain serotonin.

**Size of ear hole:** At sacrifice the ear that was clipped on arrival was cut and placed in formalin. It was then mounted on a glass slide, scanned, and the resulting image analysed with ImageJ. The scale was set once using the size of the slide. The hole was delimited using the freeform drawing tool.

**Pre-processing of organismal phenotypes.** All collected phenotypes were quantitative. Individual phenotypes were normalized using the parametric Box-Cox power transformation [48], using R's MASS package[49] and a model where strain of the animal, strain of the cage-mate, and an interaction between the two were included as fixed effects to inform the Box-Cox fit. Note that this approach does not regress out these genetic effects.

On the Box-Cox transformed data for each phenotype, we tested for association with experimental covariates and regressed out significantly associated covariates that were unrelated to the measure of interest. The covariates that were accounted for are reported in S1 Table. Significance of SGE calculated by including covariates in the models rather than regressing them out are reported in S3 Table.

**RNA extraction and sequencing.** Total RNA was extracted using a Qiagen Qiacube and RNeasy mini kit according to the manufacturers protocols and was loaded on a 96-well plate,
with samples randomly allocated to wells, for shipment to the sequencing centre. The quantity and quality of the samples was assessed using Agilent high sensitivity R6K screen tape. One sample (B42) with high degradation was excluded. Two samples (B39, D30) had very low RNA integrity numbers (RIN, 3.7 and 4.2 respectively), and one sample (B9) had low 28S/18S peak height ratio but they were included in library preparation and sequencing. All libraries were prepared according to the dUTP strand-specific protocol, multiplexed in one 85-plex and sequenced in 5 HiSeq 2500 rapid runs to obtain paired-end, 51 bp reads.

**Alignment of sequencing reads and quantification of gene expression.** We used SNPs and indels identified in D2 (Sanger Mouse Project[50], files mgp.v3.snps.rsIDdbSNPv137.vcf and mgp.v3.indels.rsIDdbSNPv137.vcf, both on GRCm38) to construct a D2 genome and the corresponding gene annotation file using the software SeqNature[51] (version 1.2). Reads were aligned using Tophat[52] (Tophat2 version 2.0.11, Bowtie version 2.2.2.0, and Samtools version 0.1.19.0), considering either the B6 genome (GRCm38) or the D2 genome. High quality alignments were selected using the following criteria: both reads of the pair had to uniquely map to the same chromosome and within 2.3 Mb of each other (length of the longest gene). Expression levels were quantified at the gene level using HTSeq[53] (version 0.6.1), considering high quality alignments only and UCSC annotations (GRCm38/mm10).

**Pre-processing and quality control.** Raw expression counts were adjusted for library size using the R package DESeq2[54] (version 1.2.10). Out of 23,420 genes, we considered 13,271 genes with at least 10 (library size normalized) reads in at least 40% of the samples. Gene expression levels were normalized using the Box-Cox transformation as done for the organis- mal phenotypes.

We identified and excluded genes whose expression was affected by the order of sacrifice: at the end of our experiment, we sacrificed each cage mate in turn with an interval of about 15 minutes between the first and the second cage mate. The cage mate that was sacrificed second was in the room and we found evidence that the expression of a small number of genes was activated during this interval (i.e. their expression was correlated with order of sacrifice within cage mate pair). Some of the affected genes were well-known immediate early genes. In our experiment, this effect was confounded with potential social genetic effects as B6 mice in BD cages were always sacrificed first and D2 mice second, while in BB cages some of the B6 mice were necessarily sacrificed first and others second, and the same applies for D2 mice in DD cages. To avoid possible false positive SGE, we took a conservative approach and excluded all genes for which the order of sacrifice within pair, in B6 mice from BB cages and D2 mice from DD cages, was associated with gene expression levels (P < 0.05 P, > 1% variance explained; 375 genes in total).

We used a principal component analysis (PCA) to identify outlier samples. The first and second principal components explained 22% and 9% of variation in gene expression, and were used to define outliers (S9 Fig). A total of 7 samples were excluded (B39, D13, D30, D9, B30, D16, D39), two of which also had very low RIN. The final gene expression dataset consisted of 11 pairs of B6 cage mates, 3 B6 mice whose B6 cage mates were not included in the final sam- ple, 15 pairs of B6/D2 cage mates, 1 B6 mouse whose D2 cage mate was not included in the final sample, 9 pairs of D2 cage mates, and 4 D2 mice whose cage mates were not included in the final sample. This amounts to 25, 16, 15, and 22 mice in the BB, BD, DB, and DD groups respectively.

In addition to removing outlier samples, we used two approaches to adjust for observed and hidden covariates affecting the expression profiles. First, for significance tests of SGE on gene expression levels, we used PEER to account for hidden sources of variation[55]; we ran PEER while including strain, strain of the cage mate, and interaction between them as observed covariates so as not to remove the genetic signal (see below), and used the expression
residuals calculated from the first 5 PEER factors for statistical tests. The p-values obtained with PEER residuals were calibrated whereas other normalization strategies led to severely deflated test statistics, supporting that the PEER factors accounted for unwanted variation such as hidden batch effects.

Second, for variance partitioning and model selection, we regressed out observed technical covariates only as removal of hidden factors can reduce the residual noise. We considered several technical covariates: initial RNA concentration (as an indicator of the size of tissue dissected, since the same volume of buffer was added to all), RIN, ratio between heights of 28S and 18S peaks, proportion of first reads mapping—as they should—to the strand opposite to the feature (i.e. strand specificity of the protocol, calculated by running htseq-count a second time with the--stranded = yes option), insert size, yield, proportion of reads mapped to the (B6 or D2) genome, proportion of mapped reads that are of high quality (as defined above), proportion of high quality reads that are exonic, number of alignments with 0, 1, or 2 errors (as calculated by SAMStat[56] version 1.09), and the following measures calculated from the output of FastQC[57] (version 0.10.1): average per base sequence quality, average per base sequence (G nucleotide used) content, average per base GC content, average per sequence GC content, average per base N content, average per sequence quality scores, average sequence duplication levels, presence or absence of overrepresented sequences. We calculated the correlation of these covariates with the first principal components, and found that *proportion of selected read pairs that are exonic* was significantly correlated with the first principal component (Pearson correlation, R² = 0.66, P = 8.10⁻¹²). Therefore we regressed the effect of this covariate out of the expression of each gene, and used the residuals for variance partitioning and model selection.

**Statistical modelling.** Conceptually we distinguish between the **focal individual** (f), whose phenotype is the **phenotype of interest** and the dependent variable in the model, and the focal individual's **social partners**—here its cage mate (cm)—which may influence the phenotype of the focal individual through a number of unknown **mediating phenotypes**. No prior knowledge or measurement of these mediating phenotypes is required as they are captured through their genetic and environmental components (social genetic effects, SGE, and social environmental effects, SEE, respectively).

Genetic effects considered in this experiment include:

* direct genetic effects (DGE, effects of the strain of the focal animal on its own phenotype—the phenotype of interest). This component is only relevant when both B6 and D2 focal mice are analysed together.

* social genetic effects (SGE, effects of the strain of the cage mate on the focal animal's phenotype; SGE act through mediating phenotypes of the cage mate).

* an interaction between DGE and SGE (whereby the effect of the strain of the cage mate depends on the strain of the focal individual). This component is only relevant when both B6 and D2 focal mice are analysed together.

Genetic effects were modelled as fixed effects.

"Environmental" (i.e. unexplained or residual) variance may arise from:

* direct environmental effects (DEE, effects of the "environment" of the focal animal on its own phenotype)

* social environmental effects (SEE, effects of the environment of the cage mate on the focal animal's phenotype; SEE are the environmental component of the mediating phenotypes, so whenever social effects are considered both SGE and SEE need to be accounted for)
* cage effects (which affect the phenotype of interest of all mice sharing a cage in the same way)

Environmental effects were modelled as random effects.

The model of phenotypic variation that includes all of the genetic and environmental effects listed above (note: a simpler model was effectively used for analysis, as explained next) is the following:

$$y_f = A_{D,f} + A_{S,cm} + A_{D,f} : A_{S,cm} + e_{D,f} + Z_f e_S + W_f c$$

where $y_f$ is the phenotype of interest of the focal animal $f$, $A_{D,f}$ and $A_{S,cm}$ respectively refer to additive DGE of the focal individual and additive SGE from the cage mate $cm$; $A_{D,f}$: $A_{S,cm}$ refers to the interaction between $A_{D,f}$ and $A_{S,cm}$; $Z_f$ is a row of the matrix $Z$ that indicates cage mates (importantly $Z_{i,i} = 0$); $e_D$ is the vector of random DEE and $e_S$ the vector of random SEE. $W_f$ is a row of the matrix $W$ that indicates cage assignment and $c$ the vector of cage effects.

The joint distribution of the environmental random effects is defined as:

$$\begin{bmatrix} e_D \\ e_S \\ c \end{bmatrix} \sim \text{MVN}(0, \begin{bmatrix} \sigma_{E_D}^2 & I & \sigma_{E_{DS}}^2 & I & 0 \\ \sigma_{E_{DS}}^2 & I & \sigma_{E_S}^2 & I & 0 \\ 0 & 0 & \sigma_{C}^2 & I & 0 \end{bmatrix})$$

The covariance $\sigma_{E_{DS}}$ between the random direct and social environmental effects exists because an environmental factor may influence both the phenotype of interest and the mediating traits, in which case it will influence the phenotype of interest both in a direct ($e_D$) and an indirect ($e_S$) manner. Depending on 1) the direction of the environmental effect on the phenotype of interest, 2) the directions of the effects on the mediating traits, and 3) the signs of the correlations between mediating traits and phenotype of interest, $\sigma_{E_{DS}}$ may be positive or negative.

The environmental covariance is:

$$\begin{align*} \operatorname{cov}(e_i, e_j) &= \sigma_{E_D}^2 + \sigma_{E_S}^2 + \sigma_C^2 & \text{if } i = j \\ \operatorname{cov}(e_i, e_j) &= 2 \sigma_{E_{DS}} + \sigma_C^2 & \text{if } i \neq j \text{ and } i \text{ and } j \text{ share a cage} \\ \operatorname{cov}(e_i, e_j) &= 0 & \text{if } i \text{ and } j \text{ are in different cages} \end{align*}$$

Note that $\operatorname{cov}(e_i, e_j)$ may be negative as a result of $\sigma_{E_{DS}}$ being negative.

When group size is constant, as is the case in this dataset (group size = 2), $e_D$, $e_S$, and $c$ are not identifiable. Thus we instead used a single random environmental term, whose correlation structure is described below, as has been done in other studies of social genetic effects[25]:

$$\begin{align*} \operatorname{cov}(e_i, e_j) &= \sigma_E^2 & \text{if } i = j \\ \operatorname{cov}(e_i, e_j) &= \rho \sigma_E^2 & \text{if } i \neq j \text{ and } i \text{ and } j \text{ share a cage} \\ \operatorname{cov}(e_i, e_j) &= 0 & \text{if } i \text{ and } j \text{ are in different cages} \end{align*}$$

From the expressions of $\operatorname{cov}(e_i, e_j)$ above we see that $\sigma_E^2 = \sigma_{E_D}^2 + \sigma_{E_S}^2 + \sigma_C^2$ and $\rho = (2 \sigma_{E_{DS}} + \sigma_C^2)/\sigma_E^2$. As the determinant of the environmental covariance matrix is necessarily positive, $-1 \le \rho \le 1$.

Note that under the assumption $0 \le \rho \le 1$ the environmental component could be rewritten as a sum of cage effects and iid residuals. This assumption is often made when group size is
large (e.g. [2]) because

$$ \rho = (2 \sigma_{E_{DS}} + (n-2)\sigma_{E_S}^2 + \sigma_C^2) / \sigma_E^2 $$

where n is the group size. $(n-2)\sigma_{E_S}^2$ is likely to be positive and large for large groups, so $\rho$ is likely be positive for large groups. However, in this experiment group size is 2 so $(n-2)\sigma_{E_S}^2 = 0$ and $\rho$ may well be negative. Thus we cannot rewrite the residual term as a sum of cage effects and iid residuals.

To investigate interactions between DGE and SGE (see Methods sections Model selection and Variance partitioning below), we considered all mice (B6 and D2) as focal individuals. For all other analyses we analysed B6 focal mice and D2 focal mice separately. We did so because there was evidence that interactions exist between DGE and SGE (results reported in the columns under “All mice” in S2 Table and shown in S2 Fig) and analysing B6 and D2 focal mice separately in that case simplified results interpretation.

When analysing B6 and D2 focal mice separately, we used the following simpler model:

$$ \begin{aligned} y_f &= A_{S,cm} + e_f \\ &= \beta_S X_{cm} + e_f \end{aligned} \qquad (1) $$

where $\beta_S$ is the regression coefficient of SGE and $X_{cm}$ the strain of the cage mate cm, with strain encoded as 0 for B6 and 1 for D2.

**Model selection.** To investigate whether interactions between DGE and SGE might contribute to phenotypic variation, we sought to identify the combination of genetic and environmental terms from {DGE, SGE, interaction between DGE and SGE, DEE, SEE, covariance between DEE and SEE, and cage effects} that best modelled phenotypic variation. All mice (B6 and D2) were included as focal individuals for this analysis. Constraints were imposed on the environmental correlation parameter $\rho$ when SEE and/or cage effects were not included in the model. Thus, the models were fitted in R using either the functions lm (stats package), lme (nlme package), or gls (nlme package, corCompSymm correlation structure for the residuals), as appropriate. The best model was selected using the Akaike information criterion (AIC) calculated based on the maximum-likelihood (ML) fit using the AIC function in the R stats package. The results of this analysis are reported in S2 Fig and the columns under “All mice” in S2 Table).

**Variance partitioning.** For S2 Fig and the columns under “All mice” in S2 Table, both B6 and D2 mice were considered as focal individuals and we used the following model:

$$ \begin{aligned} y_f &= A_{D,f} + A_{S,cm} + A_{D,f} : A_{S,cm} + e_f \\ &= \beta_D X_f + \beta_S X_{cm} + \beta_{D,S} X_f X_{cm} + e_f \end{aligned} \qquad (2) $$

where $\beta_D$ is the regression coefficient of DGE, $X_f$ the strain of f, and $\beta_{D,S}$ the regression coefficient of the interaction between DGE and SGE. $X_f X_{cm}$ is the product of the numerically encoded strains of f and cm, and therefore was 0 except for D2 focal mice with D2 cage mates.

All other results were obtained by analysing B6 and D2 focal mice separately using Model (1).

All models were fitted using restricted maximum likelihood (REML). The proportion of phenotypic variance explained by each fixed effect in the model (Model (1) or Model (2)) was calculated as the ratio of the sample variance of that term (i.e. the expected value of the variance of that term, which is var($\beta_S X$) for SGE) and the sample phenotypic variance. We obtained the sample phenotypic variance by drawing from the fitted model and calculating the average variance of the draws.

Standard errors were obtained using Gaussian propagation of uncertainty.
**Statistical significance tests.** Statistical significance analyses for SGE were carried out in B6 and D2 focal mice separately (S1 Fig), based on Model (1). Significance of the SGE covariate was assessed using a likelihood ratio test. Q values[58] implemented in the R package qvalue [59] were used to adjust for multiple testing across traits and focal strain (43 traits x 2 focal strains, 86 tests in total).

**Gene set enrichment analysis.** We considered B6 and D2 focal mice separately for this analysis, using the R package topGO[60] and the annotation package org.Mm.eg.db to perform gene set enrichment analyses. We used the residuals from PEER and considered the variance explained by SGE in Model 1 as score for each gene. Enrichment P values were obtained with the Kolmogorov-Smirnov test and the *elim* algorithm that decorrelates the GO graph structure[61]. Q values[58] were used for multiple testing correction, and obtained using the R package package[59].

## Outbred mice dataset

**Description of the dataset.** We used the raw data from Valdar *et al.*[22] (genotypes and organismal phenotypes), Huang *et al.*[23] (organismal phenotype: cellular proliferation in the subgranular zone of the dentate gyrus, a measure of adult neurogenesis) and Huang *et al.*[24] (gene expression in the hippocampus), which were collected on Heterogeneous Stock (HS) mice. HS mice are descended from eight inbred progenitor strains through more than 50 generations of circular breeding with 24 families[62]. As a result, each HS mouse is genetically unique. The mice included in this dataset are from a colony established at Oxford for a few generations, and for which pedigree data are available (3,149 mice). Cage information was available for a subset of 2,448 mice, filling 549 cages in groups of 2 to 7 mice (S4 Fig). Genotypes at 13,459 single nucleotide polymorphisms were available for a subset of 1,940 of those 2,448 mice, and a pair-wise genetic similarity matrix shows that they were related at different levels (S8 Fig), including (but not only) many full siblings. Some phenotype data were available for all 2,448 mice, but the number of mice phenotyped for each organismal phenotype is reported in S5 Table.

**Inference of relationship coefficients using the single-step or H matrix.** Our analysis of this dataset is based on linear mixed models, which use a genetic covariance matrix made of pair-wise relationship coefficients. These coefficients can be estimated from genome-wide genotypes, however as 508 cage mates of phenotyped mice had not been genotyped, we estimated relationship coefficients for all 3,149 mice based on both genotype and pedigree data.

We discarded the genotypes of 18 mice because they represented likely duplicates.

We infered the pair-wise relationship coefficient of the 1,227 non-genotyped mice with all 3,149 mice using the single-step or H matrix[63–65], which was constructed as (in the inverse scale;):

$$H^{-1} = A^{-1} + \begin{pmatrix} 0 & 0 \\ 0 & G^{(*)-1} + A_{gg}^{-1} \end{pmatrix}$$

where $A^{-1}$ is the inverse of pedigree-based Wright (1921) relationship matrix, and $A_{gg}$ is the pedigree-based relationship submatrix across genotyped individuals. The pedigree consisted of the genotyped individuals and their parents, some of them common across sibships. To warrant that both $G^{(*)}$ and $A^{-1}$ refer to the same base population and to ensure that $G^{(*)}$ is a valid covariance matrix, $G^{(*)}$ was constructed as [66, 67]

$$G^{(*)} = 0.95 \left( a + b \frac{MIM}{2 \sum (p_i q_i)} \right) + 0.05 A_{gg}$$
where a and b are constants (respectively 0.015 and 0.982) based on statistics of A and G so that average relationship and average inbreeding agree for both matrices. These constants correct for the random drift and the loss on heterozygosity across generations. Matrix M contains centered genotypes based on observed allelic frequencies p and q = 1 − p for each locus as in [68]. The matrix H was obtained by inversion of H⁻¹. Program pregsf90 ([69]; available at http://nce.ads.uga.edu/wiki/doku.php) was used to construct H. Using this matrix, all estimates of genetic variance refer to the pedigree base populations, i.e. the parents of the animals in the study. Because this is only one generation back in time, we can assume that this base genetic variance is concordant with the genetic variance in the phenotyped population.

**Processing of organismal phenotypes.** We considered a subset of the experimental variables recorded by Valdar *et al.*[22] as covariates for our analysis (S5 Table). Each phenotype was normalized using a covariate-aware Box-Cox transformation[49], and the covariates were fitted as fixed effects in the models used in our analyses.

**Re-processing of expression data.** Gene expression in the hippocampus was assessed by hybridization of mRNA to the Illumina Mouse WG-6 v1 BeadArray platform. The data, as originally pre-processed by Huang *et al.* had critical shortcomings: S10A and S10B Fig show that the higher correlation that can be seen in the raw data between probes mapping to the same gene compared to random probe pairs is lost in the pre-processed data from Huang *et al.* [24]. We re-processed the data as detailed in Section 2.6 of Krohn[70]. Briefly, the BeadArray images were imported into the Gene Expression module (V 1.6.0) of the Illumina GenomeStudio (V 2010.1) without invoking any data adjustment procedures. These raw data were transferred to R using the Bioconductor package lumi[71]. A cluster dendrogram enabled identification of four outlier samples, which were excluded from further pre-processing. Two thousand control probes scattered across the Illumina chip enabled subtraction of background noise levels via lumi. Standardisation of variance and mean were carried out via lumi's variance stabilising transformation and robust spline normalisation. Then a series of four filters were applied to eliminate unreliable probes: following the pipeline developed by Barbosa-Morais *et al.*[72], probes not rated as "perfect" (50 out of 50 nucleotides match the reference genome) or "good" (48 to 49 out of 50 matches to the genome) were removed; so were probes that mapped to more than one genomic location according to BLAST queries[73] and probes that included a SNP (NCBI dbSNP Build 137); finally, only probes detected in at least 5% of the mice at a 0.95 detection level (as per GenomeStudio) were retained. Following filter application, 15,736 of the original 47K probes, mapping to 11,996 genes, were kept. Finally, ComBat[74] was used to control for batch effects (the data were assayed over several months). S10C Fig shows that our processing of the data kept the biological correlation that exists between probes mapping to the same gene.

**Statistical models.** We considered a linear additive model,

$$y_f = X_f \underline{b} + a_{Df} + e_{Df} + \underline{Z_f a_S} + \underline{Z_f e_S} + \underline{W_f c} \quad (3)$$

Here, y_f is the phenotypic value of the focal mouse f, X_f is a row of the matrix X of covariate values and b is a column vector of corresponding estimated coefficients. a_Df is the additive direct genetic effects (DGE), also called breeding value, of f. Z_f is a row of the matrix Z that indicates cage mates (importantly Z_{i,i} = 0) and a_S the column vector of additive social genetic effects (SGE). e_D refers to direct environmental effects (DEE) and e_S to social environmental effects (SEE). W_f is a row of the matrix W that indicates cage assignment and c the column vector of cage effects
Because the genetic setup is not binary (B6 / D2 in the experiment with inbred strains) but individuals are genetically unique and related to varying degrees, we modeled the components in Model (3) as random effects, while accounting for covariance between different effects. The joint distribution of all random effects is defined as:

$$ \begin{bmatrix} \underline{a}_D \\ \underline{a}_S \\ \underline{e}_D \\ \underline{e}_S \\ \underline{c} \end{bmatrix} \sim \text{MVN} (0, \begin{bmatrix} \sigma_{A_D}^2 H_1 & \sigma_{A_D} H_2 & 0 & 0 & 0 \\ \sigma_{A_D} H_2^T & \sigma_{A_S}^2 H_3 & 0 & 0 & 0 \\ 0 & 0 & \sigma_{E_D}^2 I_1 & \sigma_{E_D} I_2 & 0 \\ 0 & 0 & \sigma_{E_D} I_2^T & \sigma_{E_S}^2 I_3 & 0 \\ 0 & 0 & 0 & 0 & \sigma_C^2 I_4 \end{bmatrix}) $$

and the phenotypic covariance is:

$$ \operatorname{cov}(y_i, y_j) = \sigma_{A_D}^2 H_{1i,j} + \sigma_{A_D} \{(H_2 Z^T)_{i,j} + (ZH_2^T)_{i,j}\} + \sigma_{A_S}^2 (ZH_3 Z^T)_{i,j} + \sigma_{E_D}^2 I_{1i,j} \\ + \sigma_{E_D} \{(I_2 Z^T)_{i,j} + (ZI_2^T)_{i,j}\} + \sigma_{E_S}^2 (ZI_3 Z^T)_{i,j} + \sigma_C^2 (WI_4 W^T)_{i,j} \quad (4) $$

$H_1$, $H_2$, and $H_3$ were calculated by dividing the (single-step) $H$ matrix by, respectively, sampleVar($H$), $\sqrt{sampleVar(H) \times sampleVar(ZHZ^T)}$, and sampleVar($ZHZ^T$), where sampleVar is the sample variance of the corresponding covariance matrix: suppose that we have a vector $\underline{x}$ of random variables with covariance matrix $M$, the sample variance of $M$ is calculated as

$$ \operatorname{sampleVar}(M) = \frac{\operatorname{Tr}(PMP)}{n-1} $$

$Tr$ denotes the trace, $n$ is the sample size, and $P = I - \frac{11^T}{n}$ [75, 76].

As a result of this scaling, $\sigma_{A_D}^2$, $\sigma_{A_D}$, and $\sigma_{A_S}^2$ directly reflect the relative contribution of DGE, the covariance between DGE and SGE, and SGE.

Similarly, $I_1$, $I_2$, and $I_3$ were calculated by dividing the identity matrix $I$ by, respectively, $\operatorname{sampleVar}(I) = 1$, $\sqrt{\operatorname{sampleVar}(I) \times \operatorname{sampleVar}(ZIZ^T)}$, and $\operatorname{sampleVar}(ZIZ^T)$, so that $\sigma_{E_D}^2$, $\sigma_{E_D}$, and $\sigma_{E_S}^2$ directly reflect the phenotypic variance explained by DEE, the covariance between DEE and SEE, and SEE.

$I_4$ was calculated by dividing the identity matrix $I$ by $\operatorname{sampleVar}(WIW^T)$.

See S1 Note for a discussion of how these models differ from published SGE models.

**Variance partitioning.** For each trait, Eq (4) was used to partition phenotypic variance into its individual components. All models were fitted using LIMIX[26, 27], which we extended to handle models with correlated variance components. The model was fitted using restricted maximum likelihood (REML). We report the contribution of DGE and SGE to phenotypic variation as the proportions of phenotypic variance explained by DGE and SGE, i.e. the ratio between the sample variance explained by DGE ($\sigma_{A_D}^2$) or SGE ($\sigma_{A_S}^2$) and the sample variance of the overall phenotype covariance matrix (defined by Eq 4, with all variance components estimated).

Standard errors on variance components are estimated as the square root of the corresponding diagonal entries of the inverse of the Fisher information, which is the negative expected value of the Hessian[77].

**Simulation experiments.** Phenotypes were simulated based on the genotypes and the cage relationships of the full set of 2,448 mice. Phenotypes were drawn from model (3). We note $\rho_{A_D}$ the correlation between $a_D$ and $a_S$: $\rho_{A_D} = \sigma_{A_D} / (\sigma_{A_D} \times \sigma_{A_S})$.

We present several sets of simulations:
* one based on the estimated contributions of all parameters averaged over all organismal phenotypes: $\sigma_{A_D}^2 = 15\%$, $\sigma_{A_S}^2 = 4\%$, $\rho_{A_{DS}} = 0.5$, $\sigma_{E_D}^2 = 52\%$, $\sigma_{E_D}^2 = 0$, $\rho_{E_{DS}} = 0$, $\sigma_C^2 = 22\%$ (S6B Fig and S7B Fig)

* one based on the estimated contribution of all parameters averaged over the six phenotypes with highest SGE: $\sigma_{A_D}^2 = 5\%$, $\sigma_{A_S}^2 = 27\%$, $\rho_{A_{DS}} = 0.92$, $\sigma_{E_D}^2 = 5\%$, $\sigma_{E_D}^2 = 0$, $\rho_{E_{DS}} = 0$, $\sigma_C^2 = 51\%$ (S6C Fig and S7C Fig)

* $\sigma_{A_D}^2$ taking values {0,10,20,30,40,50} with all other parameters set to their average estimated contribution across all phenotypes (so $\sigma_{A_S}^2 = 4\%$, $\rho_{A_{DS}} = 0.5$, $\sigma_{E_D}^2 = 52\%$, $\sigma_{E_S}^2 = 0$, $\rho_{E_{DS}} = 0$, $\sigma_C^2 = 22\%$) (S6A Fig and S7A Fig)

* $\sigma_{A_S}^2$ taking values {0,10,20,30,40,50} with all other parameters set to their average estimated contribution across all phenotypes (Fig 2C, S6A Fig and S7A Fig)

* $\rho_{A_{DS}}$ taking values {-0.8,-0.4,-0.2,0.2,0.4,0.8} with all other parameters set to their average estimated contribution across all phenotypes (S6A Fig and S7A Fig)

* $\sigma_{E_D}^2$ taking values {0,10,20,30,40,50} with all other parameters set to their average estimated contribution across all phenotypes (S6A Fig and S7A Fig)

* $\sigma_{E_S}^2$ taking values {0,10,20,30,40,50} with all other parameters set to their average estimated contribution across all phenotypes (S6A Fig and S7A Fig)

* $\rho_{E_{DS}}$ taking values {-0.8,-0.4,-0.2,0.2,0.4,0.8} with all other parameters set to their average estimated contribution across all phenotypes (S6A Fig and S7A Fig)

* $\sigma_C^2$ taking values {0,10,20,30,40,50} with all other parameters set to their average estimated contribution across all phenotypes (S6A Fig and S7A Fig)

500 simulations were drawn for each combination of parameters.

**Statistical significance tests.** Statistical significance was assessed using restricted likelihood ratio tests, comparing the full model in Eq (3) to a null model without the SGE term. Because both SGE and the covariance between DGE and SGE are constrained to zero in the null model, we calculated significance by comparing the log likelihood ratio statistics to the $\chi^2$ distribution with 2 degrees of freedom, which is conservative. Q values [58] were used for multiple testing correction, and obtained using the R package qvalue [59].

## Accession numbers

Gene expression data from the experiment with inbred strains are available from ArrayExpress E-MTAB-5276. Phenotype data for the same experiment are provided as S7 Table.

## Supporting Information

**S1 Note. Differences between our models and those published in the field of animal breeding.**
(DOCX)

**S1 Fig. Experimental design and analysis approach (experiment with inbred strains).** (A) When jointly analyzing B6 and D2 focal mice (for S2 Fig and columns “all mice” of S2 Table), DGE, SGE and their interaction can be detected. All other results were obtained by considering B6 focal mice (B) and D2 focal mice (C) separately. The only genetic parameter then is the
strain of the cage mate. Environmental effects are not shown.
(TIF)

**S2 Fig. Evidence for interactions between DGE and SGE in the experiment with inbred strains.** The contribution to phenotypic variance of DGE, SGE, and their interaction is shown. All organismal phenotypes are shown (see S1 Table for a description of each measure). They are plotted in order of decreasing contribution of the interaction component. Blue stars indicate phenotypes for which a model with interaction was selected based on the Akaike information criterion (AIC).
(TIF)

**S3 Fig. Consistent direction of social genetic effects across the anxiety phenotypes that are significantly affected by SGE in D2 focal mice (experiment with inbred strains).** Box plots for the six anxiety-related phenotypes that are significantly affected by SGE in D2 mice (P < 0.05, see Table 1). The raw data are shown, for all four types of cages (B6 focal mice with B6 cage mate (BB), B6 focal mice with D2 cage mate (BD), D2 focal mice with B6 cage mate (DB), and D2 focal mice with D2 cage mate (DD)).
* The sign of the correlation between anxiety and time spent not moving in the center of the activity chamber in the first minute of the test in unclear: greater time spent in the center of the apparatus usually indicates lower anxiety levels, but the measure may reflect the mouse freezing where it was dropped off in the apparatus-which would be a sign of high anxiety.
(EPS)

**S4 Fig. Group sizes in the outbred mice dataset.**
(EPS)

**S5 Fig. Consistent results using genomic and inferred relationship coefficients (outbred mice dataset).** Comparison of the contribution of DGE (A) and SGE (B) to organismal phenotypes, estimated either using genotype and pedigree data (see Methods; all mice included, x axis) or genotype data only (subset of mice that were in fully genotyped cages, y axis).
(EPS)

**S6 Fig. Need to account for SGE to obtain unbiased DGE estimates in the outbred mice dataset.** DGE estimates were evaluated using simulated phenotypes: the difference between estimated and simulated DGE is shown on the y-axis. Each boxplot summarizes the results from 500 simulations, and the whiskers extend to the 1st and 9th deciles of the distribution. (A) Phenotypes simulated with varying levels DGE, SGE, correlation between DGE and SGE, DEE, SEE, correlation between DEE and SEE, and cage effects (top to bottom panels). (B) Phenotypes simulated based on the estimated contributions of all parameters averaged over all organismal phenotypes. (C) Phenotypes simulated based on the estimated contributions of all parameters averaged over the six organismal phenotypes with strongest SGE. The boxplots in blue correspond to results from the full model with DGE, SGE, SEE and cage effects, those in green correspond to the model with DGE and cage effects only, and those in red correspond to the model with DGE only.
(EPS)

**S7 Fig. Unbiased SGE estimates in the outbred mice dataset.** SGE estimates were evaluated using simulated phenotypes: the difference between estimated and simulated SGE is shown on the y-axis. Each boxplot summarizes the results from 500 simulations, and the whiskers extend to the 1st and 9th deciles of the distribution. (A) Phenotypes simulated with varying levels DGE, SGE, correlation between DGE and SGE, DEE, SEE, correlation between DEE and SEE, and cage effects (top to bottom panels). (B) Phenotypes simulated based on the estimated
contributions of all parameters averaged over all organismal phenotypes. (C) Phenotypes simulated based on the estimated contributions of all parameters averaged over the six organismal phenotypes with strongest SGE. Estimates were obtained from the full model with DGE, SGE, and cage effects.

(EPS)

**S8 Fig. Above-average relatedness of cage mates compared to non-cage mates in the outbred dataset.**

(EPS)

**S9 Fig. Identification of outlier gene expression profiles using principal component analysis in the experiment with inbred strains.** (A) Variance explained by the principal components. (B) Scatterplot of the first and second principal components. Samples that were excluded from all analyses (outliers) are highlighted.

(EPS)

**S10 Fig. Improved quantification of gene expression in the hippocampus following re-processing of the data (outbred mice dataset).** The quality of the quantification is evaluated based on the correlation between probes that target the same gene (orange line) and probes that target different genes (blue line). The empirical cumulative function of these correlations is shown, so a smaller area under the curve is indicative of greater correlations.

(A) Raw data from Huang *et al.*²⁶: expression values from pairs of probes that map to the same gene are more highly correlated than expression values from pairs of probes that map to different genes.

(B) Pre-processed data from Huang *et al.*²⁶: same-gene probes are only as correlated as pairs of probes mapping to different genes, suggesting a loss of biological signal.

(C) Data pre-processed by us: same-gene probes are more highly correlated than different-gene probes, suggesting biological signal is retained.

(EPS)

**S11 Fig. Independency of residual variability on group size (outbred mice dataset).** The residuals from a model with covariates fitted as fixed effects only (y axis) were plotted as a function of group size (x axis). Only the six organismal phenotypes with strongest SGE are shown.

(EPS)

**S12 Fig. Consistent results using mice in cages of 4,5 or 6 mice only (outbred mice dataset).** Comparison of the contribution of DGE (A) and SGE (B) to organismal phenotypes, estimated either using all mice (x axis) or only those mice that were in groups of size 4, 5 or 6 (y axis).

(EPS)

**S1 Table. Organismal phenotypes collected in the experiment with inbred strains.**

(XLS)

**S2 Table. SGE significance and effect sizes for all organismal phenotypes measured in the experiment with inbred strains.**

(XLSX)

**S3 Table. SGE significance for non-residualized organismal phenotypes measured in the experiment with inbred strains.**

(XLSX)
**S4 Table.** SGE significance and effect sizes for all genes expressed in the prefrontal cortex (experiment with inbred strains).

(XLS)

**S5 Table.** Description of all organismal phenotypes in the outbred mice dataset, and significance and effect sizes of SGE.

(XLSX)

**S6 Table.** SGE significance and effect sizes for all genes expressed in the hippocampus (outbred mice dataset).

(XLSX)

**S7 Table.** Raw phenotype data collected on inbred strains.

(XLS)

## Acknowledgments

We thank the High-Throughput Genomics Group at the Wellcome Trust Centre for Human Genetics for the generation of the RNA Sequencing data.

## Author Contributions

**Conceptualization:** AB OS.

**Formal analysis:** AB FPC JK AL.

**Investigation:** AB MKM CJB JFI JC.

**Methodology:** AB MKM FPC AL RWW OS.

**Resources:** AB JML RWW OS.

**Writing – original draft:** AB OS.

**Writing – review & editing:** AB MKM FPC JFI JC JK AL RWW OS.

## References

1. Moore AJ, Brodie ED III, Wolf JB. Interacting phenotypes and the evolutionary process: I. Direct and indirect genetic effects of social interactions. Evolution; international journal of organic evolution. 1997:1352–62.

2. Bergsma R, Kanis E, Knol EF, Bijma P. The contribution of social effects to heritable variation in finishing traits of domestic pigs (Sus scrofa). Genetics. 2008; 178(3):1559–70. Epub 2008/02/05. PubMed Central PMCID: PMC2391867. doi: 10.1534/genetics.107.084236 PMID: 18245326

3. Wolf JB, Wade MJ. What are maternal effects (and what are they not)? Philosophical Transactions of the Royal Society B: Biological Sciences. 2009; 364(1520):1107–15.

4. Wilson A, Coltman D, Pemberton J, Overall A, Byrne K, Kruuk L. Maternal genetic effects set the potential for evolution in a free-living vertebrate population. Journal of evolutionary biology. 2005; 18(2):405–14. doi: 10.1111/j.1420-9101.2004.00824.x PMID: 15715846

5. de Albuquerque G. Estimates of direct and maternal genetic effects for weights from birth to 600 days of age in Nelore cattle. Journal of Animal breeding and Genetics. 2001; 118(2):83–92.

6. Department of Biological Sciences University of South Carolina Timothy A. Mousseau Associate Professor C, University CWFAPF. Maternal Effects As Adaptations: Oxford University Press, USA; 1998.

7. Strandberg E, Jacobsson J, Saetre P. Direct genetic, maternal and litter effects on behaviour in German shepherd dogs in Sweden. Livestock Production Science. 2005; 93(1):33–42.

8. Champagne FA, Meaney MJ. Stress during gestation alters postpartum maternal care and the development of the offspring in a rodent model. Biological psychiatry. 2006; 59(12):1227–35. doi: 10.1016/j.biopsych.2005.10.016 PMID: 16457784
9. Brinker T, Bijma P, Visscher J, Rodenburg TB, Ellen ED. Plumage condition in laying hens: genetic parameters for direct and indirect effects in two purebred layer lines. Genetics, selection, evolution: GSE. 2014; 46:33. Epub 2014/06/03. PubMed Central PMCID: PMC4073196. doi: 10.1186/1297-9686-46-33 PMID: 24885199

10. Peeters K, Eppink TT, Ellen ED, Visscher J, Bijma P. Indirect genetic effects for survival in domestic chickens (Gallus gallus) are magnified in crossbred genotypes and show a parent-of-origin effect. Genetics. 2012; 192(2):705-13. doi: 10.1534/genetics.112.142554 PMID: 22851648

11. Duijvesteijn N, Knol E, Bijma P. Direct and associative effects for androstenone and genetic correlations with backfat and growth in entire male pigs. Journal of animal science. 2012; 90(8):2465-75. doi: 10.2527/jas.2011-4625 PMID: 22367075

12. Arango J, Misztal I, Tsuruta S, Culbertson M, Herring W. Estimation of variance components including competitive effects of Large White growing gilts. Journal of animal science. 2005; 83(6):1241-6. PMID: 15890801

13. Khaw H. Indirect genetic effects for harvest weight in Nile tilapia (Oreochromis niloticus). Age. 2014; 6330(166):58.97.

14. Petfield D, Chenoweth SF, Rundle HD, Blows MW. Genetic variance in female condition predicts indi-rect genetic variance in male sexual display traits. Proceedings of the National Academy of Sciences. 2005; 102(17):6045-50.

15. Alemu SW, Bijma P, Møller SH, Janss L, Berg P. Indirect genetic effects contribute substantially to heri-table variation in aggression-related traits in group-housed mink (Neovison vison). Genet Sel Evol. 2014; 46:30. doi: 10.1186/1297-9686-46-30 PMID: 24884874

16. Wilson A, Morrissey M, Adams M, Walling C, Guinness F, Pemberton J, et al. Indirect genetics effects and evolutionary constraint: an analysis of social dominance in red deer, Cervus elaphus. Journal of evolutionary biology. 2011; 24(4):772-83. doi: 10.1111/j.1420-9101.2010.02212.x PMID: 21288272

17. Sartori C, Mantovani R. Indirect genetic effects and the genetic bases of social dominance: evidence from cattle. Heredity. 2013; 110(1):3-9. doi: 10.1038/hdy.2012.56 PMID: 22968222

18. Peirce JL, Lu L, Gu J, Silver LM, Williams RW. A new set of BXD recombinant inbred lines from advanced intercross populations in mice. BMC genetics. 2004; 5(1):7.

19. Bredy TW, Wu H, Crego C, Zellhoefer J, Sun YE, Barad M. Histone modifications around individual BDNF gene promoters in prefrontal cortex are associated with extinction of conditioned fear. Learning & memory. 2007; 14(4):268-76.

20. Caldji C, Diorio J, Anisman H, Meaney MJ. Maternal behavior regulates benzodiazepine/GABAA recep-tor subunit expression in brain regions associated with fear in BALB/c and C57BL/6 mice. Neuropsychopharmacology. 2004; 29(7):1344-52. doi: 10.1038/sj.npp.1300436 PMID: 15085086

21. Kerns RT, Ravindranathan A, Hassan S, Cage MP, York T, Sikela JM, et al. Ethanol-responsive brain region expression networks: implications for behavioral responses to acute ethanol in DBA/2J versus C57BL/6J mice. The Journal of neuroscience. 2005; 25(9):2255-66. doi: 10.1523/JNEUROSCI.4372-04.2005 PMID: 15745951

22. Valdar W, Solberg LC, Gauguier D, Burnett S, Klenerman P, Cookson WO, et al. Genome-wide genetic association of complex traits in heterogeneous stock mice. Nature genetics. 2006; 38(8):879-87. doi: 10.1038/ng1840 PMID: 16832355

23. Huang G-J, Smith AL, Gray DH, Cosgrove C, Singer BH, Edwards A, et al. A genetic and functional rela-tionship between T cells and cellular proliferation in the adult hippocampus. 2010.

24. Huang G-J, Shifman S, Valdar W, Johannesson M, Yalcin B, Taylor MS, et al. High resolution mapping of expression QTLs in heterogeneous stock mice in multiple tissues. Genome research. 2009; 19(6):1133-40. doi: 10.1101/gr.088120.108 PMID: 19376938

25. Bijma P, Muir WM, Ellen ED, Wolf JB, Van Arendonk JA. Multilevel selection 2: estimating the genetic parameters determining inheritance and response to selection. Genetics. 2007; 175(1):289-99. doi: 10.1534/genetics.106.062729 PMID: 17110493

26. Lippert C, Casale FP, Rakitsch B, Stegle O. LIMIX: genetic analysis of multiple traits. bioRxiv. 2014.

27. Casale FP, Rakitsch B, Lippert C, Stegle O. Efficient set tests for the genetic analysis of correlated traits. Nature methods. 2015.

28. Bijma P. Estimating indirect genetic effects: precision of estimates and optimum designs. Genetics. 2010; 186(3):1013-28. doi: 10.1534/genetics.110.120493 PMID: 20713688

29. Gardner M, Steinberg L. Peer influence on risk taking, risk preference, and risky decision making in ado-lescence and adulthood: an experimental study. Developmental psychology. 2005; 41(4):625-35. doi: 10.1037/0012-1649.41.4.625 PMID: 16060809
30. Groh DR, Jason LA, Keys CB. Social network variables in alcoholics anonymous: a literature review. Clinical psychology review. 2008; 28(3):430–50. PubMed Central PMCID: PMC2289871. doi: 10.1016/j.cpr.2007.07.014 PMID: 17719158

31. Haeffel GJ, Hames JL. Cognitive Vulnerability to Depression Can Be Contagious. Clinical Psychological Science. 2014; 2(1):75–85.

32. Kumpulainen K. Psychiatric conditions associated with bullying. International journal of adolescent medicine and health. 2008; 20(2):121–32. PMID: 18714551

33. Russek LG, Schwartz GE. Feelings of parental caring predict health status in midlife: A 35-year follow-up of the Harvard Mastery of Stress Study. J Behav Med. 1997; 20(1):1–13. PMID: 9058175

34. Hutson J, Niall M, Evans D, Fowler R. Effect of salivary glands on wound contraction in mice. 1979.

35. Detillion CE, Craft TK, Glasper ER, Prendergast BJ, DeVries AC. Social facilitation of wound healing. Psychoneuroendocrinology. 2004; 29(8):1004–11. doi: 10.1016/j.psyneuen.2003.10.003 PMID: 15219651

36. Ebrecht M, Hextall J, Kirtley L-G, Taylor A, Dyson M, Weinman J. Perceived stress and cortisol levels predict speed of wound healing in healthy male adults. Psychoneuroendocrinology. 2004; 29(6):798–809. doi: 10.1016/S0306-4530(03)00144-6 PMID: 15110929

37. Manolio TA, Collins FS, Cox NJ, Goldstein DB, Hindorff LA, Hunter DJ, et al. Finding the missing heritability of complex diseases. Nature. 2009; 461(7265):747–53. Epub 2009/10/09. PubMed Central PMCID: PMC2831613. doi: 10.1038/nature08494 PMID: 19812666

38. Eaves L. A model for sibling effects in man. Heredity. 1976; 36(2):205–14. PMID: 1063734

39. Silventoinen K, Sammalisto S, Perola M, Boomsma DI, Cornes BK, Davis C, et al. Heritability of adult body height: a comparative study of twin cohorts in eight countries. Twin research. 2003; 6(05):399–408.

40. Hyttinen V, Kaprio J, Kinnunen L, Koskenvuo M, Tuomilehto J. Genetic Liability of Type 1 Diabetes and the Onset Age Among 22,650 Young Finnish Twin Pairs A Nationwide Follow-Up Study. Diabetes. 2003; 52(4):1052–5. PMID: 12663480

41. Wirdefeldt K, Gatz M, Reynolds CA, Prescott CA, Pedersen NL. Heritability of Parkinson disease in Swedish twins: a longitudinal study. Neurobiology of aging. 2011; 32(10):1923. e1–e8.

42. Conley D, Siegal ML, Domingue BW, Harris KM, McQueen MB, Boardman JD. Testing the key assumption of heritability estimates based on genome-wide genetic relatedness. Journal of human genetics. 2014; 59(6):342–5. doi: 10.1038/jhg.2014.14 PMID: 24599117

43. Muir WM. Incorporation of competitive effects in forest tree or animal breeding programs. Genetics. 2005; 170(3):1247–59. doi: 10.1534/genetics.104.035956 PMID: 15911590

44. Wolf JB, Brodie ED III, Cheverud JM, Moore AJ, Wade MJ. Evolutionary consequences of indirect genetic effects. Trends in ecology & evolution. 1998; 13(2):64–9.

45. Dawkins R. The Extended Phenotype: The Long Reach of the Gene. Oxford University Press; 1999.

46. Kema IP, Schellings AM, Hoppenbrouwers CJ, Rutgers HM, de Vries EG, Muskiet FA. High performance liquid chromatographic profiling of tryptophan and related indoles in body fluids and tissues of carcinoid patients. Clinica chimica acta; international journal of clinical chemistry. 1993; 221(1–2):143–58. Epub 1993/11/30. PMID: 7512001

47. Gamache P, Ryan E, Svendsen C, Murayama K, Acworth IN. Simultaneous measurement of monamines, metabolites and amino acids in brain tissue and microdialysis perfusates. Journal of chromatography. 1993; 614(2):213–20. Epub 1993/05/05. PMID: 7686176

48. Box GE, Cox DR. An analysis of transformations. Journal of the Royal Statistical Society Series B (Methodological). 1964:211–52.

49. Venables W, Ripley B. Modern applied statistics with S. Fourth edition ed: Springer; 2002.

50. Keane TM, Goodstadt L, Danecek P, White MA, Wong K, Yalcin B, et al. Mouse genomic variation and its effect on phenotypes and gene regulation. Nature. 2011; 477(7364):289–94. Epub 2011/09/17. PubMed Central PMCID: PMC3276836. doi: 10.1038/nature10413 PMID: 21921910

51. Munger SC, Raghupathy N, Choi K, Simons AK, Gatti DM, Hinerfeld DA, et al. RNA-Seq alignment to individualized genomes improves transcript abundance estimates in multiparent populations. Genetics. 2014; 198(1):59–73. Epub 2014/09/23. PubMed Central PMCID: PMC4174954. doi: 10.1534/genetics.114.165886 PMID: 25236449

52. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome biology. 2013; 14(4): R36. Epub 2013/04/27. PubMed Central PMCID: PMC4053844. doi: 10.1186/gb-2013-14-4-r36 PMID: 23618408

53. Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015; 31(2):166–9. Epub 2014/09/28. PubMed Central PMCID: PMC4287950. doi: 10.1093/bioinformatics/btu638 PMID: 25260700
54. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology. 2014; 15(12):550. Epub 2014/12/18. PubMed Central PMCID: PMC4302049. doi: 10.1186/s13059-014-0550-8 PMID: 25516281

55. Stegle O, Parts L, Piipari M, Winn J, Durbin R. Using probabilistic estimation of expression residuals (PEER) to obtain increased power and interpretability of gene expression analyses. Nature protocols. 2012; 7(3):500-7. doi: 10.1038/nprot.2011.457 PMID: 22343431

56. Lassmann T, Hayashizaki Y, Daub CO. SAMStat: monitoring biases in next generation sequencing data. Bioinformatics. 2011; 27(1):130-1. Epub 2010/11/20. PubMed Central PMCID: PMC3008642. doi: 10.1093/bioinformatics/btq614 PMID: 21088025

57. Andrews S. FastQC: A quality control tool for high throughput sequence data. Available from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

58. Storey JD. The positive false discovery rate: a Bayesian interpretation and the q-value. Annals of statistics. 2003;2013-35.

59. Bass JDSwcfAJ DAaRD. qvalue: Q-value estimation for false discovery rate control. R package version 220. 2015.

60. Alexa A, Rahnenfuhrer J. topGO: enrichment analysis for gene ontology. R package version. 2010; 2 (0).

61. Alexa A, Rahnenführer J, Lengauer T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 2006; 22(13):1600-7. doi: 10.1093/bioinformatics/btl140 PMID: 16606683

62. Demarest K, Koyner J, McCaughran Jr J, Cipp L, Hitzemann R. Further characterization and high-resolution mapping of quantitative trait loci for ethanol-induced locomotor activity. Behavior genetics. 2001; 31(1):79-91. PMID: 11529277

63. Legarra A, Aguilar I, Misztal I. A relationship matrix including full pedigree and genomic information. J Dairy Sci. 2009; 92(9):4656-63. doi: 10.3168/jds.2009-2061 PMID: 19700729

64. Christensen OF, Lund MS. Genomic prediction when some animals are not genotyped. Genetics Selection Evolution. 2010; 42(1):1.

65. Aguilar I, Misztal I, Johnson D, Legarra A, Tsuruta S, Lawlor T. Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score. Journal of Dairy Science. 2010; 93(2):743-52. doi: 10.3168/jds.2009-2730 PMID: 20105546

66. Vitezica ZG, Aguilar I, Misztal I, Legarra A. Bias in genomic predictions for populations under selection. Genet Res (Camb). 2011; 93(5):357-66.

67. Christensen OF, Madsen P, Nielsen B, Ostersen T, Su G. Single-step methods for genomic evaluation in pigs. Animal. 2012; 6(10):1565-71. doi: 10.1017/S1751731112000742 PMID: 22717310

68. VanRaden P. Efficient methods to compute genomic predictions. Journal of dairy science. 2008; 91 (11):4414-23. doi: 10.3168/jds.2007-0980 PMID: 18946147

69. Aguilar I, Misztal I, Legarra A, Tsuruta S. Efficient computation of the genomic relationship matrix and other matrices used in single-step evaluation. Journal of Animal Breeding and Genetics. 2011; 128 (6):422-8. doi: 10.1111/j.1439-0388.2010.00912.x PMID: 22059575

70. Krohn J. Genes contributing to variation in fear-related behaviour: Oxford, UK; 2013.

71. Du P, Kibbe WA, Lin SM. Iumi: a pipeline for processing Illumina microarray. Bioinformatics. 2008; 24 (13):1547-8. Epub 2008/05/10. doi: 10.1093/bioinformatics/btn224 PMID: 18467348

72. Barbosa-Morais NL, Dunning MJ, Samarajiwa SA, Darot JF, Ritchie ME, Lynch AG, et al. A re-annotation pipeline for Illumina BeadArrays: improving the interpretation of gene expression data. Nucleic acids research. 2010; 38(3):e17. Epub 2009/11/20. PubMed Central PMCID: PMC2817484. doi: 10.1093/nar/gkp942 PMID: 19923232

73. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic acids research. 1997; 25 (17):3389-402. Epub 1997/09/01. PubMed Central PMCID: PMC146917. PMID: 9254694

74. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007; 8(1):118-27. doi: 10.1093/biostatistics/kxj037 PMID: 16632515

75. Kang HM, Sul JH, Service SK, Zaitlen NA, Kong S-y, Freimer NB, et al. Variance component model to account for sample structure in genome-wide association studies. Nature genetics. 2010; 42(4):348-54. doi: 10.1038/ng.548 PMID: 20208533

76. Searle SR. Matrix algebra useful for statistics: Wiley; 1982.

77. Murphy KP. Machine Learning: A Probabilistic Perspective: The MIT Press; 2012. 1096 p.