Etienne Patin
email: epatin@pasteur.fr

Milena Hasan

Jacob Bergstedt

Vincent Rouilly

Valentina Libri

Alejandra Urrutia

Cécile Alanio

Petar Scepanovic

Christian Hammer

Friederike Jönsson

Benoît Beitz

Hélène Quach

Yoong Wearn Lim

Julie Hunkapiller

Magge Zepeda

Cherie Green

Barbara Piasecka

Claire Leloup

Lars Rogge

François Huetz

Isabelle Peguillet

Olivier Lantz

Magnus Fontes

James P Di

Stéphanie Thomas

Jacques Fellay

Darragh Duffy

Lluís Quintana-Murci

Matthew L Albert
email: albertm7@gene.com

Natural variation in innate immune cell parameters is preferentially driven by genetic factors

The enumeration and characterization of circulating immune cells provide key indicators of human health and disease. To identify the relative impact that environmental and genetic factors have on variation of innate and adaptive immune cell parameters in homeostatic conditions, we combined standardized flow cytometric analysis of blood leukocytes and genome-wide DNA genotyping in 1,000 healthy, unrelated individuals of western European ancestry. We show that smoking, together with age, sex and latent cytomegalovirus infection, are the main non-genetic factors affecting human variation in immune cell parameters. Genome-wide association studies of 166 immunophenotypes identified 15 loci that are enriched in disease-associated variants. Finally, we demonstrate that innate cell parameters are more strongly controlled by genetic variation than adaptive cell parameters, which are primarily driven by environmental exposures. Our data establish a resource that generates new hypotheses in immunology and highlight the role of innate immunity in the susceptibility to common autoimmune diseases.

Introduction

The immune system plays an essential role in maintaining homeostasis in individuals challenged by microbial infections, a physiological mechanism conceptualized by the French physician Claude Bernard in 1865, when he defined the notion of "milieu intérieur" [START_REF] Bernard | Introduction à l'étude de la médecine expérimentale[END_REF] . Host-pathogen interactions trigger immune responses through the activation of specialized immune cell populations, which may eventually result in pathogen clearance. The study of immune cell populations circulating in the blood provides a view into innate cells that are transiting between the bone marrow and tissues, and adaptive cells that are recirculating through lymphoid organs. Clinical studies of patients with past or chronic latent infections have reported profound perturbations of subsets of circulating immune cells due to altered trafficking, selective expansion or attrition [START_REF] Altfeld | Innate immunity against HIV-1 infection[END_REF][START_REF] Orme | The balance between protective and pathogenic immune responses in the TB-infected lung[END_REF] . However, several studies have suggested that extensive differences in white blood cell composition also exist among healthy individuals [START_REF] Tollerud | The Influence of Age, Race, and Gender on Peripheral Blood Mononuclear-Cell Subsets in Healthy Nonsmokers[END_REF][START_REF] Reichert | Lymphocyte Subset Reference Ranges in Adult Caucasians[END_REF] . The evaluation of the naturally occurring variation of immune cell parameters, together with its environmental and genetic determinants, could accelerate hypothesis generation in basic immunology, and ultimately improve the characterization of pathological states.

Population immunology approaches, which compare the immune status across a large number of healthy individuals, have highlighted the predominant effect of intrinsic factors such as age and sex on human blood cell composition [START_REF] Liston | Shaping Variation in the Human Immune System[END_REF] . Several activated and memory T cell subpopulations increase with age [START_REF] Goronzy | Successful and maladaptive T cell aging[END_REF] , which may partially result from diminished thymic activity [START_REF] Sauce | Altered thymic activity in early life : how does it affect the immune system in young adults ?[END_REF] and explain reduced vaccination efficacy in the elderly [START_REF] Furman | Apoptosis and other immune biomarkers predict influenza vaccine responsiveness[END_REF] . Seasonal fluctuations in B cells, regulatory T (T reg ) cells and monocytes [START_REF] Aguirre-Gamboa | Differential Effects of Environmental and Genetic Factors on T and B Cell Immune Traits[END_REF] and a strong effect of cohabitation on human immune profiles [START_REF] Carr | The cellular composition of the human immune system is shaped by age and cohabitation[END_REF] have been observed, suggesting that environmental exposures also drive immune variation. For instance, latent cytomegalovirus (CMV) infection, detected in 40% to >90% of the general population [START_REF] Boeckh | Cytomegalovirus: pathogen, paradigm, and puzzle[END_REF] , has been associated with an increased number of effector memory T cells [START_REF] Wertheimer | Aging and Cytomegalovirus Infection Differentially and Jointly Affect Distinct Circulating T Cell Subsets in Humans[END_REF] , which could in turn alter immune responses to heterologous infection [START_REF] Furman | Cytomegalovirus infection enhances the immune response to influenza[END_REF] . However, the respective impact of age, sex and CMV infection on both innate and adaptive cells, as well as the precise nature of the environmental factors affecting immune variation, are largely unknown.

Recent technological advances in flow cytometry, combined with genome-wide DNA genotyping, now allow the dissection of the genetic basis of interindividual variation in immune cell parameters. A seminal genome-wide association study identified 13 genetic loci strongly associated with the proportion of different leukocyte subpopulations, in a cohort of 249 Sardinian families [START_REF] Orrù | Genetic variants regulating immune cell levels in health and disease[END_REF] . Another study reported the deep immunophenotyping of ~1,800 independent traits in 245 healthy twin pairs, identifying 11 independent genetic loci that accounted for up to 36% of the variation of 19 different traits [START_REF] Roederer | The Genetic Architecture of the Human Immune System: A Bioresource for Autoimmunity and Disease Pathogenesis[END_REF] . A third study estimated the genetic heritability of 95 different immune cell frequencies in 105 healthy twin pairs, and suggested that variation in immune cells is largely explained by non-heritable factors [START_REF] Brodin | Variation in the Human Immune System Is Largely Driven by Non-Heritable Influences[END_REF] . Finally, four novel loci were associated to B and T cell traits in a cohort of 442 healthy donors, in a study that dissected both non-genetic and genetic factors affecting immune cell traits mediating adaptive immunity [START_REF] Aguirre-Gamboa | Differential Effects of Environmental and Genetic Factors on T and B Cell Immune Traits[END_REF] . Together, these studies provided valuable insights into the contribution of genetic factors to inter-individual differences in adaptive immune cell populations, but largely neglected several major innate cell types in circulation. An integrated evaluation of the nature and respective impact of intrinsic, environmental and genetic factors driving human variation in both innate and adaptive immunity is thus lacking.

Here, we report the use of standardized flow cytometry to comprehensively establish the white blood cell composition of 1,000 healthy, unrelated individuals of western European ancestry, which compose the Milieu Intérieur cohort. We confirm with this broad resource that age, sex, CMV seropositivity and smoking have major, independent effects on innate and adaptive immune cell parameters. We identified, through a genome-wide association study, 15 loci associated with parameters of circulating leukocyte subpopulations, 12 of which are novel. Finally, we show that cellular mediators of innate and adaptive immunity are differentially affected by non-genetic and genetic factors under homeostatic conditions.

Results

Variation of immune cell parameters in the general population

The Milieu Intérieur cohort includes 500 men and 500 women, stratified across five decades of age from 20 to 69 years. Subjects were surveyed for a number of demographic variables, including past infections, vaccination and surgical histories and health-related habits (Supplementary Table 1).

Detailed inclusion and exclusion criteria used to define "healthy" subjects recruited into the cohort have been previously reported [START_REF] Thomas | The Milieu Intérieur study -An integrative approach for study of human immunological variance[END_REF] .

To describe natural variation of both innate and adaptive immune cells in the 1,000 subjects, we used ten 8-color immunophenotyping flow cytometry panels (Supplementary Figs. 1-10 and Supplementary Table 2; Online Methods), which allowed us to report a total of 166 distinct immunophenotypes (Supplementary Table 3). Our resource includes 75 (46%) and 91 (54%) immunophenotypes obtained in innate and adaptive immune cells, respectively. Innate cells were defined as those lacking somatic recombination of the genome [START_REF] Vivier | Innate or Adaptive Immunity? The Example of Natural Killer Cells[END_REF] , and included granulocytes (neutrophils, basophils and eosinophils), monocytes, natural killer (NK) cells, dendritic cells and innate lymphoid cells (ILCs) (Fig. 1). Adaptive cells were defined by their dependence on RAG1/2 activity and included T cells (γδ T, MAIT, NKT, T reg and T H cells) and B cells. The immunophenotypes in both innate and adaptive immune cells included 76 absolute counts of circulating cells, 87 expression levels of cell-surface protein markers (quantified by the mean fluorescence intensity, or MFI), and 3 ratios of cell counts or MFI (Supplementary Fig. 11 and

Supplementary Table 3).

To reduce technical variation introduced by sample temperature fluctuations and pre-analytical procedures, we strictly followed a standardized protocol for tracking and processing samples [START_REF] Hasan | Semi-automated and standardized cytometric procedures for multi-panel and multiparametric whole blood immunophenotyping[END_REF] . We verified that measured immunophenotypes were highly reproducible using technical replicates (Supplementary Figs. S12 andS13 and Supplementary Table 3), demonstrating the high precision of the data. We nevertheless identified two technical batch effects that impacted flow cytometric analyses. One effect corresponded to the hour at which the blood sample was drawn from fasting subjects (Supplementary Fig. 14a), which may possibly be explained by the spike in cortisol at the time of waking [START_REF] Patterson | Cortisol Patterns Are Associated with T Cell Activation in HIV[END_REF] . The second effect corresponded to temporal variation of immunophenotypes over the one-year sampling period, which did not follow the periodic distribution observed for cellular traits under seasonal fluctuations [START_REF] Carr | The cellular composition of the human immune system is shaped by age and cohabitation[END_REF] , and primarily affected MFI measures (Supplementary Fig. 14b).

We corrected for these batch effects in all subsequent analyses (Supplementary Fig. 15; Online Methods), and provide the distribution, ranges and statistics of all batch-corrected immune cell counts (Supplementary Table 3), thereby facilitating comparisons with cytometry data collected as part of routine clinical practice. This resource can be accessed through a user-friendly web application (http://104.236.137.56:3838/LabExMICytometryBrowser_ShinyApp/), which can be queried based on personal characteristics, such as age or sex.

Owing to the hierarchical structure of immune cell differentiation (i.e., cellular lineages emerge from common progenitor cells), a substantial portion of the immune cell counts measured in this study were highly correlated (Supplementary Fig. 16). These correlations were not directly attributable to the influence of factors such as age or sex, which were regressed out in this analysis. We observed correlations between circulating levels of ILC and NK populations, reflecting their common developmental pathway and dependence on γ c cytokines [START_REF] Serafini | Transcriptional regulation of innate lymphoid cell fate[END_REF] . Likewise, MAIT cells and CCR6 + CD8 + T cells were also correlated, owing to the former being the major subset of CCR6 + T cells in circulation [START_REF] Dusseaux | Human MAIT cells are xenobiotic-resistant, tissue-targeted, CD161hi IL-17secreting T cells[END_REF] . Finally, we identified a strong correlation between the number of T reg and conventional CD4 + T cells, validating previous experimental work that defined an IL-2-driven self-regulatory circuit that integrates the homeostasis of these cell populations [START_REF] Amado | IL-2 coordinates IL-2-producing and regulatory T cell interplay[END_REF] .

Impact of age, sex and CMV infection on innate and adaptive cell parameters

Prior studies have shown that two intrinsic factors, age and sex, are responsible for inter-individual variation in white blood cell composition [START_REF] Liston | Shaping Variation in the Human Immune System[END_REF][START_REF] Goronzy | Successful and maladaptive T cell aging[END_REF][START_REF] Aguirre-Gamboa | Differential Effects of Environmental and Genetic Factors on T and B Cell Immune Traits[END_REF][START_REF] Furman | Cytomegalovirus infection enhances the immune response to influenza[END_REF][START_REF] Pennell | Sex affects immunity[END_REF][START_REF] Furman | Systems analysis of sex differences reveals an immunosuppressive role for testosterone in the response to influenza vaccination[END_REF][START_REF] Astle | The Allelic Landscape of Human Blood Cell Trait Variation and Links to Common Complex Disease[END_REF] . We used linear mixed models to quantify the respective impact of each of these intrinsic factors on variation in innate and adaptive cell composition. We observed a significant effect of age on 35% of immune cell parameters (adjusted P<0.01; Fig. 2a and Supplementary Fig. 17a), among which only 29% were measured in innate cells. We detected a general decline in the number of ILC and plasmacytoid dendritic cells (pDCs) and an increase in the number of CD16 hi monocytes with increasing age (Fig. 2a), which might contribute to the altered immune response to viral infections in elderly persons and age-associated inflammation [START_REF] Furman | Cytomegalovirus infection enhances the immune response to influenza[END_REF][START_REF] Della Bella | Peripheral blood dendritic cells and monocytes are differently regulated in the elderly[END_REF][START_REF] Puchta | TNF Drives Monocyte Dysfunction with Age and Results in Impaired Antipneumococcal Immunity[END_REF] . We found a modest increase in the number of memory T cells with age, supporting the view that the observed expansion of these cell populations in elderly subjects is not due to aging per se, but to CMV seropositivity [START_REF] Wertheimer | Aging and Cytomegalovirus Infection Differentially and Jointly Affect Distinct Circulating T Cell Subsets in Humans[END_REF] , which we accounted for in the model. Our analyses also revealed that naive CD8 + T cells decrease more than twice as rapidly with age as compared to naive CD4 + T cells, at a rate of 3.6 % (99% FCR-adjusted Confidence Interval (99%CI): [3.0%, 4.1%]) and 1.6 % (99%CI: [1.1%, 2.1%]) per year, respectively (Fig. 2a-c), supporting the view that CD8 + T cells are more susceptible to concentrations of homeostatic cytokines and/or that the production of CD4 + T cells is preferentially enhanced in the human thymus [START_REF] Vrisekoop | Sparse production but preferential incorporation of recently produced naïve T cells in the human peripheral pool[END_REF] .

Although sex differences have been previously reported for various immune responses and diseases [START_REF] Pennell | Sex affects immunity[END_REF] , studies examining circulating cellular parameters have reported inconsistent results, owing to both differences between flow cytometry procedures and relatively small, underpowered or poorlystratified study cohorts. We report a significant impact of sex on 16% of measured immunophenotypes (adjusted P<0.01, Fig. 2d and Supplementary Fig. 17b), of which 38% were measured in innate cells. We found a higher number of activated NK cells in men, as compared to women. By contrast, MAIT cells were systematically increased in women, across all age decades (Fig. 2e-f), collectively suggesting a lasting effect of early hormonal differences on immune cell development and biology.

Environmental exposures are also known to drive immune variation, among which persistent CMV infection is one of the strongest candidates [START_REF] Liston | Shaping Variation in the Human Immune System[END_REF][START_REF] Wertheimer | Aging and Cytomegalovirus Infection Differentially and Jointly Affect Distinct Circulating T Cell Subsets in Humans[END_REF][START_REF] Furman | Cytomegalovirus infection enhances the immune response to influenza[END_REF][START_REF] Brodin | Variation in the Human Immune System Is Largely Driven by Non-Heritable Influences[END_REF] . We observed a significant effect of latent CMV infection on 13% of immune cell parameters (Fig. 2g and Supplementary Fig. 17c), of which more than 75% were measured in adaptive cells. We confirm that CMV triggers a major change in the number of memory T cells, which is independent from age effects [START_REF] Wertheimer | Aging and Cytomegalovirus Infection Differentially and Jointly Affect Distinct Circulating T Cell Subsets in Humans[END_REF][START_REF] Brodin | Variation in the Human Immune System Is Largely Driven by Non-Heritable Influences[END_REF] . In particular, CMV seropositivity associated with a 12.5-fold (99%CI: [8.8, 17.6]) higher number of CD4 + effector memory RA T cells (T EMRA ), and a 4.6-fold (99%CI: [3.5, 6.0]) higher number of CD8 + T EMRA cells (Fig. 2g-i). However, we did not find evidence that CMV infection impacts the number of naive or central memory (T CM ) T cell compartments. Supporting this observation, the total number of CD8 + and CD4 + T cells increased in parallel with the expanded number of memory T cells, thus suggesting independent regulation of the naive and T EM and/or T EMRA cell pools. CMV seropositive donors also presented lower numbers of circulating NKT and MAIT cells (Fig. 2g). Together, our broad resource provides a comprehensive quantification of the respective impact that age, sex and CMV infection have on immune cell parameters. In doing so, our results suggest a stronger impact of these factors on adaptive cells, relative to innate cells.

Tobacco smoking extensively alters innate and adaptive cell numbers

Capitalizing on the detailed lifestyle and demographic data obtained for the Milieu Intérieur cohort, we evaluated the influence of additional environmental factors on immune cell parameters, controlling for the defined effects of age, sex and CMV serological status. A total of 39 variables were chosen for analysis and tested for association with each immunophenotype. These include socio-economic characteristics, past infections, health-related habits and surgery and vaccination history (Supplementary Fig. 18 and Supplementary Table 1). We identified a unique environmental factor that significantly alters circulating numbers of immune cells: active tobacco cigarette smoking, which affects 36% of measured immunophenotypes (Fig. 3a and Supplementary Fig. 19), of which 36% were measured in innate cells.

We observed a 23% (99%CI: [11%, 37%]) increase in the number of circulating CD45 + cells, and a 26% (99%CI: [10%, 45%]) increase in the number of conventional lymphocytes in smokers as compared to non-smokers (Fig. 3b). Previous studies suggested that smokers have alterations in circulating cell populations due to diminished adherence of leukocytes to blood vessel walls, possibly as a result of lower antioxidant concentrations [START_REF] Tsuchiya | Smoking a Single Cigarette Rapidly Reduces Combined Concentrations of Nitrate and Nitrite and Concentrations of Antioxidants in Plasma[END_REF] . Furthermore, we found in active smokers a significant increase of 43% (99%CI: [17%, 76%]) and 41% (99%CI: [15%, 71%]) of activated and memory T reg cells, respectively, a pattern that was also observed to a lesser extent in past smokers (Fig. 3b-d). Active smokers also showed decreased numbers of NK cells, ILCs, γδ T cells and different subsets of MAIT cells (Fig. 3b). These findings are consistent with a study showing that smoking triggers local release of IL-33 by the lung epithelium [START_REF] Kearley | Cigarette Smoke Silences Innate Lymphoid Cell Function and Facilitates an Exacerbated Type I Interleukin-33-Dependent Response to Infection[END_REF] , in turn engaging the IL-33 receptor, ST2, on both innate and non-classical lymphocytes [START_REF] Monticelli | Articles Innate lymphoid cells promote lung-tissue homeostasis after infection with influenza virus[END_REF] . Collectively, these findings reveal that active smoking has a profound impact on immune cell parameters, which is similar in magnitude to that of age, and affects both innate and adaptive cells.

Genome-wide association study of 166 immune cell parameters

To identify common genetic variants affecting inter-individual variation in immune cell parameters, the Milieu Intérieur cohort was genotyped at 945,213 SNPs, enriched in exonic SNPs (Online Methods). After quality control (Supplementary Fig. 20), genotype imputation was performed and yielded a total of 5,699,237 highly accurate SNPs, which were tested for association with the 166 immunophenotypes using linear mixed models. The models were adjusted for the genetic relatedness among subjects and any non-genetic variable identified as predictive of each specific immunophenotype by stability selection based on elastic net regression (Supplementary Table 3; Online Methods). We confirmed our power to identify medium-effect genotype-phenotype associations by simulations, and by empirically replicating well-known genetic associations with nonimmune traits, such as eye and hair color or uric acid and cholesterol levels (Online Methods).

With respect to immune traits, we found 14 independent genetic loci associated with 42 out of 166 immunophenotypes (25%), at a conservative genome-wide significant threshold of P<1.0x10 -10 (Fig. 4a, Table 1, Supplementary Fig. 21, Supplementary Tables 4 and5). We then conducted conditional GWAS, by adjusting these 42 immunophenotypes on the 14 leading associated variants (Table 1), and found an additional independent locus reaching genome-wide significance (Supplementary Fig. 22 and Supplementary Table 6). Genome-wide significant associations were replicated in an independent cohort of 75 European-descent donors, for all immune traits measured in this replication cohort (P<0.05; Table 1; Online Methods). Also, we confirmed that our immune cell measurements were stable, as all genome-wide significant associations were confirmed for immunophenotypes measured in a new blood draw taken in 500 of the 1,000 subjects of the Milieu Intérieur cohort, sampled 7 to 44 days after the initial visit (P<10 -3 ; Table 1). We also provide a list of 26 suggestive association signals (P<5.0x10 -8 ), including a number of biologically relevant candidate genes (Supplementary Table 6). The associated genetic loci were enriched in SNPs associated by GWAS with diseases (31% observed vs. 5% expected, resampling P=0.0032), most of which were autoimmune diseases, including rheumatoid arthritis, Vogt-Koyanagi-Harada syndrome and atopic dermatitis (Supplementary Table 4). These findings highlight the importance of loci altering immune cell populations in the context of ultimate organismal traits affecting human health.

Genetic associations primarily identify immune cell-specific protein QTLs

Of the 42 immunophenotypes for which a significant genetic association was detected, 36 (86%) were MFI, which measures the cell-specific expression of protein markers conventionally used to determine the differentiation or activation state of leukocytes. For 28 of these 36 MFI measurements (78%), the genetic association was observed between the protein MFI and SNPs located in the vicinity of the gene encoding the corresponding protein (Table 1 and Supplementary Fig. 21), i.e., local protein QTLs (local-pQTLs). For instance, genetic variation close to the ENPP3 gene was associated with CD203c MFI in basophils (rs2270089, P=2.1x10 -28 ), CD24 with CD24 MFI in marginal zone B cells (rs12529793, P=3.8x10 -21 ) and CD8A with CD8a MFI in CD69 + CD16 hi NK cells (rs71411868, P=5.9x10 -58 ).

We identified two independent local-pQTLs in the FCGR gene cluster (Table 1), which encodes the most important Fc receptors for inducing phagocytosis of opsonized microbes. Genetic variation close to FCGR3A was associated here with CD16 MFI in CD16 hi NK cells (rs3845548, P=3.0x10 -87 ).

The same variants were also shown to affect the number of CD62L-myeloid cDCs in a previous study [START_REF] Orrù | Genetic variants regulating immune cell levels in health and disease[END_REF] . The second signal associated FCGR2B variation with CD32 MFI in basophils (rs61804205, P=1.7x10 -36 ), but not in eosinophils and neutrophils. Consistently, it is known that basophils express both CD32a and CD32b proteins, while eosinophils and neutrophils predominantly express CD32a [START_REF] Cassard | Fcγ Receptors Inhibit Mouse and Human Basophil Activation[END_REF] .

Conversely, a local-pQTL was identified at the SELL gene, which was associated with CD62L MFI in eosinophils and neutrophils (rs2223286, P=1.6x10 -35 and 8.8x10 -13 , respectively), but not in basophils (Fig. 4b,c).

A number of other local-pQTLs were found to be cell-specific; three different association signals were found in the HLA-DR gene region, with HLA-DR MFI in pDCs and CD14 hi monocytes (rs114973966, P=2.2x10 -56 ), in cDC1 (rs2760994, P=6.1x10 -38 ) and in cDC3 cells (rs143655145, P=2.6x10 -11 ). To verify if these signals were independent from each other, we conducted omnibus association tests on imputed HLA alleles (Online Methods). We found that the association signals in CD14 hi monocytes, pDCs and cDC1 actually resulted from different amino acid-altering variants at the same multi-allelic position 13 of the HLA-DRβ1 protein (P=2.0x10 -47 , 7.0x10 -90 and 5.3x10 -41 in CD14 hi monocytes, pDC and cDC1, respectively; Supplementary Tables 7 and8), recently shown to explain a large part of the association signal in the HLA locus for type 1 diabetes [START_REF] Hu | Additive and interaction effects at three amino acid positions in HLA-DQ and HLA-DR molecules drive type 1 diabetes risk[END_REF] . A different aminoacid variant, at position 67 of HLA-DRβ1, was identified in cDC3s (P=3.9x10 -13 ). Conditional analyses also revealed independent associations of HLA-DR cell-surface expression with two residues in class I HLA-B gene (position 97 and 194; P=3.8x10 -17 and 1.3x10 -18 ; Supplementary Tables 7 and8). Collectively, these results show that the protein expression of markers of immune cell differentiation and activation can be affected by common genetic variants, of which some are known to be implicated in human pathogenesis.

Immune cell local protein QTLs control mRNA levels of nearby genes

Although four of the 9 local-pQTLs identified by our analyses are likely explained by amino acidaltering variants in surrounding genes (Supplementary Tables 4 and7), the remaining signals do not present obvious candidate causal variants. To dissect the functional basis of these associations, we tested if the corresponding SNPs were also associated with mRNA levels of nearby genes (i.e., expression QTL, eQTL) using gene expression data obtained from the same donors [START_REF] Piasecka | Distinctive Roles of Age, Sex and Genetics in Shaping Transcriptional Variation of Human Immune Responses to Microbial Challenges[END_REF] and results from the Genotype-Tissue Expression (GTEx) Project [START_REF]The Genotype-Tissue Expression (GTEx) pilot analysis: Multitissue gene regulation in humans[END_REF] . Five of the local-pQTLs were strongly associated with the transcript levels of a surrounding gene (P<1.0x10 -5 ; Fig. 4d). The SNPs controlling the MFI of CD16 in CD16 hi NK cells and CD32 in basophils, CD62L in eosinophils, CD8a in CD69 + CD16 hi NK cells and CD203c in basophils were associated with mRNA levels of FCGR2B, SELL, CD8A, and ENPP3, respectively (Supplementary Table 4). These analyses indicate that genetic variants associated with immunophenotypes can directly affect gene expression of markers of immune cells in whole blood. This suggests that eQTL mapping in different immune cell compartments can greatly improve our knowledge of the genetic factors controlling human inter-individual variation in flow cytometric parameters.

Novel trans-acting genetic associations with immune cell parameters

We detected six loci that do not exclusively act as local-pQTLs on immunophenotypes (Table 1 and Supplementary Fig. 21). These included variants that are associated with immune cell counts, or that are genetically independent from the genes encoding immune cell markers with which they are associated (i.e., trans-pQTLs). A variant in the vicinity of the S1PR1 gene was associated with CD69 MFI in CD16 hi NK cells (rs6693121, P=4.8x10 -37 ). CD69 is known to downregulate cell-surface expression of the sphingosine-1-phosphate receptor-1 (S1P1) on lymphocytes, a mechanism that elicits egress from the thymus and secondary lymphoid organs [START_REF] Garris | Sphingosine-1-phosphate receptor 1 signalling in T cells : trafficking and beyond[END_REF] . Genetic variation in an intron of the ACOXL gene, close to BCL2L11, was associated with the absolute count of CD8a + CD56 hi NK cells (rs12986962, P=9.1x10 -19 ). BCL2L11 (also known as BIM) is an important regulator of lymphocyte apoptosis [START_REF] Pellegrini | Loss of Bim Increases T Cell Production and Function in Interleukin 7 Receptordeficient Mice[END_REF] , and is associated with chronic lymphocytic leukemia and total blood cell number [START_REF] Van Der Harst | Seventy-five genetic loci influencing the human red blood cell[END_REF] . A third association involved genetic variants close to the ACTL9 gene and the ratio of CD16 MFI in CD16 hi and CD56 hi NK cells (rs114412914, P=4.3x10 -30 ). The same variants have been also found to be associated with CD56 ++ CD16 -NK cells in another study [START_REF] Aguirre-Gamboa | Differential Effects of Environmental and Genetic Factors on T and B Cell Immune Traits[END_REF] .

Although identified here for their trans effects on markers of immune cell differentiation or activation, three trans-acting genetic associations were also local-eQTLs for nearby immune-related genes [START_REF]The Genotype-Tissue Expression (GTEx) pilot analysis: Multitissue gene regulation in humans[END_REF] (Supplementary Tables 4 and6). The MFI of CCR7 in CD4 + and CD8b + naive T cells was associated with a variant in the TMEM8A gene (rs11648403, P=3.0x10 -19 ), which also controls TMEM8A mRNA levels (P=2.5x10 -27 ). TMEM8A is expressed on the surface of resting T cells and is down-regulated after cell activation [START_REF] Motohashi | Molecular Cloning and Chromosomal Mapping of a Novel Protein Gene, M83[END_REF] , suggesting a possible functional association and/or coregulation with CCR7. Variants in the vicinity of the ALOX15 gene were associated with increased protein levels of the high-affinity IgE receptor in eosinophils (rs56170457, P=9.2x10 -14 ) and increased ALOX15 mRNA levels (P=2.7x10 -13 ). These results, together with the high expression of the ALOX15 protein and its pro-inflammatory effect in circulating eosinophils [START_REF] Feltenmark | Eoxins are proinflammatory arachidonic acid metabolites produced via the 15lipoxygenase-1 pathway in human eosinophils and mast cells[END_REF] , suggest that this lipoxygenase plays an important role in IgE-dependent allergic reactions. Finally, conditional GWAS identified an additional trans-acting association, between a variant close to the CD83 gene and HLA-DR MFI in cDC1 (rs72836542, P=2.8x10 -12 , Supplementary Fig. 22), the same variant being also identified as a local-eQTL of CD83 gene expression (P=5.4x10 -21 ). These results suggest that CD83, an early activation marker of human DCs, upregulates HLA-DR expression in activated dendritic cells.

Natural variation of innate immune cell parameters is preferentially driven by genetic factors

A large proportion of both MFI and cell number immunophenotypes that presented a genome-wide association were detected in innate immune cells (35/44, 80%), including granulocytes, monocytes, NK and dendritic cells (Table 1), while 47% of all immunophenotypes were measured in innate cells (Supplementary Table 3). Furthermore, of the adaptive cell immunophenotypes showing genetic associations, 3 of the 9 measurements (33%) were related to naive T or B cells, while naive adaptive cell parameters represented <10% of all adaptive cell measurements. These observations suggest a stronger effect of genetic variants on innate and naive adaptive cell subpopulations, relative to differentiated or experienced adaptive immune cells.

In support of this hypothesis, the presence of HLA-DR molecules, which was assessed at the surface of both innate and adaptive immune cells, was strongly associated with HLA-DR genetic variation in monocytes, NK and dendritic cells (Table 1), but not in memory CD4 + or CD8 + T CM , T EM and T EMRA cells (P>1.0x10 -6 ; Supplementary Table 5). Because we observed substantial correlations among HLA-DR + memory T cell numbers (R²≈0.3, P<0.05; Supplementary Fig. 16), we hypothesized that they were at least partly controlled by the same genetic factors, which were further examined using a multivariate GWAS (Online Methods). This refined approach detected a suggestive genetic association close the HLA-DRB1 gene with a variant (rs35743245, P=1.0x10 -8 ) in strong linkage disequilibrium with that detected in pDCs, monocytes and NK cells (r²=0.92; Supplementary Fig. 23). This finding provides proof-of-concept that immunophenotypes in both innate and adaptive cells can be controlled by the same genetic factors, but their effects are stronger in innate cells, relative to experienced adaptive cells.

We next systematically quantified the impact of genetic and non-genetic factors on innate and adaptive cells. We established, for each immunophenotype, a linear regression model that included the four most impactful non-genetic variables (Figs. 2 and3) and all genome-wide significant and suggestive variants (Table 1 and Supplementary Table 6), and estimated their respective contribution to the total variance (Online Methods). We found that a larger proportion of the variance of innate cell immunophenotypes was explained by genetic factors (Fig. 5b and5d), relative to adaptive cell immunophenotypes (Fig. 5a and5c). Inversely, the variance in adaptive cell numbers was dominated by non-genetic factors such as age and CMV serostatus (Fig. 5a). To test if these differences were significant, we used a mixed model that accounted for correlations among immunophenotypes (Online Methods). Conclusively, we estimated that the variance explained by genetics was 66% larger for innate cell measurements, relative to adaptive cells (95%CI: [13%-143%]; bootstrap P=0.012; Mann-Whitney U test: P=0.032), while the variance explained by nongenetic factors was 46% smaller for innate cell measurements (95%CI: [22%-63%]; bootstrap P=1.8x10 -3 ; Mann-Whitney U test: P=8.1x10 -3 ). When considering non-genetic factors separately, the ratio of explained variance between innate and adaptive cell measurements was the smallest for smoking (0.46, 95%CI: [0.17

Discussion

Over the last two decades, research in human immunology has employed multi-parametric cytometry to enumerate and assess the activation state of immune cells in healthy and disease conditions.

Although immune cell parameters do vary in the general population, the extent to which intrinsic, environmental and genetic factors explain this variability remained elusive. To tackle these questions, we generated a broad resource by combining standardized flow cytometry with genome-wide DNA genotyping in a demographically well-defined cohort of 1,000 healthy individuals. We confirm the strong and independent impacts of age and CMV infection on naive and memory T cell populations, respectively, and provide robust evidence for sex differences in innate and adaptive cell numbers. We show that immune homeostasis is altered upon chronic cigarette smoke exposure, which elicits both a decline of MAIT cells, possibly due to their increased migration to sites of inflammation, and an increase in the numbers of activated and memory T reg cells, suggesting a role for these immunosuppressive populations in the increased susceptibility of smokers to infection [START_REF] Stämpfli | How cigarette smoke skews immune responses to promote infection, lung disease and cancer[END_REF] . Furthermore, we found that human genetic variation substantially impacts immune cell parameters, particularly the cell-surface expression of markers conventionally used to identify leukocyte differentiation or activation. These results highlight the need to consider non-genetic and genetic features when interpreting parameters such as circulating white blood cells of patients, a critical aspect in clinical monitoring. For instance, HLA-DR expression on monocytes is routinely measured by flow cytometry to predict the clinical course of septic shock and identify patients who should benefit from immunoadjuvant therapies [START_REF] Venet | Monitoring the immune response in sepsis: a rational approach to administration of immunoadjuvant therapies[END_REF] . We identified a strong effect of HLA-DRβ1 coding variation on HLA-DR expression in CD14 hi monocytes, suggesting that prognostic tools of fatal outcome in sepsis should be tailored to patient's genetic makeup.

The most prominent result of our study is the lower number of genetic associations detected in memory T and B cells, relative to innate cells, an observation that could be explained by their strong dependence on the varying individual history of past infections. Adaptive immune cells are known to possess a much longer half-life as compared to myeloid innate cells, in mice and humans [START_REF] Kolaczkowska | Neutrophil recruitment and function in health and inflammation[END_REF][START_REF] Farber | Human memory T cells: Generation, compartmentalization and homeostasis[END_REF] .

Stimulus-induced differentiation and expansion may also result in the possible masking of genetic associations for adaptive cell types. Consistently, genetic associations in adaptive immune cells were primarily observed for immunophenotypes of naive adaptive cells. Our observations are further supported by a GWAS of 36 blood traits in 173,480 individuals, which found that the genetic heritability of monocyte and eosinophil counts was larger than that of lymphocyte counts [START_REF] Astle | The Allelic Landscape of Human Blood Cell Trait Variation and Links to Common Complex Disease[END_REF] . This is however at odds with another recent study, which concluded that adaptive immune traits are more affected by genetics, whereas innate immune traits are more affected by environment, based on the estimated genetic heritability of 23,394 immune phenotypes in 497 adult female twins [START_REF] Mangino | Innate and adaptive immune traits are differentially affected by genetic and environmental factors[END_REF] . We suggest that such deep immunophenotyping in large-scale cohorts, combined with statistical tests for differences in heritability that account for inherent correlations among phenotypes, may reveal a more balanced contribution of genetics on the natural variation of innate and adaptive immune cell traits.

Our findings that genetic factors preferentially controls variation in innate immune cells have other important consequences. A previous study of 105 healthy twin pairs concluded that variation in cell population frequencies is largely driven by non-heritable influences [START_REF] Brodin | Variation in the Human Immune System Is Largely Driven by Non-Heritable Influences[END_REF] . We find instead that genetic variation explains a large part of the variance of immune cell parameters, particularly MFIs (i.e., cell-surface expression of protein markers) measured in innate cells. This discrepancy may stem from the fact that this previous study considered only a fraction of innate myeloid and lymphoid populations [START_REF] Casanova | Disentangling Inborn and Acquired Immunity in Human Twins[END_REF] , and possibly because of its limited power due to a moderate sample size. Also, our results suggest that the genetic control of cell-surface expression of immune cell markers is stronger than that of cell counts, and the former were not assessed in most previous population immunology studies [START_REF] Aguirre-Gamboa | Differential Effects of Environmental and Genetic Factors on T and B Cell Immune Traits[END_REF][START_REF] Orrù | Genetic variants regulating immune cell levels in health and disease[END_REF][START_REF] Brodin | Variation in the Human Immune System Is Largely Driven by Non-Heritable Influences[END_REF] .

Finally, the mapping of genetic loci that control immune cell parameters identified cell-specific pQTLs that are enriched in genetic variants associated with human diseases and traits. For example, we identified the position 13 of the HLA-DRβ1 protein as a predictor of HLA-DR expression at the cell-surface of pDCs and monocytes, which in turn is strongly associated with type 1 diabetes [START_REF] Hu | Additive and interaction effects at three amino acid positions in HLA-DQ and HLA-DR molecules drive type 1 diabetes risk[END_REF] , suggesting the implication of innate immunity in the disease [START_REF] Astle | The Allelic Landscape of Human Blood Cell Trait Variation and Links to Common Complex Disease[END_REF] . Furthermore, the expression of CD56 and CD16 in NK cells is controlled by genetic variants close to the ACTL9 gene, which were shown to be associated with atopic dermatitis [START_REF] Paternoster | Multi-ancestry genome-wide association study of 21,000 cases and 95,000 controls identifies new risk loci for atopic dermatitis[END_REF] , suggesting a possible involvement of NK cells in this pathology [START_REF] Von Bubnoff | Natural killer cells in atopic and autoimmune diseases of the skin[END_REF] . More generally, genetic variants found to modulate innate immune cell parameters, in this and previous studies [START_REF] Aguirre-Gamboa | Differential Effects of Environmental and Genetic Factors on T and B Cell Immune Traits[END_REF][START_REF] Orrù | Genetic variants regulating immune cell levels in health and disease[END_REF][START_REF] Roederer | The Genetic Architecture of the Human Immune System: A Bioresource for Autoimmunity and Disease Pathogenesis[END_REF] , have been directly implicated in the aetiology of several autoimmune disorders, such as inflammatory bowel disease, ulcerative colitis and atopic dermatitis. Together, these findings illustrate the value of our approach, which mapped novel genetic associations to specific cell populations and cellular states, providing new insights into the mechanisms underlying disease pathogenesis. Further evaluations of the natural variability in cellular mediators of immunity, together with the elucidation of their environmental and genetic determinants, will facilitate a detailed dissection of the immune system in human health and disease. 2 and3 (a) Levels of association (i.e., -log 10 (q-values)) between 39 non-genetic factors and adaptive and innate cell counts, at a false discovery rate (FDR) < 1%. Except when their effects were specifically measured, immunophenotypes were regressed on age, sex, CMV status, batch effects and genomewide significant SNPs (Table 1) The effect of smoking on MFI can be found in Supplementary Fig. 19. and3) and genetic factors (independent significant and suggestive GWAS hits, a Other immunophenotypes correspond to any measured immunophenotype in the Milieu Intérieur cohort that was also significantly associated with the 594 b Replication was performed in an independent cohort of 75 European-descent Americans. Only panels 4 and 7 could be used, due to sample limitations.

Effects were in the same direction as in the primary cohort.

c P-values for biological replicates were estimated based on immunophenotypes measured from a new blood draw taken ~17 days after the initial visit, in 500 subjects of the Milieu Intérieur cohort. d EAF is the frequency of the effect allele, which was defined as the allele with a positive effect on the immunophenotype.

Cytometric analyses

Ten 8-color flow cytometry panels were developed. Details on staining antibodies can be found in Supplementary Table 2. A unique lot of each antibody was used for the entire study. Each antibody was selected and titrated as described earlier [START_REF] Hasan | Semi-automated and standardized cytometric procedures for multi-panel and multiparametric whole blood immunophenotyping[END_REF] . Gating strategies are described in Supplementary Figures 12345678910. The acquisition of cells was performed using two MACSQuant analyzers (Serial numbers 2420 & 2416), each fit with identical three lasers and ten detector optical racks (FSC, SSC and eight fluorochrome channels). Calibration of instruments was performed using MacsQuant calibration beads (Miltenyi, ref. 130-093-607). Flow cytometry data were generated using MACSQuantify™ software version 2.4.1229.1 and saved as .mqd files (Miltenyi). The files were converted to FCS compatible format and analyzed by FlowJo software version 9.5.3. A total of 313 immunophenotypes were exported from FlowJo. These included 110 cell proportions, 106 cell counts, 89 MFI and 8 ratios. We excluded from subsequent analyses all cell proportions, 35 immunophenotypes that were measured several times on different panels and were exported for quality controls, and two MFI that were measured with a problematic clone (Supplementary Table 3). A total of 166 flow cytometry measurements were thus analysed, including 76 cell counts, 87 MFI and 3 ratios (Supplementary Table 3). Problems in flow cytometry processing, such as abnormal lysis or staining, were systematically flagged by trained experimenters, which resulted in 8.70% missing data among the 166,000 measured values.

Outlier removal

Despite the exclusion of flagged problematic values, a limited number of outlier values were observed. As the goal of this study was to identify common non-genetic and genetic factors controlling immune cell levels, we removed these outlier values. Outliers were detected using a distance-based algorithm instead of a parametric method (e.g., removal based on a number of standard deviations from the mean), because of the substantial and highly variable skewness of the the total range of the sample, while a value in the lower tail was considered an outlier if that distance was more than 15% of the total range of the sample. To choose these threshold values, we simulated 10,000 log-normal distributions with a skewness similar to that of the flow cytometry measurements.

We then searched for threshold values so that simulated values outside of these ranges were observed in less than 5% of the distributions. Outliers were only looked for in the 50 highest and lowest values. This threshold was chosen to make sure that we do not miss any effect on immunophenotypes of common genetic variants (minor allele frequency>5%), or that of one of 39 continuous or common categorical non-genetic factors studied here. All values more extreme than the points labelled as outliers were also labelled outliers. A total of 24 values was removed at this stage.

Batch effects on flow cytometry measurements

Two batch effects on flow cytometry measurements were considered: the hour at which blood samples were drawn (from 8h to 11h in the morning) and the day at which samples were processed (8 to 12 samples per day, from September 2012 to august 2013). The hour of blood draw effect was evaluated with linear regression on all immunophenotypes. We observed that hour of blood draw impacts a limited number of cell counts, mainly CD16 hi NK cells (Supplementary Fig. 14a). The sampling day effect was evaluated by estimating its variance component on all immunophenotypes.

Visual inspection was used to determine whether temporal fluctuations -observed for those immunophenotypes with a large variance explained -were seasonal or not. We observed that sample processing day has a substantial impact on MFI. Fluctuations in MFI across time were strongly discontinuous, suggesting technical issues possibly related to the compensation matrix, rather than seasonal effects (Supplementary Fig. 14b).

Inclusion and imputation of candidate non-genetic factors

A large number of demographic variables were available in the Milieu Intérieur cohort [START_REF] Thomas | The Milieu Intérieur study -An integrative approach for study of human immunological variance[END_REF] . These included infection and vaccination history, childhood diseases, health-related life habits, and socio-Table 1), based on the fact that they are intrinsic factors (i.e., age, sex), or measure the exposure of individuals to exogenous factors, and thus may not be affected by the immunophenotypes themselves. These variables were filtered based on their distribution (i.e., categorical variables with only rare levels, such as infrequent vaccines, were excluded) and on their levels of dependency with other variables (e.g., height and BMI). The dependency matrix among the 39 non-genetic variables, together with batch variables, was obtained based on the generalized R 2 measures for pairwise fitted generalized linear models. If the response was a continuous variable, we used a Gaussian linear model. If the response was binary, we used logistic regression. Categorical variables were used only as predictors. Missing values were imputed using the random forest-based R package missForest [START_REF] Behar | The genome-wide structure of the Jewish people[END_REF] .

Impact of candidate non-genetic factors on immunophenotypes

To analyse the impact of non-genetic factors on immunophenotypes, we fitted a linear mixed model for each of the 166 immunophenotypes and each of the 39 non-genetic treatment variables. A total of 6,474 models were therefore fitted using the lme4 R package [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF] . All models were fitted to complete cases. Due to lack of a priori knowledge on how the non-genetic variables impact the immunophenotypes, we did not attempt to make a full causal structural equation model for all variables. Instead, we chose to keep the amount of controls in the models small to increase interpretability of the results, and to make the study easier to reproduce. We included age, sex and CMV seropositivity as fixed-effect controls for all models (Fig. 3 and Supplementary Fig. 19), except when they were the treatment variable to be tested (Fig. 2 and Supplementary Figs. 17). The intrinsic factors, i.e., age and sex, were included as covariates because they are known to have an impact on immunophenotypes [START_REF] Liston | Shaping Variation in the Human Immune System[END_REF][START_REF] Goronzy | Successful and maladaptive T cell aging[END_REF][START_REF] Aguirre-Gamboa | Differential Effects of Environmental and Genetic Factors on T and B Cell Immune Traits[END_REF][START_REF] Furman | Cytomegalovirus infection enhances the immune response to influenza[END_REF][START_REF] Pennell | Sex affects immunity[END_REF][START_REF] Furman | Systems analysis of sex differences reveals an immunosuppressive role for testosterone in the response to influenza vaccination[END_REF][START_REF] Astle | The Allelic Landscape of Human Blood Cell Trait Variation and Links to Common Complex Disease[END_REF] , as well as on many of the other environmental exposures, and are therefore possible confounders. CMV seropositivity was included because it has been shown to strongly affect some immunophenotypes [START_REF] Liston | Shaping Variation in the Human Immune System[END_REF][START_REF] Wertheimer | Aging and Cytomegalovirus Infection Differentially and Jointly Affect Distinct Circulating T Cell Subsets in Humans[END_REF][START_REF] Furman | Cytomegalovirus infection enhances the immune response to influenza[END_REF][START_REF] Brodin | Variation in the Human Immune System Is Largely Driven by Non-Heritable Influences[END_REF] . We also controlled for genome-wide significant SNPs for corresponding immunophenotypes (Table 1). Genetic variants were included to reduce the residual variance of the models and to make the inferences more robust. To correct for the batch effect related to the day of sample processing, we included it as a random effect for all models: we included imputation was performed by IMPUTE v.2 [START_REF] Howie | A Flexible and Accurate Genotype Imputation Method for the Next Generation of Genome-Wide Association Studies[END_REF] , considering 1-Mb windows and a buffer region of 1 Mb.

Out of the 37,895,612 SNPs obtained after imputation, 37,164,442 were imputed. We removed 26,005,463 imputed SNPs with information ≤ 0.8, 43,737 duplicated SNPs, 955 monomorphic SNPs, and 449,903 SNPs with missingness >5% (individual genotype probabilities < 0.8 were considered as missing data). After quality-control filters, a total of 11,395,554 high-quality SNPs were further filtered for minor allele frequencies >5%, yielding a final set of 5,699,237 SNPs for association analyses.

Genome-wide association analyses

Prior to the genome-wide association study, we transformed immunophenotypes using a different procedure than that used for the analysis of non-genetic factors. This is because we tested for association between immunophenotypes and millions of genetic variants, among which some have an unbalanced genotypic distribution (i.e., SNPs with a low minor allele frequency), which makes this analysis more sensitive to deviations from distributional assumptions. Our primary aim was therefore to use transformations that make the GWAS as robust as possible against such deviations. Also, we map loci associated with immunophenotypes based on P-values, so it was less important to keep effect sizes on the same scale, in contrast with the analysis of non-genetic factors, for which we favoured the interpretability of effect sizes. A unit value was first added to all phenotypes with zero values. The transformations were then chosen based on an AIC measure using the Jacobian-adjusted Gaussian likelihood, among three possible choices of increasing skewness: identity transformation, squareroot-tranformation and log-transformation. We kept the amount of possible transformations low to minimize the amount of added unmodelled stochasticity. The added unit value was kept only for immunophenotypes for which the log-transformation was chosen.

After transformation, a second round of outlier removal was done, to remove extreme values on the new scale. The thresholds for the lower and higher tail were 20%, obtained as for the first step of outlier removal (see description of the distance-based outlier removal algorithm above), but on the Gaussian scale. The immunophenotypes were then imputed using the missForest R package [START_REF] Behar | The genome-wide structure of the Jewish people[END_REF] , as missing data is not allowed by the subsequent analyses. We finally adjusted all immunophenotypes for the batch effect of processing days. We used the ComBat non-parametric empirical-Bayes framework [START_REF] Johnson | Adjusting batch effects in microarray expression data using empirical Bayes methods[END_REF] , instead of the mixed model described above (see section "Impact of candidate nongenetic factors on immunophenotypes" above), because the GEMMA mixed model used to conduct GWAS (see below) includes only the random effect capturing genetic relatedness. ComBat adjusts for batch effects by leveraging multivariate correlations among response variables. We did not include variables of interest in the ComBat model (none of the non-genetic variables were significantly different across sample processing days, with the exception of smoking (regression P=0.002)).

To reduce the residual variance of GWAS models and make the inferences more robust [START_REF] Mefford | The Covariate's Dilemma[END_REF] , we sought to adjust models for covariates selected among 42 variables. These included the 39 non-genetic variables (Supplementary Table 1), the hour of blood draw variable, and the two first principal components of a PCA based on genetic data (Supplementary Fig. 20b). Covariates were selected by stability selection [START_REF] Meinshausen | Stability selection[END_REF][START_REF] Shah | Variable selection with error control : another look at stability selection[END_REF] , with elastic net regression as the selection algorithm. A selection algorithm uses a cost function that drives regression parameters of non-predictive variables to zero, unlike leastsquare regressions. The elastic net method was used in particular because it has lower variance than stepwise methods and overcomes limitations of the LASSO method related to correlated variables [START_REF] Hastie | Elements of statistical learning[END_REF] .

To perform stability selection, we estimated, for each of the ∈ 1, … ,42 variables, the probability p i = P(β i = 0) that the elastic net regression parameter β i of variable i equals zero. Specifically, we first took 50 subsamples of half of the data, performed variable selection on each subsample, and estimated p i as the number of subsamples in which β i > 0, divided by the total number of subsamples. The variables were then chosen to be controls in the GWAS models by thresholding the probability ^i. It has been shown that this procedure, with the right threshold and under certain assumptions, controls the false discovery rate of selected variables [START_REF] Shah | Variable selection with error control : another look at stability selection[END_REF] . The procedure is more stable than selecting variables by, for instance, stepwise regression or elastic net without stability selection, and thus adds less unmodelled variability to the estimates. Still, because this approach does select predictive variables for each individual response variable, it adds more variance to the model selection, relative to models in which only age, sex, CMV infection and smoking would be systematically included. However, controlling for the selected variables is expected to generate more parsimonious models (i.e., the inclusion of unnecessary covariates could reduce power [START_REF] Wakefield | Bayesian and Frequentist Regression Methods[END_REF] ), and to decrease the risk of type 1 errors (e.g., some of the many rare genetic variants that are tested could associate, by chance, with an immunophenotype when the model does not fulfil inference assumptions due to a specific, unmodelled covariate).

The univariate genome-wide association study was conducted for each imputed, transformed and batch-effect corrected immunophenotype using the linear mixed model implemented in GEMMA [START_REF] Zhou | Efficient multivariate linear mixed model algorithms for genome-wide association studies[END_REF] , adjusting on selected covariates. GEMMA is an efficient mixed model that controls for genetic relatedness among donors and allows for multivariate analyses. Genetic relatedness matrices (GRM) were estimated for each chromosome separately, using the 21 other chromosomes, to exclude from the GRM estimation potentially associated SNPs (i.e., "leave-one-chromosome" approach; see [START_REF] Yang | Advantages and pitfalls in the application of mixed-model association methods[END_REF] ). A conditional GWA analysis was also carried out for each of the 14 immunophenotypes that showed the strongest genome-wide significant signals ("main immunophenotypes" in Table 1), by including as a covariate in GEMMA the genotypes of the most strongly associated variant. A multivariate GWAS was conducted on a set of 6 candidate immunophenotypes (i.e., number of HLA-DR+ memory T cells), using GEMMA linear mixed model adjusted on covariates that were selected for at least one of the six traits. For all genome-wide association analyses, a conservative genome-wide significant threshold of P<10 -10 was used, to account for testing multiple SNPs and immunophenotypes.

Power estimation

We used simulations to estimate the minimum effect of a variant that we could detect with 95% power by our GWAS. Namely, we sampled 100,000 times a SNP in our data, and simulated an immunophenotype by adding to a randomly sampled immunophenotype the effect k of that SNP, k being drawn from a uniform distribution of bounds 0 and 1 (k is expressed in unit of phenotype standard deviations, as in 'scheme 1' of ref [START_REF] Zhang | Mixed linear model approach adapted for genome-wide association studies[END_REF] ). We then ran the GEMMA mixed model on the simulated data, and estimated the probability that the variant was detected, assuming our genomewide significant threshold of P<10 -10 . We found that we have 95% power to detect a SNP with a medium effect of 0.6 phenotype standard deviation. We also confirmed empirically the power to identify medium-effect genotype-phenotype associations in the Milieu Intérieur cohort by replicating well-known genetic associations with non-immune traits, including OCA2-HERC2 genes with eye and

  -1.25]), followed by age (0.63, 95%CI: [0.42-0.95]), CMV infection (0.71, 95%CI: [0.51-0.99]), and sex (0.95, 95%CI: [0.60-1.51]). Taken together, our results indicate that genetic factors account for a substantial fraction of human variation in immune cell parameters, with their influence being stronger in innate immune cells, relative to adaptive immune cell phenotypes.

Figure 1

 1 Figure 1 Immune cell counts and cell-surface markers measured in the Milieu Intérieur cohort. Panel numbers refer to the cytometric analyses performed, grouped based on cellular lineage (Supplementary Figs. 1-10 and Supplementary Tables2 and 3). The expression of phenotypic

  ). The expression of phenotypic markers of differentiation or activation was quantified based on their mean fluorescent intensity (MFI), indicated per panel. Interconnecting lines illustrate cellular lineages or differentiation states. Red and blue squares indicate immunophenotypes significantly associated in this study with nongenetic or genetic factors, respectively.

Figure 2

 2 Figure 2 Respective effects of age, sex and CMV infection on innate and adaptive cell counts in 1,000 healthy individuals. Significant multiplicative effects (adjusted P<0.01) of (a-c) increasing age, (d-f) female sex and (g-i) CMV seropositivity on circulating levels of immune cells. (a, d, g) Effect sizes were estimated in a linear mixed model with a log-transformed immunophenotype as response, controlling for batch effects and genome-wide significant SNPs, and then transformed to the original scale. Adaptive and innate immune cells are represented in grey and black, respectively. The 99% confidence intervals (99%CIs) were false coverage-adjusted. (b, e, h) Regression lines were fitted using local polynomial regression. (b) Impact of age on naive CD8b + (in dark green) and CD4 + (in light green) T cells. (e) Impact of age and sex on the absolute count of MAIT cells. Females are represented in pink and men in blue. (h) Impact of age and CMV serostatus on CD4 + EMRA T cells. CMV+ individuals are represented in red and CMV-in orange. (c) Flow cytometry plots of naive CD8b + and CD4 + T cells for representative persons in their 20s and their 60s. (f) Flow cytometry plots of EMRA CD4 + T cells in representative CMV-and CMV+ subjects. (i) Flow cytometry plots of MAIT cells in representative woman and man. The significant effects of age, sex and CMV seropositivity on MFI can be found in Supplementary Fig.17.

17

 17

  .

Figure 3

 3 Figure 3 Effects of smoking on innate and adaptive immune cell counts in 1,000 healthy individuals.

  . (b) Significant multiplicative effects (adjusted P<0.01) of active and past smoking on circulating levels of immune cells. The multiplicative effect sizes were estimated in a linear mixed model with a log-transformed immunophenotype as response, controlling for age, sex, CMV serostatus, batch effects and genome-wide significant SNPs, and then transformed to the original data scale. 99%CIs were false coverage-adjusted. Adaptive and innate immune cells are represented in grey and black, respectively. (c) Impact of age and smoking on the number of circulating T reg cells. Brown indicates active smokers, orange indicates past smokers and yellow indicates non-smokers. Regression lines were fitted using local polynomial regression. (d) Flow cytometry plots of HLA-DR expression in T reg cells of representative non-smoker and active smoker.

Figure 4 Figure 5

 45 Figure 4 Genome-wide significant associations with 166 immunophenotypes measured in 1,000 healthy individuals. (a) Manhattan plots of genome-wide significant associations with variants acting locally (local-pQTLs, in blue) or not (cell count QTLs or trans-pQTLs, in yellow) on immunophenotypes. The gray line indicates the genome-wide significance threshold (P<10 -10 ). Zoomed Manhattan plots for all hits are shown in Supplementary Figure 21. (b) Differential expression of the CD62L protein marker in granulocytes of representative individuals homozygous for the major (T/T, in dark colors) and minor (C/C, in light colors) rs2223286 alleles. (c) Cell-specific CD62L expression is shown for age-matched individuals homozygous for the major (open distribution with solid line) or minor (shaded distribution with dotted line) rs2223286 alleles. (d) ZoomedManhattan plots of genetic associations between SNP rs2223286 in the SELL gene and cell-surface expression on CD62L in eosinophils or SELL mRNA levels in whole blood. Each point is a SNP,

Table 1 and Supplementary Table 6).

 1

	Locus	FACS panel	Immunophenotype	Other immunophenotypes a P-value	Replication P-value b	P-value for biological replicates c	Identified by a previous study	Effect size (SE)	Chr	Position	Candidate variant	Effect allele d	Other allele	EAF d Candidate gene	Distance to TSS (kb)
	1	4	CD69 in CD16 hi NK cells	CD69 + CD16 hi NK cells; CD69 in CD8a + and CD69 + CD16 +	4.8 x 10 -37 6.3 x 10 -4 2.0 x 10 -16	-	0.14 (0.01)	1 101744633 rs6693121	A	C	0.40	S1PR1	41.0
				NK cells CD16 in CD56 hi NK
	2	4	CD16 in CD16 hi NK cells	cells; HLA-DR in CD16 hi , CD8a + CD16 + NK cells CD16 + and CD69 +	3.0 x 10 -87 7.1 x 10 -7 2.6 x 10 -41 Orrù et al., Cell 2013	22.77 (1.04)	1 161507448 rs3845548	C	T	0.87 FCGR3A	12.4
	3	7	CD32 in basophils	-	1.7 x 10 -36 3.6 x 10 -7 1.6 x 10 -18	-	11.23 (0.86)	1 161653737 rs61804205	C	T	0.10 FCGR2B	20.8
	4	7	CD62L in eosinophils	CD62L in neutrophils 1.6 x 10 -35 3.7 x 10 -2 1.4 x 10 -8	-	542.78 (42.08)	1 169665632 rs2223286	C	T	0.33	SELL	0.0
	5	4	CD8a in CD69 + CD16 hi NK cells	CD8a in CD16hi, CD56 hi , CD69 + CD56 hi , CD8 + CD56 hi , CD8a + CD16 hi and HLA-DR + CD16 hi NK	5.9 x 10 -58 5.9 x 10 -2 3.4 x 10 -24 Orrù et al., Cell 2013	0.44 (0.03)	2	87026807 rs71411868	A	G	0.76	CD8A	0.0
				cells
	6	4	Number of CD8a + CD56 hi NK cells	CD56 hi NK cells; CD69+ CD56 hi NK cells; CD56 + ILC	9.1 x 10 -19 2.7 x 10 -2 2.5 x 10 -9	-	1.57 (0.18)	2 111808558 rs12986962	A	G	0.62	ACOXL / BCL2L11	0.0
	7	8	HLA-DR in cDC3	-	2.6 x 10 -11	-	3.1 x 10 -10	-	0.11 (0.02)	6	32340176 rs143655145	T	C	0.19 HLA-DRA	67.4

Table 1

 1 Genome-wide signals of association with immunophenotypes in the Milieu Intérieur cohort. 593

Acknowledgments:

This work benefited from support of the French government's Program Investissement d'Avenir, managed by the Agence Nationale de la Recherche (ANR, reference 10-LABX-69-01). We thank the Center for Translational Science, Institut Pasteur; and the OMNI Biomarker Development-Flow Cytometry Biomarker group, Genentech for their expert support. J.B. is a member of the LCCC Linnaeus Center and the ELLIIT Excellence Center at Lund University and is supported by the ELLIIT Excellence Center.

Data Availability:

The SNP array data that support the findings of this study have been deposited in the European Genome-Phenome Archive (EGA) with the accession code EGAS00001002460. The flow cytometric data can be downloaded as an R package (XX) and explored with the Shiny web application available in http://milieu_interieur_cytoGWAS.pasteur.fr.

Code Availability:

The code developed to identify non-genetic factors that impact immunophenotypes and quantify their effects has been made available in http://github.com/JacobBergstedt/mmi.

Author contributions:

Author contributions were as follows: Conceptualization, E.P., L.Q.-M., M.L.A.; Methodology, M.H., V.L., A.U., F.J., B.B., C.L., F.H., L.R., I.P., O.L., J.P.D.; Software, J.B., E.P., V.R., P.S., C.H., B.P., J.F.; Validation, M.H., V.R., F.J., Y.W.L., M.L.A.; Formal analysis, E.P., J.B., V.R., P.S., C.G., C.H., B.P., J.F.; Investigation, E.P., M.H., J.B., V.L., A.U., C.A., F.J., H.Q., M.Z., B.B., C.L., F.H., L.R., O.L., J.P.D., M.L.A.; Data curation, E.P., M.H., J.B., V.R., V.L., A.U., B.P., C.L., L.R., I.P., O.L., J.P.D.; Writing -Original Draft, E.P., J.B., C.A., D.D., M.L.A.; Writing -Review & Editing, E.P., J.B., M.H., C.A., L.R., O.L., M.F., J.F., J.P.D., L.Q.-M., M.L.A.; Supervision, E.P., M.H., M.F., J.F.,

Competing financial interests:

A.U., C.H., Y.W.L., J.H., M.Z., C.G. and M.L.A. are employees of Genentech Inc., a member of The Roche Group.

Online Methods

A summary of the Online Methods can be found in the Life Sciences Reporting Summary.

The Milieu Intérieur cohort

The 1,000 healthy donors of the Milieu Intérieur cohort were recruited by BioTrial (Rennes, France), and included 500 women and 500 men, and 200 individuals from each decade of life, between 20 and 69 years of age. Donors were selected based on stringent inclusion and exclusion criteria, detailed elsewhere [START_REF] Thomas | The Milieu Intérieur study -An integrative approach for study of human immunological variance[END_REF]

Human material and staining protocol

Whole blood samples were collected from the 1,000 healthy, fasting donors on Li-heparin, every working day from 8 to 11AM, from September 2012 to August 2013, in Rennes, France. Tracking procedures were established in order to ensure delivery to Institut Pasteur, Paris, within 6 hours of blood draw, at a temperature between 18°C and 25°C. To check the stability of our flow cytometry measures through time, a second blood sample was drawn for half of the cohort during a second visit, ~17 days on average after the first visit, ranging from 7 to 44 days. After receipt, samples were kept at room temperature prior to sample staining. Details on staining protocols can be found elsewhere [START_REF] Hasan | Semi-automated and standardized cytometric procedures for multi-panel and multiparametric whole blood immunophenotyping[END_REF] .

Reproducibility testing and assay development

For optimization studies and panel development, whole blood samples were collected from healthy volunteers enrolled at the Institut Pasteur Platform for Clinical Investigation and Access to Research Bioresources (ICAReB) within the Diagmicoll cohort. The biobank activity of ICAReB platform is NF S96-900 certified. The Diagmicoll protocol was approved by the French Ethical Committee (CPP) a constant for each day and assumed that all constants were drawn from the same normal distribution. This procedure models correlation among subjects processed during the same day. We also included the hour of blood draw as a fixed-effect control for all models.

The distributions of the immunophenotypes have variable skewness. We considered normal, lognormal and negative binomial response distributions, and chose to model all immunophenotypes as lognormal based on diagnostic plots, AIC measures and our aim to have comparable results across immunophenotypes and facilitate the interpretation of effect sizes. A total of 46 immunophenotypes had zero values. A unit value was added to those before log-transformation.

For each model, we tested the hypothesis that the regression parameter for the treatment variable was zero by an F-test with the Kenward-Roger approximation. This test has better small-and medium-sample properties than the traditional chi-square-based likelihood ratio test for mixed models [START_REF] Kenward | Small sample inference for fixed effects from restricted maximum likelihood[END_REF] and can readily be applied using the pbkrtest R package [START_REF] Halekoh | A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models -The R Package pbkrtest[END_REF] . We assumed that our sample size was large enough for this test to be appropriate and chose therefore not to do parametric bootstrapping. We considered all 6,474 tests as one multiple testing family and we used the false discovery rate (FDR) as error rate. An effect was considered significant if the adjusted P-value was smaller than 0.01. If a test was significant, confidence intervals were constructed using the profile likelihood method in such a way that the false coverage rate was controlled at a level of 0.01. The false coverage rate measures the rate of confidence intervals that do not cover the true parameter and is needed if confidence intervals are selected based on a criterion that makes these intervals especially interesting, for instance significant hypothesis tests [START_REF] Benjamini | False discovery rate-adjusted multiple confidence intervals for selected parameters[END_REF] . FCR-adjusted confidence intervals are always wider than regular intervals. All these analyses were done, and can be reproduced, with the mmi R package (http://github.com/JacobBergstedt/mmi).

Genome-wide DNA genotyping

The 1,000 subjects of the Milieu Intérieur cohort were genotyped at 719,665 SNPs by the HumanOmniExpress-24 BeadChip (Illumina, California). SNP call rate was higher than 97% in all donors. To increase coverage of rare and potentially functional variation, 966 of the 1,000 donors were also genotyped at 245,766 exonic SNPs by the HumanExome-12 BeadChip (Illumina, California). HumanExome SNP call rate was lower than 97% in 11 donors, which were thus removed from this dataset. We filtered out from both datasets SNPs that: (i) were unmapped on dbSNP138, (ii) were duplicated, (iii) had a low genotype clustering quality (GenTrain score < 0.35), (iv) had a call rate < 99%, (v) were monomorphic, (vi) were on sex chromosomes and (vii) were in Hardy-Weinberg disequilibrium (HWE P < 10 -7 ). These SNP quality-control filters yielded a total of 661,332 and 87,960 SNPs for the HumanOmniExpress and HumanExome BeadChips, respectively. The two datasets were then merged, after excluding triallelic SNPs, SNPs with discordant alleles between arrays (even after allele flipping), SNPs with discordant chromosomal position, and SNPs shared between arrays that presented a genotype concordance rate < 99%. Average concordance rate for the 16,753 SNPs shared between the two genotyping platforms was 99.9925%, and individual concordance rates ranged from 99.80% to 100%, validating that no problem occurred during DNA sample processing. The final dataset included 732,341 QC-filtered genotyped SNPs.

Genetic relatedness and structure

Possible pairs of genetically related subjects were detected using an estimate of the kinship coefficient and the proportion of SNPs that are not identical-by-state between all possible pairs of subjects, obtained with KING [START_REF] Manichaikul | Robust relationship inference in genome-wide association studies[END_REF] . Genetic structure was visualized with the Principal Component Analysis (PCA) implemented in EIGENSTRAT [START_REF] Patterson | Population Structure and Eigenanalysis[END_REF] . For comparison purposes, the analysis was performed on 261,827 independent SNPs and 1,723 individuals, which include the 1,000 Milieu Intérieur subjects together with a selection of 723 individuals from 36 populations of North Africa, the Near East, western and northern Europe [START_REF] Behar | The genome-wide structure of the Jewish people[END_REF] .

Genotype imputation

Prior to imputation, we phased the final SNP dataset with SHAPEIT2 58 using 500 conditioning haplotypes, 50 MCMC iterations, 10 burn-in and 10 pruning iterations. SNPs and allelic states were then aligned to the 1,000 Genomes Project imputation reference panel (Phase1 v3.2010/11/23). We removed SNPs that have the same position in our data and in the reference panel but incompatible alleles, even after allele flipping, and ambiguous SNPs that have C/G or A/T alleles. Genotype Ethical agreement was obtained for all gGAP donors. Samples were received at room temperature and processed 1 h after blood draw. Prior to staining, the blood was washed with PBS 1X. Except for the CD32 antibodies, the antibodies for population identification were titrated using the same clones and providers as in the primary study (Supplementary Table 2). Cell labelling were performed manually in deep-well plates. Data acquisition was performed within one hour using a calibrated FacsCantoII (Becton Dickinson). We selected panels 4 and 7 for the replication study, because 10 of the 16 GWAS hits were identified with these panels, and because of sample limitations. Immunophenotypes were transformed based on models chosen in the primary cohort. The GEMMA linear mixed model was used to test for replication, with age and sex as covariates and a GRM estimated from 1,960,432 autosomal SNPs obtained by the Illumina HumanOmni1-Quad v1.0 array.

Gene expression assays

NanoString nCounter®, a hybridization-based multiplex assay, was used to measure gene expression in non-stimulated whole blood of the 1,000 Milieu Intérieur subjects, with the Human Immunology v2 Gene Expression CodeSet. This data is described in detail in a separate work [START_REF] Piasecka | Distinctive Roles of Age, Sex and Genetics in Shaping Transcriptional Variation of Human Immune Responses to Microbial Challenges[END_REF] . Expression probes that bind to cDNAs in which at least 3 known common SNPs segregate in humans were removed from the analyses (i.e., HLA-DQB1, HLA-DQA1, HLA-DRB1, HLA-B and C8G). After quality control filters, mRNA levels were available for 986 individuals at 90 candidate genes, i.e., immunity-related genes in a 1-Mb window around the genome-wide significant and suggestive associations identified in this study. For each sample, probe counts were log 2 transformed, normalized and adjusted for batch effects. eQTL mapping was performed in a 1-Mb window around corresponding association signals, using the linear mixed model implemented in GenABEL [START_REF] Aulchenko | GenABEL : an R library for genome-wide association analysis[END_REF] . All models were adjusted on the proportion of eight major cell populations, including neutrophils, CD19 + B cells, CD4 + T cells, CD8 + T cells, CD4 + CD8 + T cells, CD4 -CD8 -T cells, NK cells, and CD14 + monocytes, to account for the effect of heterogeneous blood cell composition on gene expression.

Decomposition of the proportion of variance explained

We analysed each of the 166 batch-corrected and transformed immunophenotypes (see section "Genome-wide association analyses" of Online Methods) with a linear regression model including the four most impactful non-genetic factors (Fig. 2), i.e., age, sex, CMV seropositivity status and smoking, and both genome-wide significant (P<10 -10 ) and suggestive (P<5x10 -8 ) genetic factors. The contribution of each of these variables to the variance of each immunophenotype was calculated by averaging over the sums of squares in all orderings of the variables in the linear model, using the lmg metric in the relaimpo R package [START_REF] Grömping | Relative Importance for Linear Regression in R: The Package relaimpo[END_REF] .The averaging over orderings was done to avoid bias due to correlations among predictors.

The difference in contribution to explained variance between innate and adaptive immunophenotypes was tested using linear mixed models, where we used the log-transformed proportions of variance of each immunophenotype explained by age, sex, CMV serostatus, smoking or genetics as different response variables, and indicator variables for the immunophenotype being innate or adaptive, and being a count or an MFI. The sum of the individual contributions of associated genetic variants was used to estimate the overall contribution of genetics. Since some of the immunophenotypes are correlated, their proportion of variance explained are also correlated. To account for this, we included a random effect term whose covariance matrix was modelled as a variance component multiplied by the sample correlation matrix among the immunophenotypes. Due to the small sample size, hypothesis testing was done by building a null distribution of likelihood ratios using the parametric bootstrap. The models were fitted using the R package lme4qtl (http://github.com/variani/lme4qtl).

Because the distribution of variance explained by genetics was zero-inflated, we also tested for differences in the proportion of variance explained by non-genetic and genetic factors between innate and adaptive cell measurements with a non-parametric Mann-Whitney U test. Because the Mann-Whitney U test cannot account for correlations among immune cell measurements, we conducted this test on a subset of immunophenotypes that were selected to be uncorrelated (h<0.6 with the protoclust R package). Fifty immunophenotypes were kept, including 19 adaptive and 31 innate cell measures, among which the median Pearson's r was 0.039.