Skip to main content


We're creating a new version of this page. See preview

  • Research article
  • Open Access
  • Open Peer Review

Human newborn bacille Calmette–Guérin vaccination and risk of tuberculosis disease: a case-control study

  • 1,
  • 2,
  • 3,
  • 3,
  • 3,
  • 3,
  • 3,
  • 3,
  • 3,
  • 3,
  • 3,
  • 4,
  • 5,
  • 3,
  • 2,
  • 3Email author and
Contributed equally
BMC Medicine201614:76

  • Received: 12 November 2015
  • Accepted: 23 April 2016
  • Published:
Open Peer Review reports



An incomplete understanding of the immunological mechanisms underlying protection against tuberculosis (TB) hampers the development of new vaccines against TB. We aimed to define host correlates of prospective risk of TB disease following bacille Calmette–Guérin (BCG) vaccination.


In this study, 5,726 infants vaccinated with BCG at birth were enrolled. Host responses in blood collected at 10 weeks of age were compared between infants who developed pulmonary TB disease during 2 years of follow-up (cases) and those who remained healthy (controls).


Comprehensive gene expression and cellular and soluble marker analysis failed to identify a correlate of risk. We showed that distinct host responses after BCG vaccination may be the reason: two major clusters of gene expression, with different myeloid and lymphoid activation and inflammatory patterns, were evident when all infants were examined together. Cases from each cluster demonstrated distinct patterns of gene expression, which were confirmed by cellular assays.


Distinct patterns of host responses to Mycobacterium bovis BCG suggest that novel TB vaccines may also elicit distinct patterns of host responses. This diversity should be considered in future TB vaccine development.


  • Tuberculosis
  • Vaccine
  • Correlates of risk
  • Systems biology


Newborn vaccination with bacille Calmette–Guérin (BCG) protects infants and young children against disseminated forms of tuberculosis (TB), and to some extent, against pulmonary TB [1]. As BCG’s protective efficacy in other age groups is variable and mostly poor [1], new vaccines against TB are needed. An incomplete understanding of immunological mechanisms underlying protection against TB hampers vaccine discovery and development. Our aim was to define host correlates of prospective risk of TB disease following BCG vaccination. We proposed that this knowledge would lead to greater insight into protective mechanisms against the disease, which, in turn, would inform vaccine development.

We applied a systems biology approach to analyze unique prospective samples, and showed that in healthy infants who displayed different patterns of immune responses to the BCG vaccine, distinct mechanisms may underlie susceptibility to future TB disease, suggesting that the same intervention may not benefit all.


Study participants, blood collection and follow-up

We enrolled infants at the South African Tuberculosis Vaccine Initiative (SATVI) field site, near Cape Town, South Africa. This area has a high TB disease incidence in children < 2 years of age (>1,000/100,000/year). This study was nested within a randomized controlled trial (RCT) [2], which aimed to determine whether intradermal or percutaneous delivery of Japanese BCG at birth resulted in equivalent protection against TB. For our correlates of risk study, we enrolled 5,724 infants, randomly, from the 11,680 infants enrolled in the RCT, when they were 10 weeks of age, as previously described [3] (Fig. 1). The following were exclusion criteria: mothers known to be HIV-infected, BCG not received within 24 h of birth, significant perinatal complications, any acute or chronic disease in the infant, clinically apparent anemia in the infant, and household contact of the infant with any person with TB, or any person who was coughing.
Fig. 1
Fig. 1

Cohort of infants vaccinated with BCG at birth. At 10 weeks of age, blood was collected from HIV-negative, HIV-unexposed infants with no active or chronic illnesses (including suspected TB), and with no household exposure to an adult who was coughing, or who had TB disease. Infants were then followed for 2 years. Community-wide surveillance systems identified all children exposed to adults with TB, or children with suspected TB disease. Among these children, “definite” TB cases were defined by presence of clinical signs and symptoms of lung disease plus a sputum (induced, or early morning gastric aspirate) culture positive for M. tuberculosis, while “probable” TB cases were defined by absence of a positive culture in the presence of strong epidemiological, clinical and chest roentgenographic evidence of TB disease. Two groups of controls were identified: “household” controls were exposed to an adult in the household with TB but were found not to have TB, whereas “community” controls were infants who were either investigated for TB and found not to have disease, or infants chosen at random from the rest of the cohort. For functional assays, up to 29 definite cases and 110 controls (household controls, n = 55, and community controls, n = 55) were included in different analyses. Primary analysis of transcriptional profiling was restricted to those cases and controls included in functional assay analysis for whom PBMC were available

At 10 weeks of age, blood was collected. Infants were then followed for 2 years. Community-wide surveillance systems, described for the parent study [2], identified all children exposed to adults with TB, or children with suspected TB disease. These children were admitted to a dedicated case verification ward for evaluation [2]. Investigations included clinical and epidemiologic evaluation, two induced sputa and gastric aspirates for mycobacterial culture, chest roentgenography and a tuberculin skin test [2]. Cases included in the current study were either culture positive for Mycobacterium tuberculosis (primary analysis), or were not culture positive but there was strong clinical, roentgenographic and epidemiological evidence of TB disease (validation analysis). Infants who were evaluated in our case verification ward and were found not to have TB were included as controls (primary analysis). Some of the control infants included in the validation analysis lived in the same endemic area but did not meet the criteria for TB investigation (community controls).

Human participation was according to the US Department of Health and Human Services and good clinical practice guidelines. This included protocol approval by the University of Cape Town Research Ethics Committee (reference 016/2001), and written informed consent from the parent or legal guardians.

Gene expression analysis

PBMC isolation, storage, and in vitro incubation for RNA isolation and amplification

Peripheral blood mononuclear cells (PBMC) were isolated from peripheral blood of 10-week-old infants by density gradient centrifugation and cryopreserved. Later, PBMC were thawed, rested, and 106 cells were incubated in culture with either BCG (Danish strain, Statens Serum Institut, reconstituted as previously described [4]) at a multiplicity of infection (MOI) of 0.18, or with medium alone. PBMC were incubated for 12 h at 37 °C and 5 % CO2. Following the incubation period, cells were harvested and RNA extracted using the QIAamp RNA Blood Mini Kit (Qiagen), according to the manufacturer’s instructions. A median of 0.92 μg (range 0.24–1.56 μg) RNA was obtained from 106 PBMC. Extracted RNA was cryopreserved at −80 °C. Later, RNA was thawed and a median 339 ng (range 88–595 ng) amplified using the Illumina RNA Amplification Kit (Ambion), as previously described [4]. The median RNA yield after amplification was 15.3 μg (range 2.35–32.8 μg).

DNA microarray

Amplified RNA (750 ng per array) was hybridized to the Illumina HumanRef-8 Expression BeadChip (version 2), according to the manufacturer’s instructions. Arrays were scanned with an Illumina bead array reader confocal scanner. Array data is now housed at the gene expression database repository Gene Expression Omnibus (GEO) submission number GSE20716 [5].

Microarray data analysis

Raw Illumina probe data were exported from BeadStudio and screened for quality. Scanned array images were visually inspected for artifacts or mishandling. Diagnostic plots, including density plots, box plots and heatmaps of between-array distances, were used to assess hybridization quality within the whole data set. One array failed quality screening. Pre-processing and statistical analysis were conducted using the R statistical language and various software packages from Bioconductor (, version 3.1.2) [6]. Quantile normalization was applied, followed by a log2 transformation. Bioconductor’s genefilter package was used to filter out genes with low expression (genes retained with expression values greater than 100 or 200 units in at least two samples) and insufficient variation in expression as measured by the interquartile range (IQR; genes retained with IQR higher than 0.2, 0.3 and 0.4) across all samples tested.

Unsupervised clustering and class assignment analysis

Unsupervised cluster analysis was performed, using R’s cluster package, which utilizes hierarchical clustering. This hierarchical clustering approach has been successfully used to identify clinically relevant cancer sub-types in several studies [7, 8]. PBMC were incubated with BCG or with media alone for 12 h prior to measuring gene expression. Expression from the unstimulated condition was subtracted from that in the BCG-stimulated condition; unsupervised clustering was used to identify potential sub-groups of infants. Using the genefilter package with four different filtering criteria, we selected different numbers of probes to assess the stability of clustering: probes retained with expression values greater than 100 in at least two samples and IQR range in expression across samples > 0.2 (9,878 probes), probes retained with expression values greater than 200 in at least two samples and IQR range in expression across samples > 0.2 (7,902 probes), probes retained with expression values greater than 200 in at least two samples and IQR range in expression across samples > 0.3 (5,306 probes), and probes retained with expression values greater than 200 in at least two samples and IQR range in expression across samples > 0.4 (3,077 probes). The selected probes were used as input to test the stability and robustness of the clustering process. Two major unchanged cluster groups of 18 (cluster 1) and 15 (cluster 2) core infants were observed across the four different pre-clustering conditions. The filtering was performed before the unstimulated condition was subtracted from that in the BCG-stimulated condition. There were 11 infants with an unstable cluster membership as they switched membership from cluster 1 to cluster 2 dependent upon input gene number (Additional file 1: Figure S1). The 33 donor samples were then used to train a two-class binary classifier using pamr, a class prediction tool implemented in R, to perform a class assignment for the 11 unstable infants. The 11 remaining donors were classified, six in cluster 1 group and five in cluster 2 group. Similar class assignment was also obtained using a k-Nearest Neighbors classification algorithm, with k = 3.

Identification of differentially expressed genes

Bioconductor’s limma package was used to identify differentially expressed genes ( This approach estimates the fold change (FC) between predefined groups by fitting a linear model and using an empirical Bayes method to moderate standard errors of the estimated log-fold changes for expression values from each gene. The p values from the resulting comparison were adjusted for multiple testing according to the method of Benjamini and Hochberg [9]. This method controls the false discovery rate (FDR), which was set to 0.05 (5 % FDR). Limma analysis was then used to identify genes significantly differentially expressed between cluster 1 and cluster 2 infants (5 % FDR). We also used limma to identify genes significantly differentially expressed in a one-versus-all analysis comparing each group (cluster 1 cases, cluster 1 controls, cluster 2 cases and cluster 2 controls) with a pool of the three other groups.

Pathway analysis

Pathway analysis was performed using Gene Set Enrichment Analysis (GSEA), a non-parametric annotation-driven statistical analysis method [10]. Analysis of groups of genes encompassed into distinct signal transduction or transcriptional pathways ultimately yields greater insight into biological processes than simple exploration of lists of differentially expressed genes [11, 12]. The basic principle of GSEA is to use a list of genes ranked by FC to determine whether a known biological pathway or sets of individual genes are significantly enriched in one phenotype when compared to another. To evaluate this degree of enrichment the GSEA method calculates an enrichment score (ES) based on the Kolmogorov–Smirnov statistic. The statistical significance of a gene set’s ES is estimated by an empirical phenotype-based permutation test procedure. To account for multiple hypotheses testing, GSEA normalizes the ES for each gene set to account for variation in set sizes and calculates a FDR corresponding to each normalized ES.

We systematically tested gene sets from the Molecular Signatures Database (MSigDB; C2 collection (c2.cp.v3.0.symbols.gmt), collected from various online canonical pathways databases, including KEGG, BioCarta, Reactome and NetPath, to which we added a collection of immune-related gene sets described in Chaussabel et al. [13].

An over-representation test was performed using Fisher’s exact test. Statistical significance of p < 0.05 (5 % FDR) was used to identify gene signature enrichment from a list of genes selected as differentially expressed between two groups of samples.

Whole blood intracellular cytokine assay

At 10 weeks of age, heparinized blood was collected from all infants. One milliliter was immediately incubated with BCG (1.2 × 106 organisms/mL), as previously described [3]. Medium alone served as negative control, while staphylococcal enterotoxin B (SEB; Sigma-Aldrich, 10 μg/mL) was used as a positive control. The co-stimulatory antibodies anti-CD28 and anti-CD49d (BD Biosciences, 1 μg/mL each) were added to all conditions, as this results in enhancement of specific T cell responses. Blood was incubated for 7 h at 37 °C. At this time, plasma supernatants were harvested and cryopreserved at −80 °C. Brefeldin A was then added, followed by incubation for an additional 5 h. Cells were harvested, fixed and cryopreserved, as previously described. Later, cells were thawed, permeabilized, and stained for T cell surface markers and for intracellular TNF-α, IFN-γ and IL-2, as previously described. Frequency of CD4+ and CD8+ T cells expressing these cytokines was detected by a LSRII flow cytometer (BD Biosciences), using FACSDiva 6.1 software.

Soluble cytokine/chemokine measurement

Plasma was thawed once and levels of 29 cytokines/chemokines determined with the human cytokine LINCOplex 29-bead array assay kit (LINCO Research, Millipore). The following cytokines/chemokines were measured: interleukin-1 alpha (IL-1α), IL-1β, IL-1RA, IL-2, IL-4, IL-5, IL-6, IL-7, IL-8, IL-10, IL-12p40, IL-12p70, IL-13, IL-15, IL-17, soluble CD40 ligand (sCD40L), epidermal growth factor (EGF), eotaxin, fractalkine, granulocyte-colony stimulating factor (G-CSF), granulocyte-macrophage colony-stimulating factor (GM-CSF), interferon-gamma (IFN-γ), 10 kDa interferon gamma-induced protein (IP-10), monocyte chemoattractant protein-1 (MCP-1), macrophage inflammatory protein-1 alpha (MIP-1α), MIP-1β, transforming growth factor-alpha (TGF-α), tumor necrosis factor-alpha (TNF-α) and vascular endothelial growth factor (VEGF). The manufacturer’s instructions were followed. Fluorescence was detected using a Luminex 100 IS machine (xMAP Technology, Luminex Corporation), using Luminex software.

Lymphocyte proliferation assay

Cryopreserved PBMC were thawed in culture medium (12.5 % v/v AB+ serum in RPMI) containing 10 μg/mL DNase (Sigma-Aldrich). After washes, cells were stained with 1 μg/mL Oregon Green (Molecular Probes) and rested overnight at 37 °C in 5 % CO2. PBMC (2 × 105/well) in 200 μL culture medium (12.5 % v/v AB serum in RPMI) were incubated with BCG at a MOI of 0.01 for 6 days, at 37 °C in 5 % CO2. SEB (0.5 μg/mL) served as positive control, while medium only served as negative control. Cells were harvested with 2 mM EDTA (Sigma-Aldrich) in PBS and stained with 1 μg/mL LIVE/DEAD Fixable Violet Dead Cell Stain (Invitrogen). Cells were fixed with FACS Lysing Solution (BD Biosciences) and cryopreserved. Later, cells were thawed, washed, and stained for 1 h at 4 °C with anti-CD3 Qdot605 (clone UCHT1) and anti-CD8 Cy5.5PerCP (SK1), both from BD Biosciences. Cells were acquired as described above. Anti-mouse kappa beads (BD Biosciences), stained with the respective fluorescent-labeled antibody, were used to configure compensation settings.

Cytotoxic marker assay

Thawed PBMC (2 × 105 cells/well) in 200 μL culture medium were incubated with BCG at a MOI of 0.1, at 37 °C in 5 % CO2. Controls included SEB (0.05 μg/mL) or medium alone. After 72 h, cells were harvested, stained with the violet viability dye, fixed and cryopreserved, as described above. Later, cryopreserved, fixed cells were thawed, washed, permeabilized with Perm/Wash solution (BD Biosciences) and stained with the following antibodies for 1 h at 4 °C: anti-CD3 Qdot605 (clone UCHT1), anti-CD8 Cy5.5PerCP (SK1), anti-granzyme B Alexa 700 (GB11), anti-perforin FITC (δG9; all from BD Biosciences) and anti-granulysin PE (eBioDH2, eBioscience), acquired as described above.

Quantification of myeloid and lymphoid cell populations in PBMC

Cryopreserved PBMC from case and control infants were thawed and washed, and immediately stained. Cells were stained for 30 min at 4 °C with the following antibodies: CD45 FITC (clone 2D1), CD19 PE (clone HIB19), CD3 APC (clone UCHT1; all from BD Biosciences), CD11c PerCP-cy5.5 (clone Bu15), HLA-DR Alexa Fluor 700 (clone L243; both from BioLegend), CD14 Qdot605 (clone Tuk4) and LIVE/DEAD Fixable Violet Dead Cell Stain (both from Invitrogen). Samples were acquired on a LSRII flow cytometer and the analysis was performed using FlowJo (version 9.3.1, Tree Star). The frequencies of CD14+ monocytes, HLA-DR+ CD11c + mDC, CD3+ and CD19+ lymphocytes were expressed as relative proportions of live CD45+ leukocytes.

Analysis of non-expression data

Flow cytometric data were analyzed using FlowJo (version 8.8.4, Tree Star). Whole blood intracellular cytokine assay was analyzed as previously described [3]. For the proliferation assay, only infants who had a frequency of Oregon Greenlow CD8− T cells, following SEB incubation, that was greater than the median frequency plus 3 median absolute deviations (MAD) of the negative control, were included in the analysis. For the cytotoxic marker assay, results from infants who had a frequency of granzyme B-expressing CD8+ or CD8− T cells, following SEB incubation, that was greater than the median frequency plus 3 MAD of the negative control, were included.

Cytokine data were imported directly into Microsoft Excel for analysis, from the Luminex software. Duplicate standard curves were generated for each cytokine from the standards; these curves were used to calculate cytokine levels. Since many analytes were already secreted at broadly variable levels in unstimulated samples, we did not subtract values from unstimulated samples, and we report concentrations upon stimulation with BCG, unless otherwise indicated.

For univariate assessment of differences between infant groups, a Mann–Whitney U test was performed, using Prism (GraphPad Software, Inc.).



Our systems biology approach combined unbiased transcriptional profiling with hypothesis-driven characterization of innate and T cell responses thought to be important for TB control. Blood was collected, processed and stored at 10 weeks of age from a random 5,726 healthy infants, who belonged to a parent cohort of infants routinely vaccinated with BCG at birth [2] (refer to Fig. 1 for detail). During 2 years of follow-up, infants who developed pulmonary TB were identified (cases), as well as those who did not develop TB disease (controls) [2].

For functional assays, up to 29 definite cases and 110 controls (household controls, n = 55, and community controls, n = 55) were included in different analyses.

Primary analysis of transcriptional profiling was restricted to those cases and controls included in functional assay analysis for whom PBMC were available. Two groups of cases were identified: for primary analysis, infants with culture positive pulmonary TB (“definite” cases, n = 26 were included in analysis), and for validation analysis, infants with strong epidemiological, clinical and roentgenographic evidence of TB disease, without a positive culture—the latter is the most common clinical presentation of childhood TB (“probable” cases, n = 20). Two groups of controls were identified: for primary analysis, infants with household exposure to an adult with TB disease and found not to have TB (“household” controls, n = 18), and for validation analysis, infants similarly defined plus infants who were either investigated for TB and found not to have disease, or infants chosen at random from the rest of the cohort (“community” controls, n = 25). The definition of primary and validation cohorts differed because of limited sample sizes.

Correlates of prospective risk of TB disease could not be identified

For transcriptional profiling, PBMC isolated at 10 weeks of age from cases and controls were incubated with BCG or medium alone for 12 h. RNA was extracted and analyzed by DNA microarrays, as previously described [4]. No difference could be shown between cases and controls when gene expression from PBMC incubated with medium only was subtracted from that in PBMC incubated with BCG (Additional file 1: Figure S2, Additional file 2: Table S1; two genes were differentially expressed with a nominal p value less than 0.05); neither could differences be found when BCG-induced gene expression alone was examined independently (Additional file 2: Table S1; 25 genes were different with a nominal p value less than 0.05), nor when gene expression in medium alone was examined independently (Additional file 2: Table S1; 32 genes were different with a nominal p value less than 0.05).

Multiple additional assays were completed to compare cases and controls, focusing on immune activation thought to be important for control of TB. In whole blood incubated with BCG, there was no difference in frequency of specific Th1 cells [3], nor release of pro- and anti-inflammatory mediators (Additional file 1: Figure S3). In PBMC incubated with BCG, there was no difference in proliferation of T cells (Additional file 1: Figure S4), nor T cell expression of cytotoxic molecules (Additional file 1: Figure S5).

Finally, we tested whether frequencies of myeloid and lymphoid cellular subsets in PBMC differed, because of multiple reports of associations between ratios of these cells and various infectious diseases and vaccination outcomes [1417]. No difference was shown (Additional file 1: Figure S6).

Two distinct patterns of gene expression when cases and controls were examined together

We hypothesized that divergent host responses following BCG vaccination were responsible for the inability to show correlates of risk of TB disease on a population basis. We therefore examined gene expression of definite TB cases and household contacts in the primary cohort. Expression data from PBMC incubated with media alone was subtracted from that of PBMC incubated with BCG for 12 h. Unsupervised analysis revealed two major clusters: one containing 24 infants—cluster 1—and one containing 20 infants—cluster 2 (Fig. 2a, b). Clustering was not associated with being a case or a control (Fig. 2b, Additional file 1: Figure S2), and was stable, as altering input gene number did not impact cluster membership for the majority of infants (Additional file 1: Figure S1). We found no clinical variable that differed between the two clusters, including vaccination route, weeks of gestation, gender, ethnicity and birth weight (not shown); there were also no experimental variables that differed, including RNA quantity and quality (not shown); and there was no significant difference between the time to TB diagnosis in cluster 1 and cluster 2 case infants (data not shown).
Fig. 2
Fig. 2

Identification of two clusters of gene expression in infants. Samples from the primary cohort of definite cases and household contacts were examined together. a Unsupervised clustering analysis was performed as described in “Methods”, resulting in two major clusters of infants. The heatmap shows results of supervised hierarchical clustering analysis (Pearson correlation), using genes differentially expressed between infants belonging to the two clusters (Additional file 1: Figures S1 and S2); the normalized probe intensity obtained from PBMC stimulated with media only was subtracted from the normalized probe intensity obtained from PBMC stimulated with BCG for 12 h. b Representation of cases and controls among the two clusters. c Gene expression within biological pathways, identified by GSEA, that differed between cluster 1 (blue) and cluster 2 (green). Representative GSEA pathways were ranked by FDR q value (and p value < 0.05). d Frequencies of BCG-specific CD4+ T cells among infants from cluster 1 (blue circles) or cluster 2 (green circles). Antigen-specific cells were identified by cytokine expression following incubation of whole blood with BCG for 12 h. Bars depict medians and IQR; the Mann–Whitney U test was used to assess differences. e Quantification of molecules released in whole blood after a 7-h incubation with BCG. BCG-specific levels were calculated by subtracting levels in plasma from whole blood incubated with costimulatory antibodies alone from those in BCG-stimulated blood. Bars depict medians and IQR; the Mann–Whitney U test was used to assess differences

Overall, 461 genes were differentially expressed between the two clusters (FC > 1.3 and FDR < 5 %; Fig. 2, Additional file 3: Table S2). Differential expression of these genes could also identify two clusters of infants in the validation cohort of probable TB cases and community controls (Additional file 1: Figure S7A, Additional file 4: Table S3).

GSEA showed pathways that were unique to each cluster. The primary cohort demonstrated enrichment (FDR < 5 %) of gene sets associated with interferons, T cells, cytokines and JAK-STAT signaling in cluster 1, while gene sets associated with inflammation distinct from those in cluster 1, plus gene sets of myeloid cells, were enriched in cluster 2 (Fig. 2c, Additional file 5: Table S4). Again, the validation cohort confirmed these findings (Additional file 1: Figure S7B).

Next, we assessed whether results from hypothesis-driven cellular and soluble marker assays used to interrogate differences between cases and controls (above) would align with the distinct gene expression pathway findings. This was the case, as cluster 1 infants showed a higher frequency of BCG-specific Th1 cells and greater Th1-polarizing cytokine release, compared with cluster 2 (Fig. 2d, e).

Collectively, these results suggest that infants routinely vaccinated with BCG vaccination at birth display two distinct patterns of immune activation.

Clustering is independent of BCG stimulation ex vivo

Next, we wished to learn whether clustering of gene expression was dependent on antigenic stimulation. Using expression data from PBMC from the primary cohort incubated with media alone, we showed that 670 genes were differentially expressed between cluster 1 and 2 infants in the absence of antigen stimulation (Limma analysis, FDR < 0.05 and FC > 1.3; Additional file 1: Figure S8). Using these 670 genes, the majority of individual participants fell into the same clusters as reported above, generated when expression data from PBMC incubated from media alone was subtracted from that of PBMC incubated with BCG, with only three exceptions (Additional file 1: Figure S8). Similarly, in the validation cohort, 3,221 transcripts were differentially expressed (FDR < 0.05 and FC > 1.3) between cluster 1 and cluster 2; clustering was found to be independent of BCG antigen stimulation (data not shown).

Pathway analysis (GSEA), using results from the unstimulated samples, in the primary cohort showed that gene sets representing myeloid cells and anti-inflammatory responses were prominent in cluster 1 (Fig. 3b). In contrast, pathways commonly expressed in cluster 2 were associated with T cells (Fig. 3b). Again, similar clustering was also shown in unstimulated samples from the validation cohort (data not shown).
Fig. 3
Fig. 3

Pathway analysis, using gene expression data from unstimulated PBMC, in the two clusters of infants. Samples from the primary cohort were analyzed. A one group versus all other groups analysis approach was used to identify genes differentially expressed in each of the four groups, when compared to a pool of the three other groups (Additional file 4: Table S3). Then, genes were ranked according to the moderated t test and GSEA used to identify pathways enriched in each of the four groups. Finally, pathways uniquely and commonly expressed within each group were identified, using a FDR adjusted p value cutoff of < 0.05. The Venn diagram (a) shows both shared and unique gene sets, which are identified in (b). Gene sets are sized by normalized enrichment score (NES)

Taken together, these findings suggest that after BCG vaccination of newborns, underlying differences in the two clusters exist that are independent of specific antigen stimulation, and which may determine outcome of antigen exposure.

Distinct cellular and metabolic pathways among cases, and among controls, from each of the two clusters

Next, to learn about the possible association of clustering with risk of TB disease, we used pathway analysis (GSEA) of gene expression from unstimulated samples in the primary cohort to identify pathways uniquely expressed in each of four groups: a set of cases and controls from each cluster (Fig. 3, Additional file 5: Table S4). Pathways associated with oxidative phosphorylation and activated MAP kinase were uniquely enriched in cluster 1 cases, while those of amino acid metabolism and glycolysis were enriched in cluster 1 controls (Fig. 3, Additional file 6: Table S5). Pathways associated with T cell activation, including glucose metabolism, IL-7, IL-12 and ribosomal protein translation [18, 19], were the most prominent in cluster 2 cases (Fig. 3, Additional file 6: Table S5). We could not identify pathways unique to cluster 2 controls (Fig. 3, Additional file 6: Table S5). Similar patterns were found in the validation cohort for cluster 1 cases and controls (Additional file 1: Figure S7, Additional file 7: Table S6), but no pathway enrichment was shown in cluster 2 infants of the validation cohort.

It should be noted that in unstimulated PBMC for both the primary and validation cohorts, myeloid and inflammatory gene signatures decreased upon BCG antigen stimulation. This decline was greater in extent in cluster 1, compared with cluster 2 (Additional file 1: Figure S9); therefore, an antigen-specific decrease in myeloid and inflammatory signatures is associated with a rise in Th1 cytokine release.

Together, these findings confirm different host responses in different groups of infants following BCG vaccination, and suggest that distinct mechanisms could underlie TB disease risk in different groups of infants.

High monocyte to T cell ratios and higher specific T cell responses may be associated with risk of TB disease, in one cluster

To explore whether observed differences in gene expression were associated with relative abundance of cellular subsets within PBMC, we quantified frequencies of myeloid cells (monocytes and myeloid dendritic cells) and lymphoid cells (T and B cells) with flow cytometry (Additional file 1: Figure S6). Cluster 1 infants had higher frequencies of CD14+ monocytes and lower frequencies of CD3+ T cells in PBMC, compared with cluster 2 infants (Fig. 4a), which aligned with relatively higher myeloid pathway activation seen in cluster 1 (Fig. 3).
Fig. 4
Fig. 4

Cellular phenotype and function in case and control infants from cluster 1 and cluster 2. Samples from the primary cohort were analyzed. a Frequencies of monocytes and T cells, measured by flow cytometry (Additional file 1: Figure S6) in unstimulated PBMC from infants in cluster 1 (blue circles, n = 23) and cluster 2 (green circles, n = 18). b Comparison of the cellular phenotype of definite TB case infants from cluster 1 (red closed circles, n = 15) and from cluster 2 (red open circles, n = 10) with that of pooled household controls (black closed circles for cluster 1, n = 8, and black open circles for cluster 2, n = 8). CD14+ monocytes, CD3+ T cells and the ratio of monocytes divided by lymphocytes are shown. c Frequencies of CD4+ T cells expressing either IFN-γ, IL-2 or TNF-α, or all three cytokines together, measured by intracellular cytokine staining and flow cytometry, following whole blood stimulation with BCG for 12 h. For this analysis, definite TB case infants from cluster 1 (red closed circles, n = 12) and from cluster 2 (red open circles, n = 7) were compared to pooled household controls (black closed circles for cluster 1, n = 5, and black open circles for cluster 2, n = 8). Bars depict medians and IQR; the Mann–Whitney U test was used to assess differences in all analyses

Since our sample sizes were small, we pooled the controls from each cluster for further flow cytometric analyses. This was justified because only small differences in gene expression between the two groups of controls were shown—only seven genes were differentially expressed (Additional file 6: Table S5); in contrast, the two groups of cases showed 646 differentially expressed genes (Additional file 7: Table S6). Cases in cluster 1 displayed higher monocyte to T cell ratios, compared with the pool of controls and with cluster 2 cases (Fig. 4b). No differences were observed for the other cell subsets, including B cells (data not shown).

We then showed that cluster 1 cases had higher frequencies of BCG-specific CD4+ T cells producing type 1 cytokines (Fig. 4c), and that BCG induced greater upregulation of granulysin and perforin (Additional file 1: Figure S10A), compared with pooled controls. CD8+ and γδ TCR+ T cell responses to BCG, as well as BCG-specific CD4+ and CD8+ T cell proliferation, did not differ between groups (Additional file 1: Figure S10B and data not shown).

These results suggest that in a major cluster of BCG vaccinated infants, the presence of high frequencies of monocytes, or of a specific type 1 cellular immune response, both regarded as important for control of mycobacteria, might be associated with higher risk of developing TB disease.

M2 alternatively activated and M1 classically activated monocytes in different clusters of cases

Next, we explored gene expression pathways enriched in different clusters of case infants in greater detail. Genes in the M1 or M2 monocyte pathways showed an association with ratio of myeloid to lymphoid cellular frequency in case infants (Fig. 5a). Pathways in cluster 1 cases included oxidative phosphorylation, which is the mechanism that activated anti-inflammatory monocytes, or M2, use to gain energy for function and survival [18]. A detailed analysis revealed enrichment for expression of genes characteristic of both M1—classically activated monocytes—and M2 in cluster 1 cases, compared with the other group of cases and the controls (Fig. 5b). Cluster 2 cases showed enrichment of a sub-set of M1 genes (Fig. 5b). Cluster 1 controls showed modest enrichment of M1 genes, possibly related to higher frequencies of monocytes observed in cluster 1 infants (Figs. 4a, 5b).
Fig. 5
Fig. 5

Evidence of differential M1 and M2 monocyte activation in cases and controls from different clusters. Samples from the primary cohort were analyzed. a Heatmap showing genes associated with AMPK, oxidative phosphorylation, translation, M1 monocytes and M2 monocytes, which were differentially expressed between cluster 1 cases and the three other groups, or cluster 2 cases and the three other groups. Cases from clusters 1 and 2 are shown. To build the heatmap, the infants were first ranked by increasing expression intensity for each gene. Then, the mean-rank, across the set of genes, for each infant was used to order infants from the lowest to the highest mean-rank. A Spearman correlation was used to assess the significance (p value < 0.0008) of the association between the ordering of the infants and the monocyte to T cell ratio. b The p value table from an over-representation test performed using Fisher’s exact test to identify M1 (right column) and M2 (left column) gene signature enrichment amongst up-regulated in each of the four groups of infants. (c) Concentrations of pro- and anti-inflammatory molecules were measured by Luminex in plasma from whole blood stimulated with BCG for 7 h. Bars depict medians and IQR; the Mann–Whitney test was used to assess differences in all analyses

We evaluated whether whole blood production of cytokines associated with M1 and M2 would follow gene expression patterns. Production of EGF in unstimulated blood, and of EGF, VGEF and IL-1Ra—all associated with M2 macrophages [20]—in BCG-stimulated blood modestly confirmed gene expression patterns (Fig. 5c). Conversely, IFN-γ-driven expression of IP-10, typically from M1 monocytes, was similar across groups, again supporting gene expression results (Fig. 5c). Pro-inflammatory IL-6, MIP-1α and MIP-1β, typically from M1 monocytes, were slightly elevated in cluster 1 cases, compared with cluster 2 cases (Fig. 5c).

These findings suggest that monocyte phenotype may additionally contribute to differences in immunogenicity in the two infant clusters.


Our primary aim was to determine correlates of risk of TB disease following BCG vaccination of newborns. Despite combining unbiased and hypothesis-driven approaches, no such correlates could be identified. This may be due to the observation that, on a population level, infants vaccinated with BCG at birth demonstrated two distinct clusters of immune response, as measured by global transcriptional response to in vitro stimulation with BCG. Within these clusters, infants at risk of TB disease showed different gene expression, suggesting that unique correlates of risk may not be found within clusters. Our sample size was not adequate to do this analysis definitively, but was sufficient for learning about distinct immune responses within clusters. Immediate implications of this observation are that future studies aimed at defining correlates of risk of, or of protection against, TB disease should plan for adequate sample sizes to address heterogeneity, and that analysis strategies should account for heterogeneity.

We observed that infants with the highest or lowest ratios of monocytes to T cells were at risk of developing TB disease. We have also used full differential blood count data available from study cohorts in Durban, South Africa, to show that both high and low monocyte to lymphocyte ratios were associated with TB disease risk in HIV-exposed infants, post-partum mothers and adults initiating anti-retroviral therapy [1416]. It has previously been shown that hematopoietic stem cells can directly respond to IFN-γ mobilizing myeloid cells in response to mycobacterial infections [21]. In the absence of M. tuberculosis infection, it is therefore possible that persistence of live attenuated BCG from vaccination at birth may have driven myeloid cell expansion seen in case infants with higher monocyte to T cell ratios. We cannot exclude that host genetic determinants, other infections, other immunizations, or environmental determinants, such as nutrition, might also be important; we did not find influence of other infections or environmental co-variates in our infants. Future studies will address these variables in a systematic fashion.

We further showed that case infants with higher monocyte to T cell ratios have a transcriptional profile associated with an alternatively activated macrophage phenotype (M2), modestly confirmed by soluble marker measurement data. The M1/M2 pathways used in this study included those defined but not yet included in the MSigDB database [22], and pathways associated with the development of M2 phenotype, such as activated MAP kinase signaling and platelets [18, 23]. The latter two functions have been shown to drive monocyte differentiation into giant cells, which can suppress the mycobacteria-specific immune response [18, 23]. We have previously shown that healthy BCG vaccinated infants have M1 macrophage phenotype [4, 24], suggesting that the usual host response to BCG immunization involves classically activated monocytes (M1 activation was seen in all cases and controls in the current study, with M2 in addition in the case group referred to above). Typically, M1 macrophages are associated with killing of mycobacteria, whereas M2 macrophages are associated with tissue repair and bacterial persistence [25, 26]. Therefore, activation of M2 cells may be involved in pathogenesis of risk of TB disease.

It is widely accepted that T cells are critical for control of M. tuberculosis infection, and to prevent progression to active TB disease [2729]. The importance of CD4+ T cells is evident in persons with HIV-1-associated depletion of CD4+ T cell numbers and function; while not on antiretroviral therapy, these individuals are at > tenfold greater risk of progression from infection to disease [30]. Humans with mutations of genes in the Th1 response axis are highly susceptible to mycobacterial disease [31]. As a consequence, IFN-γ expression has been considered as a correlate of induced immunity, or “vaccine take”, in clinical trials of novel TB vaccines. However, we have recently shown that vaccine-induced antigen-specific Th1 responses did not correlate with risk of TB disease, following BCG vaccination of newborns: frequencies of BCG-specific IFN-γ-expressing, polyfunctional (co-expressing IFN-γ, TNF-α and IL-2) or IL-17-expressing CD4+ T cells, as well as IFN-γ-expressing CD8+ T cells, did not correlate with risk [3]. Surprisingly, in the current study higher frequencies of BCG-specific Th1 and polyfunctional cells were detected in a cluster of infants at risk of TB disease; the same infants also had M2 monocyte activation. Together, these findings suggest that a strong Th1 response in a context of heightened alternative activation of monocytes may be detrimental to the host.


Our findings suggest that distinct patterns of host responses to Mycobacterium bovis BCG are present in infants vaccinated with BCG. Vaccines for intracellular pathogens are designed to boost Th1 immune responses, which may be an effective strategy for immune response cluster 2 infants, but may not be effective for protection of cluster 1 infants. It is tempting to hypothesize that partial efficacy observed with candidate vaccines for TB, HIV and malaria [3234] can be ascribed to differential host responses within populations. Indeed, we have shown that the monocyte/lymphocyte ratio impacts on the efficacy of the malaria vaccine RTS,S, which protects children with lower ratio (similar to immune responses in our cluster 2), but not higher ratio (similar to immune responses in our cluster 1) in a phase II efficacy trial in Kilifi, Kenya [17].



Bacille Calmette–Guérin


Epidermal growth factor


Enrichment score


Fold change


False discovery rate


Granulocyte-colony stimulating factor


Granulocyte-macrophage colony-stimulating factor


Gene Set Enrichment Analysis






10 kDa interferon gamma-induced protein


Interquartile range


Median absolute deviation


Monocyte chemoattractant protein-1


Macrophage inflammatory protein


Multiplicity of infection


Peripheral blood mononuclear cells


Randomized controlled trial


South African Tuberculosis Vaccine Initiative


Soluble CD40 ligand


Staphylococcal enterotoxin B




Transforming growth factor-alpha


Tumor necrosis factor-alpha


Vascular endothelial growth factor



We thank study participants and their families, the community of Cape Winelands East district, and South African Tuberculosis Vaccine Initiative (SATVI) personnel.


National Institutes of Health (RO1-AI065653 to WAH), European and Developing Countries Clinical Trials Partnership (senior fellowship to WAH), Aeras, Sequella Foundation (GK) and the European Union’s Horizon 2020 project “TBVAC2020” (grant no. 643381, HAF).

The BCG study team included: Deon Minnies, Rachel Tanner, Cheong Kwet Choy Kwong Chung, Mark J Cameron, Jean-Philippe Goulet, Renaud Gaujoux, Jane Hughes, Sebastian Gelderbloem, Karen Iloni, Marie Buchanan, Linda van der Merwe, Andre Burger, Lily Denation, Sylvia Mlanjeni, Stefanie Abrahams, Ronelle Shepherd, Marijke Geldenhuys, Cathal Seoighe, Andreia Soares, Hoyam Gamieldien, Mzwandile Sidibana, Hassan Mahomed, Helen McShane and Adrian VS Hill.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

London School of Hygiene & Tropical Medicine, London, UK
Department of Pathology, Case Western Reserve University, Cleveland, OH, USA
South African Tuberculosis Vaccine Initiative (SATVI), Institute of Infectious Disease and Molecular Medicine, Division of Immunology, Department of Pathology, University of Cape Town, Cape Town, South Africa
Grupo de Inmunología Celular e Inmunogenética, Sede de Investigación Universitaria, Universidad de Antioquia, Medellín, Colombia
Public Health Research Institute, Rutgers Biomedical and Health Sciences, Newark, NJ, USA


  1. Mangtani P, Abubakar I, Ariti C, Beynon R, Pimpin L, Fine PE, et al. Protection by BCG vaccine against tuberculosis: a systematic review of randomized controlled trials. Clin Infect Dis. 2014;58(4):470–80.View ArticlePubMedGoogle Scholar
  2. Hawkridge A, Hatherill M, Little F, Goetz MA, Barker L, Mahomed H, et al. Efficacy of percutaneous versus intradermal BCG in the prevention of tuberculosis in South African infants: randomised trial. BMJ. 2008;337:a2052.View ArticlePubMedPubMed CentralGoogle Scholar
  3. Kagina BM, Abel B, Scriba TJ, Hughes EJ, Keyser A, Soares A, et al. Specific T cell frequency and cytokine expression profile do not correlate with protection against tuberculosis after bacillus Calmette-Guerin vaccination of newborns. Am J Respir Crit Care Med. 2010;182(8):1073–9.View ArticlePubMedPubMed CentralGoogle Scholar
  4. Fletcher HA, Keyser A, Bowmaker M, Sayles PC, Kaplan G, Hussey G, et al. Transcriptional profiling of mycobacterial antigen-induced responses in infants vaccinated with BCG at birth. BMC Med Genomics. 2009;2:10.View ArticlePubMedPubMed CentralGoogle Scholar
  5. Reddy TB, Riley R, Wymore F, Montgomery P, DeCaprio D, Engels R, et al. TB database: an integrated platform for tuberculosis research. Nucleic Acids Res. 2009;37(Database issue):D499–508.View ArticlePubMedPubMed CentralGoogle Scholar
  6. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5(10):R80.View ArticlePubMedPubMed CentralGoogle Scholar
  7. Beer JH, Haeberli A, Vogt A, Woodtli K, Henkel E, Furrer T, et al. Coagulation markers predict survival in cancer patients. Thromb Haemost. 2002;88(5):745–9.PubMedGoogle Scholar
  8. Sorlie T, Perou CM, Tibshirani R, Aas T, Geisler S, Johnsen H, et al. Gene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implications. Proc Natl Acad Sci U S A. 2001;98(19):10869–74.View ArticlePubMedPubMed CentralGoogle Scholar
  9. Benjamini Y, Hochberg Y. Controlling the false discovery rate—a practical and powerful approach to multiple testing. J Royal Stat Soc Ser B Methodol. 1995;57(1):289–300.Google Scholar
  10. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.View ArticlePubMedPubMed CentralGoogle Scholar
  11. Qin L, Gilbert PB, Corey L, McElrath MJ, Self SG. A framework for assessing immunological correlates of protection in vaccine trials. J Infect Dis. 2007;196(9):1304–12.View ArticlePubMedGoogle Scholar
  12. Haining WN, Wherry EJ. Integrating genomic signatures for immunologic discovery. Immunity. 2010;32(2):152–61.View ArticlePubMedGoogle Scholar
  13. Chaussabel D, Quinn C, Shen J, Patel P, Glaser C, Baldwin N, et al. A modular analysis framework for blood genomics studies: application to systemic lupus erythematosus. Immunity. 2008;29(1):150–64.View ArticlePubMedPubMed CentralGoogle Scholar
  14. Naranbhai V, Hill AV, Abdool Karim SS, Naidoo K, Abdool Karim Q, Warimwe GM, et al. Ratio of monocytes to lymphocytes in peripheral blood identifies adults at risk of incident tuberculosis among HIV-infected adults initiating antiretroviral therapy. J Infect Dis. 2014;209(4):500–9.View ArticlePubMedPubMed CentralGoogle Scholar
  15. Naranbhai V, Kim S, Fletcher H, Cotton MF, Violari A, Mitchell C, et al. The association between the ratio of monocytes:lymphocytes at age 3 months and risk of tuberculosis (TB) in the first two years of life. BMC Med. 2014;12(1):120.View ArticlePubMedPubMed CentralGoogle Scholar
  16. Naranbhai V, Moodley D, Chipato T, Stranix-Chibanda L, Nakabaiito C, Kamateeka M, et al. The association between the ratio of monocytes: lymphocytes and risk of tuberculosis among HIV-infected postpartum women. J Acquir Immune Defic Syndr. 2014;67(5):573–5.View ArticlePubMedPubMed CentralGoogle Scholar
  17. Warimwe GM, Fletcher HA, Olotu A, Agnandji ST, Hill AV, Marsh K, et al. Peripheral blood monocyte-to-lymphocyte ratio at study enrollment predicts efficacy of the RTS,S malaria vaccine: analysis of pooled phase II clinical trial data. BMC Med. 2013;11:184.View ArticlePubMedPubMed CentralGoogle Scholar
  18. O’Neill LA, Hardie DG. Metabolism of inflammation limited by AMPK and pseudo-starvation. Nature. 2013;493(7432):346–55.View ArticlePubMedGoogle Scholar
  19. Asmal M, Colgan J, Naef F, Yu B, Lee Y, Magnasco M, et al. Production of ribosome components in effector CD4+ T cells is accelerated by TCR stimulation and coordinated by ERK-MAPK. Immunity. 2003;19(4):535–48.View ArticlePubMedGoogle Scholar
  20. Allavena P, Sica A, Garlanda C, Mantovani A. The Yin-Yang of tumor-associated macrophages in neoplastic progression and immune surveillance. Immunol Rev. 2008;222:155–61.View ArticlePubMedGoogle Scholar
  21. Baldridge MT, King KY, Boles NC, Weksberg DC, Goodell MA. Quiescent haematopoietic stem cells are activated by IFN-gamma in response to chronic infection. Nature. 2010;465(7299):793–7.View ArticlePubMedPubMed CentralGoogle Scholar
  22. Martinez FO, Gordon S, Locati M, Mantovani A. Transcriptional profiling of the human monocyte-to-macrophage differentiation and polarization: new molecules and patterns of gene expression. J Immunol. 2006;177(10):7303–11.View ArticlePubMedGoogle Scholar
  23. Feng Y, Dorhoi A, Mollenkopf HJ, Yin H, Dong Z, Mao L, et al. Platelets direct monocyte differentiation into epithelioid-like multinucleated giant foam cells with suppressive capacity upon mycobacterial stimulation. J Infect Dis. 2014;210(11):1700–10.View ArticlePubMedPubMed CentralGoogle Scholar
  24. Ka MB, Daumas A, Textoris J, Mege JL. Phenotypic diversity and emerging new tools to study macrophage activation in bacterial infectious diseases. Front Immunol. 2014;5:500.View ArticlePubMedPubMed CentralGoogle Scholar
  25. Labonte AC, Tosello-Trampont AC, Hahn YS. The role of macrophage polarization in infectious and inflammatory diseases. Mol Cells. 2014;37(4):275–85.View ArticlePubMedPubMed CentralGoogle Scholar
  26. Schaale K, Brandenburg J, Kispert A, Leitges M, Ehlers S, Reiling N. Wnt6 is expressed in granulomatous lesions of Mycobacterium tuberculosis-infected mice and is involved in macrophage differentiation and proliferation. J Immunol. 2013;191(10):5182–95.View ArticlePubMedGoogle Scholar
  27. Flynn JL, Chan J. Immunology of tuberculosis. Annu Rev Immunol. 2001;19:93–129.View ArticlePubMedGoogle Scholar
  28. Dorhoi A, Reece ST, Kaufmann SH. For better or for worse: the immune response against Mycobacterium tuberculosis balances pathology and protection. Immunol Rev. 2011;240(1):235–51.View ArticlePubMedGoogle Scholar
  29. Walzl G, Ronacher K, Hanekom W, Scriba TJ, Zumla A. Immunological biomarkers of tuberculosis. Nat Rev Immunol. 2011;11(5):343–54.View ArticlePubMedGoogle Scholar
  30. Lawn SD, Myer L, Edwards D, Bekker LG, Wood R. Short-term and long-term risk of tuberculosis associated with CD4 cell recovery during antiretroviral therapy in South Africa. AIDS. 2009;23(13):1717–25.View ArticlePubMedPubMed CentralGoogle Scholar
  31. Casanova JL, Holland SM, Notarangelo LD. Inborn errors of human JAKs and STATs. Immunity. 2012;36(4):515–28.View ArticlePubMedPubMed CentralGoogle Scholar
  32. Mahomed H, Ehrlich R, Hawkridge T, Hatherill M, Geiter L, Kafaar F, et al. TB incidence in an adolescent cohort in South Africa. PLoS One. 2013;8(3):e59652.View ArticlePubMedPubMed CentralGoogle Scholar
  33. Bejon P, White MT, Olotu A, Bojang K, Lusingu JP, Salim N, et al. Efficacy of RTS, S malaria vaccines: individual-participant pooled analysis of phase 2 data. Lancet Infect Dis. 2013;13(4):319–27.View ArticlePubMedPubMed CentralGoogle Scholar
  34. Haynes BF, Gilbert PB, McElrath MJ, Zolla-Pazner S, Tomaras GD, Alam SM, et al. Immune-correlates analysis of an HIV-1 vaccine efficacy trial. N Engl J Med. 2012;366(14):1275–86.View ArticlePubMedPubMed CentralGoogle Scholar


© Fletcher et al. 2016