Variation in the human immune system is largely driven by non-heritable influences

Estimating heritable and non-heritable influences

Heritability for each parameter was estimated by comparing observed MZ and DZ covariance matrices to the expected values based on a structural equation model that partitioned the observed variance into three components, heritable (A), shared (C) and unique (E) non-heritable factors. This model is based on the assumptions that, I) heritable factors correlate perfectly between MZ twins (rMZ=1) but only to 50% between DZ twins (rDZ=0.5), and II) that shared non-heritable influences are equally similar (rMZ = rDZ) between MZ and DZ twin pairs (Experimental procedures, section 7). For each measurement, we subtracted the technical variance estimate from the E-component prior to normalization to correct for noise (Experimental procedures, section 8). We also corrected all measurements for the effects of age (Dorshkind et al., 2009) and gender (Furman et al., 2014), by regressing out such effects and using only residual variance for estimating heritability. Finally, we performed jack-knife bootstrapping tests to obtain 95% confidence intervals (Experimental procedures, section 7). Importantly, as our model estimates heritability by comparing MZ and DZ twins, heritable influences include genomic and shared epigenetic traits (Bell and Spector, 2011), and non-heritable influences include environmental factors and stochastic epigenetic changes (Fraga et al., 2005).

We first ran a simulation experiment to verify that our cohort size of 210 twins (78 MZ and 27 DZ pairs) would be enough to test our hypothesis that most immunological traits are explained more by non-heritable than heritable influences. We found this to be the case, and we estimate 20% heritability to be our detection limit, under which we cannot distinguish small heritable influences from zero (Figure S1).

Most cell population frequencies and serum proteins are dominated by non-heritable influences

While it is well known that the frequencies of different types of immune cells in blood often vary widely between individuals, in most cases it is not known how much of this can be attributed to heritable or non-heritable factors respectively. To address this question, we used antibodies against cell surface markers to quantify 95 different cell subset frequencies, but used the 72 most non-redundant ones and estimated the influence of heritable and non-heritable factors on their variation (Experimental Procedures, section 3). Among these a few had very strong influences from heritable factors, especially naïve, CD27⁺ and central memory CD4⁺ T-cells (Figure 1B and Table S3), but for most, non-heritable influences were clearly dominant. In fact, for 61% of all cell populations, the influence of heritable factors was undetectable (<20% of the total variation)(Figure 1B and Table S3). This was true of both adaptive (T- and B-cells) and innate cell types (granulocytes, monocytes and NK-cells).

Serum cytokines and chemokines also have important functions as immune mediators and biomarkers of disease (Villeda and Wyss-Coray, 2013) and thus we measured 51 serum proteins, but eliminated eight that were often at or below the limits of detection (Experimental Procedures, section 5). This left 24 cytokines, 10 chemokines, 6 growth factors and 3 other serum proteins for which we estimated the influences of heritable and non-heritable factors (Figure 1C, Table S4). Some cytokines were particularly heritable, such as IL-12p40 (Figure 1C, Table S4). Interestingly, variants in the IL12B gene that contribute to the IL-12p40 protein have been associated with immune-mediated diseases such as psoriasis (Nair et al., 2009) and asthma (Morahan et al., 2002). In the latter condition, the susceptibility locus was also associated with a reduced serum concentration of IL-12p40 (Morahan et al., 2002). For many other measurements, such as IL-10 and a group of chemokines, the heritable influence was low (Figure 1C, Table S4).

Homeostatic cytokine responses are largely heritable, while most other cell responses are highly non-heritable

Since these serum proteins often regulate immune cells, we assessed the responses of eight different cell populations stimulated in vitro with seven different cytokines for the phosphorylation of three important and used transcription factors STAT1, 3 and 5 using phospho-specific antibodies in flow cytometry (Krutzik and Nolan, 2006). We performed a total of 192 different measurements, but focused on the 24 baseline measurements and the 65 strongest induced responses (Experimental Procedures, section 4). Baseline measurements were generally driven by non-heritable factors, with possible minor contributions from heritable factors (Figure 2A). The important homeostatic cytokines IL-2 and IL-7, known to stimulate the proliferation and differentiation of T-cells, were found to induce STAT5 phosphorylation in both CD4⁺ and CD8⁺ T-cell populations and these responses were highly heritable (Figure 2A, Table S5). In contrast, most signaling responses such as interferon-induced STAT1 phosphorylation, and in particular the IL-6, IL-21 and regulatory IL-10 induced phosphorylation of STAT3, were dominated by non-heritable influences (Figure 2A, Table S5). In total, 69% of all signaling responses had no detectable heritable influence (e.g. < 20%)(Figure 2A, Table S5). This lack of heritability was not related to the strength of responses and explained by a bias towards weak and variable responses (Figure S2).

Open in a separate window

Figure 2

Heritable factors explain only a fraction of the variation for most immune measurements

(A) Heritability estimates for immune cell signaling states, upon stimulation with the indicated cytokines. Only unstimulated controls and induced responses > 1.5-fold are shown, see also Table S5. (B) The overall distribution of heritability for all 204 measurements. (C) The maximum number of measurements with heritability < 0.2 across 1000 synthetic datasets with the same MZ/DZ-ratio as in our twin cohort is < 40%, significantly less than our results of 58% of measurements with heritability < 20% (grey bar), p-value < 0.001.

Non-heritable influences are major factors determining immune variation

Taken together, these results show that variation in blood cell frequencies and functions and soluble factors is largely driven by non-heritable factors, with 58% of all measurements having less than 20% of their total variance explained by heritable influences (Figure 2B). There was no relationship between the absolute degree of measurement variability in the cohort and estimated heritability (Figure S3), and we could also rule out any underestimation of heritability due to the skewed ratio of MZ/DZ twin pairs in our cohort by a resampling test. Briefly, by creating 1000 synthetic datasets with uniform heritability and the same MZ/DZ ratio as in our cohort, we estimated heritability and found that none of the 1000 datasets ever had more than 40% of measurements with an estimated heritability < 0.2 (p-value < 0.001) (Figure 2C), thus suggesting that the low heritability estimates obtained are not a result of study design or overall measurement variation in the cohort.

Heritable and non-heritable measurements are interrelated in the immune network

Our analysis also allowed us to analyze the interrelationships between the different components of the immune system. To construct such a network model, we calculated a precision matrix derived from a Spearman correlation matrix (Liu et al., 2012). A matrix of this type captures partial correlations between variables and avoids spurious, indirect interactions that might occur in standard correlation analyses. By penalizing non-zero entries in this matrix (Friedman et al., 2008), we could pursue what is referred to as a sparse (rather than dense) network model, making it more interpretable. After removing unconnected nodes, and validating the edges by a permutation test (Experimental procedures, section 9), this model consists of 126 nodes and 142 edges (Figure 3A and Table S6). An interactive version is available online (http://www.brodinlab.com/twins.html). We found that heritable nodes (yellow) were generally connected to non-heritable nodes (purple) throughout the network (Figure 3A). One example shows how the weakly heritable cytokine IL-10 and CD161⁻ CD45RA⁺ regulatory T-cells are connected to the strongly heritable frequency of naïve CD4⁺ T-cells (Figure 3B). We found that all hubs in the network were dominated by non-heritable influences, like the network as a whole, showing that heritable factors are not isolated by themselves, but are buffered by connected non-heritable ones (Figure 3A–B). This may explain why the many gene polymorphisms found (for example CTLA4 (Gregersen et al., 2009)), outside of the HLA locus that are associated with immune mediated disease, only explain a small fraction of the total disease risk (Todd, 2010).

Open in a separate window

Figure 3

Network of dependencies between immune measurements

(A) Undirected network model of the healthy human immune system showing 126 nodes (measurements), connected by 142 undirected edges illustrating conditional measurement dependencies. Nodes are colored by their estimated heritability and sized by their number of edges.(B) Sub-network exemplifying direct relationships between heritable and non-heritable nodes. Solid edges represent positive relationships and dashed edges represent negative relationships. The edge weight represents the strength of relationships. See also Table S6.

With age, genetically identical twins diverge as a consequence of non-heritable influences

As a major source of non-heritable influence is likely to be environmental, particularly microbial exposure, we reasoned that such influences would increase with time. To this end, we compared twin-twin correlations for all immune measurements between the oldest (> 60 years, median: 72 years) and the youngest (< 20 years, median: 13 years) MZ pairs in our cohort. Here we also note that twins in the younger (<20 years) cohort are in most cases living together, while the older (>60 years) twins, have lived apart for decades, so concordance can also be a result of either shared environment and/or shared exposure, in addition to genetic similarity. For several cell population frequencies, we found much reduced correlations with age (Figure 4A). In the most striking example, the frequency of eosinophils between the youngest MZ twins correlated very strongly at 0.72, but was highly uncorrelated at 0.13 between the oldest MZ twin pairs (Figure 4A). Similarly, several serum proteins showed remarkably reduced correlations between older as compared to younger MZ twins (Figure 4B). In particular the chemokine CXCL10/IP10 showed a strong correlation (0.79) between the youngest MZ twins but was greatly reduced (0.18) in the older MZ twins (Figure 4B). Similar patterns were found for many cell signaling responses (data not shown), suggesting that this immune divergence between genetically identical twins with age is a common phenomenon, consistent with a major role for environmental exposure in driving variation, although some epigenetic changes could also contribute (Fraga et al., 2005).

Open in a separate window

Figure 4

Increased variability in the immune system with age

(A) Twin-twin correlations (Spearman’s rank) for all cell frequencies within the youngest MZ twin pairs (≤ 20yrs, median: 13.5, n=25 pairs), and the oldest MZ twin pairs (≥ 60yrs, median: 72yrs, n=16). (B) Twin-twin correlations (Spearman’s rank) for all serum protein concentrations within the youngest MZ twin pairs (≤ 20yrs, median: 13yrs, n=26), and the oldest MZ twin pairs (≥ 60yrs, median: 73yrs, n=13).

Cytomegalovirus infection has a broad influence on immune variation

As we postulate that microbial exposure is a likely driver of immune variation with age, a particularly interesting example is cytomegalovirus (CMV), a lifelong viral infection which has striking effects on the immune phenotypes of both humans and rhesus macaques (Sylwester et al., 2005). In our twin cohort, 16 MZ pairs were discordant for CMV seropositivity, and we compared their twin-twin correlations for all measurements to those of 26 CMV concordant (negative) MZ-pairs. Here we found that the CMV discordant MZ twins showed greatly reduced correlations for many immune cell frequencies such as effector CD8⁺ and gamma-delta T-cells (Figure 5A), as compared to CMV negative MZ-twins. The same was true for cell signaling responses, especially in response to IL-10 and IL-6 stimulation (Figure 5B), as well as the concentrations of these same cytokines in serum (Figure 5C). In general, the influence of CMV was very broad, affecting 119 of all 204 measurements (58%), dispersed throughout the immune network (Figure 5D), and illustrating how at least one type of microbial exposure can dramatically modulate the overall immune profile of healthy individuals.

Open in a separate window

Figure 5

Broad non-heritable influences in the healthy immune system

(A) Twin-twin correlations (Spearman’s rank) for all cell frequency measurements made in CMV concordant negative (neg./neg.) MZ twin pairs (n=26 pairs), and CMV discordant (pos./neg.) MZ twin pairs (n=16 pairs). (B) Twin-twin correlations (Spearman’s rank) for cell signaling responses to cytokine stimulation, and (C) serum protein measurements between CMV neg/neg and CMV pos/neg MZ twin pairs. (D) 58% of all 126 nodes in the immune network model with reduced correlations in CMV pos/neg as compared to CMV neg/neg MZ-twin pairs.

2. Blood sampling, PBMC preparation and zygosity testing

Blood samples were collected in heparinized vacutainer tubes by venipuncture at day 0 (and day +28 for HAI responses after seasonal influenza vaccination) at the Clinical and Translational Research Unit, Stanford University Hospital. Whole blood collected in sodium heparin tubes was either analyzed immediately using our whole-blood flow cytometry protocol below or enriched for PBMCs using 15ml of Ficoll-Paque PLUS (GE Health care, Inc. Piscataway, NJ) and frozen at −80° C over night, transferred to liquid nitrogen and stored until further analysis. Zygosity was determined by comparing 384 SNP-loci using a discriminatory DNA-polymerase and ligase assay (GoldenGate genotyping assay, Illumina, Inc, San Diego, CA) and were performed by IGenix (Bainbridge Island, WA). Twins were considered fraternal if similarities in DNA markers were below 99.0%. Twins of the same pair was almost exclusively analyzed in the same experimental batch, irrespective of technology used in order to minimize technical variation.

3. Immune cell phenotyping by mass cytometry and flow cytometry

All experiments were performed by the Human Immune Monitoring Center at Stanford University. For years 2010 and 2011, two million thawed PBMCs were stained without prior resting or fixation using a panel of 26 different metal-tagged probes to surface antigens and DNA (Data S2A). After repeated washes, cells were analyzed by mass cytometry (CyTOF^™, Fluidigm, South San Francisco, CA) with the following instrument settings: High resolution mode analysis with cell length set to 10–75 pushes, a lower de-convolution threshold of 10 and instrument run in dual-count detection mode and noise reduction turned off. FCS3.0 files were manually analyzed using FlowJo v9.3 (TreeStar, Inc. Ashland OR) as shown in Figure S4. Year 2009 PBMC samples were similarly processed, but instead analyzed using a set of custom-made Lyoplates (BD Biosciences, San Jose, CA) covering 7 different antibody panels of fluorescently labeled antibodies (Data S2B–C). These samples were acquired using a LSRII flow cytometer (BD Biosciences) and similarly analyzed manually using the same FlowJo v9.3 (TreeStar, Inc.).

4. Immune cell signaling experiments

PBMCs were thawed in warm media, washed twice and resuspended at 0.5×10⁶ viable cells/microliter. 200 microliter of cells were plated per well in 96-well deep-well plates. After resting for 1 hour at 37° C, cells were, stimulated by addition of 50 microliter solutions of cytokines: IFNα, IFNγ, IL-6, IL-7, IL-10, IL-2, or IL-21 (Data S2D) respectively and incubated at 37°C for 15 minutes. Cells were then immediately fixed in 1.6% paraformaldehyde, permeabilized with 100% cold methanol, and kept at −80° C overnight. Each well was bar-coded by a unique combination of Pacific Orange and Alexa-750 dye concentrations (Invitrogen/Life Technologies, Inc.), the cells were washed with FACS buffer (PBS supplemented with 2% FBS and 0.1% sodium azide), and stained with our phospho-Flow antibody panel (Data S2E). Finally cells were washed and resuspended in FACS buffer and 100,000 cells per stimulation condition were collected using DIVA 6.0 software on an LSRII flow cytometer (BD Biosciences). Data analysis was performed using FlowJo v9.3 (TreeStar, Inc.) as shown in Figure S5 and the Mean Fluorescence Intensity (MFI) of the 90^th percentile was used for downstream analysis. All stimulated samples were compared to unstimulated control (baseline) samples and fold-changes calculated. Responses above 1.5-fold change as well as baseline MFI values were used for heritability estimates.

5. Serum protein quantification

Blood samples were centrifuged and stored at −80° C awaiting analysis. Human 51-plex were purchased from Affymetrix (Santa Clara, CA) and used according to the manufacturer’s recommendations with modifications as described below. Briefly, samples were mixed with antibody-linked polystyrene beads on 96-well filter-bottom plates and incubated at room temperature for 2 h followed by overnight incubation at 4°C. Room temperature incubation steps were performed on an orbital shaker at 500–600 rpm. Plates were vacuum filtered and washed twice with wash buffer, then incubated with biotinylated detection antibody for 2 h at room temperature. Samples were then filtered and washed twice as above and resuspended in streptavidin-PE. After incubation for 40 minutes at room temperature, two additional vacuum washes were performed, and the samples resuspended in Reading Buffer. Each sample was measured in duplicate. Plates were read using a Luminex 200 instrument with a lower bound of 100 beads per sample per protein. Each sample was measured in duplicate. Plates were read using a Luminex LabMap200 instrument with a lower bound of 100 beads per sample per protein per well. The Luminex LabMap200 outputs the fluorescence intensity of a given protein in a sample. For each well, we considered the median fluorescence intensity (MFI) for a serum protein as its abundance and averaged the MFI of these replicates. To ignore low abundance proteins, only measurements with mean concentrations higher than a negative control serum were included in our analysis.

6. Hemagglutination Inhibition (HAI) Assays

The HAI assay was performed on sera from day 0 and day 28-post influenza vaccination. Fold-changes day28/day0 was used for analyses. In one analysis, subjects with a pre-vaccination titer of 40 or more were excluded (Table 1). Serially diluted 25 microliter aliquots of serum samples in PBS were mixed with 25 microliter aliquots of virus, corresponding to 4 HA units, in V-bottom 96-well plates (Nunc, Rochester, NY, USA). These were then incubated for 15 min at room temperature. At the end of the incubation, 50 microliters of 0.5% chicken red blood cells (cRBC) were added and the plate incubated for one hour at room temperature and HAI activity was read as follows: I) Postive result: hemagglutination is present, the well is hazy with no cRBC button, or II) Negative result: hemagglutination is absent, the well is relatively clear with cRBC button. The HAI titer is defined as the reciprocal of dilution of last well that inhibits hemagglutination.

7. Structural equation modeling to estimate heritable and non-heritable influences

For all the measurements made, a structural equation modeling approach was used to estimate heritability (Rijsdijk and Sham, 2002). This classical twin modeling approach is based on the assumption that MZ twins are genetically identical while DZ twins share approximately 50% of their polymorphic genes, and that MZ and DZ twin pairs are equally similar with respect to their shared environmental influence. The covariance matrix for each measurement made in MZ and DZ twin pairs can then be decomposed into three parameters: I) Additive genetic parameter II) Common environmental parameter and III) Environmental parameter unique to one twin. A linear ACE model then estimates the contribution of each of these parameters by maximum likelihood. After correcting the E-parameter for technical measurement errors, as described in detail below, all parameters were scaled as a proportion of the total variance (A+C+E). All data was corrected for age and gender prior to ACE modeling. We performed resampling tests on all heritability estimates, using a jackknife resampling approach leaving one twin pair out of the calculation in each iteration, and using the mean values of all such iterations as our final estimate with 95% confidence intervals. All calculations were performed both using our own implementation of ACE-fitting in MATLAB, version 8.3, and using the openMX library version 1.4 running in R version 3.0.3. Both platforms were used in order to prevent any bias due to the software used.

8. Correction of model estimates for technical variability

For all measurements, standard samples were analyzed repeatedly (>17 times). These were either PBMC aliquots frozen at the same time and thawed for every experiment or pooled serum used as standards for Luminex and HAI assays. By calculating pooled variance estimates for these technical replicates and subtracting this from the E-component in our ACE models prior to normalization, we prevented the underestimation of heritability due to such stochastic measurement errors. This procedure overestimates the technical noise level of the actual twin samples by being collected across multiple batches, while twin samples compared were always analyzed within the same batch and after correction no relationship between technical variability and heritability estimates are seen (Figure S6A). We also assessed the biological variability over time in an unrelated cohort, by calculating coefficients of variance (CVs) from longitudinal samples drawn once yearly for 2–5 years. No relationship between longitudinal CVs and estimated heritability was found (Figure S6B).

9. Identification of pairwise dependencies between measurements and the creation of an immune network model

To identify meaningful relationships between measurements, we used the recently developed non-paranormal SKEPTIC approach (Liu et al., 2012), with a transformed Spearman/rank correlation matrix as input. To make the network model interpretable, we pursued sparse precision matrices with a graphical lasso approach penalizing non-zero entries in the matrix (Friedman et al., 2008). A zero entry in the precision matrix encodes conditional independencies of pairs of measurements given the state of all other measurements and is less sensitive to spurious indirect connections as compared to simple correlation analyses. To validate the inferred structural relationships, two tests were performed, 1) a permutation of samples on per phenotype basis was done to obtain null distributions for entries of the precision matrix. By repeatedly producing permuted samples and running our procedure on those samples, we obtain a distribution of precision matrices that is mostly dominated by very sparse, diagonal, matrices but also has occurrences of precision matrices that have non-zero entries. The non-zero entries thus obtained were false positive and hence we can estimate which entries obtained on the actual data exceed the size of these false positives. 2) Given relative infrequent occurrences of non-zero entries when fitting sparse precision matrices, we opted to estimate confidence intervals for each entry in the precision matrix. We deemed entries whose confidence interval contained 0 insignificant. To obtain these confidence intervals we performed bootstrap analysis, by resampling the real samples.

10. CMV serology

CMV serology was determined using a commercially available ELISA kit (CMV IgG, Gold Standard Diagnostics, Davis, CA) as per manufacturer’s instructions. Briefly, sera stored at −80°C were thawed to room temperature (20–25°C) and diluted 1:51 in kit diluent. Diluted samples were added to wells coated with CMV antigen from strain AD169 and incubated at room temperature for 30 minutes. Wells were washed and drained, followed by the addition of goat anti-human IgG antibodies labeled with calf alkaline phosphatase, and incubated at room temperature for 30 minutes. Wells were again washed and drained, followed by addition of p-nitrophenyl phosphate substrate, and incubated at room temperature for 30 minutes. After the addition of 0.5M trisodium phosphate stop solution, the absorbance of each well at 405nm was read and results analyzed using the manufacturer’s instructions.

The authors thank all members of the Davis and Y. Chien labs for inspiring discussions and the staff at the Human Immune Monitoring Center and the Stanford-LPCH Vaccine Program for sample collection, and data generation. Thanks also to Rob Tibshirani and Chiara Sabatti for input on statistical analyses. Twins were recruited from Twin Research Registry at SRI International, Menlo Park, CA and the authors wish to thank the following individuals: Mary McElroy, Lisa Jack, Ruth Krasnow, Jill Rubin, Dina Basin, Lucia Panini, and Marty Ritchey. The Stanford-LPCH program thanks Project Manager Sally Mackey, Research Nurses Susan Swope, Tony Trela and Cynthia Walsh, phlebotomist Michele Ugur, and CRAs Ashima Goel, Thu Quan, Kyrsten Spann, Sushil Batra, Isaac Chang and Raquel Fleischmann. The authors also wish to thank the contributions and commitment to science provided by the twins through their ongoing participation in the Registry and various research studies. Funds to the Twin Registry were provided through SRI’s Center for Health Sciences and through grants from the NIH including: DA011170, DA023063, AI057229, AI090019, and ES022153. This work was also supported by the Wenner-Gren Foundation and Sweden-America Foundation to PB, and NIH grants U19 AI090019 and U19 AI057229 to MMD as well as an NIH/NCRR CTSA award number UL1 RR025744 to the Clinical and Translational Research Unit, Stanford University Hospital.