Association between Stress Response Genes and Features of Diurnal Cortisol Curves in the Multi-Ethnic Study of Atherosclerosis: A New Multi-Phenotype Approach for Gene-Based Association Tests

Abstract

The hormone cortisol is likely to be a key mediator of the stress response that influences multiple physiologic systems that are involved in common chronic disease, including the cardiovascular system, the immune system, and metabolism. In this paper, a candidate gene approach was used to investigate genetic contributions to variability in multiple correlated features of the daily cortisol profile in a sample of European Americans, African Americans, and Hispanic Americans from the Multi-Ethnic Study of Atherosclerosis (MESA). We proposed and applied a new gene-level multiple-phenotype analysis and carried out a meta-analysis to combine the ethnicity specific results. This new analysis, instead of a more routine single marker-single phenotype approach identified a significant association between one gene (ADRB2) and cortisol features (meta-analysis p-value=0.0025), which was not identified by three other commonly used existing analytic strategies: 1. Single marker association tests involving each single cortisol feature separately; 2. Single marker association tests jointly testing for multiple cortisol features; 3. Gene-level association tests separately carried out for each single cortisol feature. The analytic strategies presented consider different hypotheses regarding genotype-phenotype association and imply different costs of multiple testing. The proposed gene-level analysis integrating multiple cortisol features across multiple ethnic groups provides new insights into the gene-cortisol association.

About the Authors:

Zihuai He

Affiliation: Department of Biostatistics, University of Michigan, Ann Arbor, United States of America

Erin K. Payne

Affiliation: Life Sciences Program, Northrop Grumman Health Division, McLean, Virginia, United States of America

Bhramar Mukherjee

* E-mail: [email protected]

Affiliation: Department of Biostatistics, University of Michigan, Ann Arbor, United States of America

Seunggeun Lee

Affiliation: Department of Biostatistics, University of Michigan, Ann Arbor, United States of America

Jennifer A. Smith

Affiliation: Department of Epidemiology, University of Michigan, Ann Arbor, United States of America

Erin B. Ware

Affiliation: Department of Epidemiology, University of Michigan, Ann Arbor, United States of America

Brisa N. Sánchez

Affiliation: Department of Biostatistics, University of Michigan, Ann Arbor, United States of America

Teresa E. Seeman

Affiliation: Department of Medicine, Division of Geriatrics, David Geffen School of Medicine, University of California Los Angeles, Los Angeles, United States of America

Sharon L. R. Kardia

Affiliation: Department of Epidemiology, University of Michigan, Ann Arbor, United States of America

Ana V. Diez Roux

Affiliation: Department of Epidemiology, Drexel University, Philadelphia, United States of America

Introduction

Cortisol concentrations follow a strong daily pattern. They are high upon awakening, reach a maximum concentration approximately half an hour later, and slowly decrease throughout the rest of the day [1]. Additionally, cortisol concentrations increase in response to stressful situations, such as public speaking [2]. Under conditions of chronic stress, prolonged increased concentrations could have detrimental downstream physiological effects. Several population-based studies have linked daily cortisol patterns to health outcomes, including elevated blood pressure, abdominal obesity, and coronary calcification [3–5]. Cortisol concentrations and various features of the cortisol daily profile have also been linked to diabetes mellitus [6] and markers of inflammation [7]. Despite evidence of associations of various risk factors with cortisol, considerable inter-individual variability in cortisol remains unexplained. This has led to increased interest in examining genetic predictors of cortisol phenotypes. Most genetic research on cortisol to date has focused on candidate gene associations, notably the glucocorticoid receptor gene (NR3C1) and the mineralocorticoid receptor gene (NR3C2) [8–9].

The increase in cortisol concentrations in response to a stressor occurs through the activation of the hypothalamic-pituitary-adrenal (HPA) axis, which has effects on the cardiovascular system, the immune system, and metabolism [10]. The cortisol metabolic pathway suggests several key genes whose variation could affect cortisol levels. A number of genes are likely involved in the physiologic responses to stress. We selected genes known to play a role in stress reactivity and/or in responses to stress. The physiological roles of these genes were deduced from a large body of research in molecular biology, animal models, and a variety of human populations. This was done solely based on existing literature. MESA data was not used in any manner to generate the list of candidate genes. The six genes of interest for this work include the glucocorticoid receptor gene (NR3C1) and the mineralocorticoid receptor gene (NR3C2) (which act in concert to regulate cortisol responses to psychological stressors [11]), the tyrosine hydroxylase gene (TH), the alpha-2A-adrenergic receptor gene (ADRA2A), the beta-2-adrenergic receptor gene (ADRB2) (all three of which are related to cateholamine synthesis and reactivity to stress [12–15]), and the serotonin transporter gene (SLC6A4) (which may modulate various responses to stress [16]). These six genes appear to be involved in the response to psychological stressors [11, 13–16] and could therefore impact cortisol levels. However, few if any population-based studies have investigated the associations between polymorphisms in these genes and cortisol levels in multiple ethnic groups.

In this study we investigate how variation in these six genes is related to seven diurnal cortisol features within and across ethnic groups. In an attempt to better characterize the biological processes involved, the seven cortisol features are defined as mathematical functions of measurements collected six times per day over three consecutive weekdays. This is an advance over prior work that has focused on single measurements of cortisol, which are likely to be imperfect markers of the real phenotype of interest (the actual shape and level of the entire cortisol curve) [17]. Moreover, jointly modeling these cortisol features can additionally improve power over a single-feature analysis, as it allows identification of genetic associations that are present in multiple phenotypes. Several methods for analyzing multiple phenotypes jointly for single marker tests have been proposed recently [18–20].

Recent studies showed the advantages of gene-level analysis over individual single nucleotide polymorphism (SNP) analyses [21, 22]. When comparing SNP-level results across multiple ethnic groups, differences in ethnicity specific linkage disequilibrium structures may result in inconsistent findings. Since gene structure (exon and intron organization) is not likely to differ across ethnic groups, making the assessment of entire genes is likely to be a better analytic approach than analyzing individual SNPs in a multi-ethnic cohort. The gene-level analysis bypasses the problem that different tagging SNPs within gene regions may index relevant variation across ethnic groups, and reduces the burden of multiple testing. However, there is no published method so far that can handle gene-based inference with multiple correlated phenotypes.

We contrasted and compared four different analytic strategies with different hypotheses, all targeted towards gene discovery, resulting in varying number of tests: 1. Standard strategy which tests association of single SNPs (approximately 500 markers located on the six candidate genes), one at a time with each cortisol feature (~500 x 7 = 3500 tests in our study); 2. Multivariate association tests of single SNPs with multiple cortisol features (~ 500 x 1 = 500 tests); 3. Gene-level (defined as all SNPs within the gene and a 5 kilobases (kb) window up- and downstream of each gene) SNP-set analysis for a single cortisol feature (6 x 7 = 42 tests); 4. A new gene-level multiple-phenotype analysis (only 6 tests) we proposed for summarizing the gene-level results of the seven cortisol features. In each approach we test for different hypotheses, for example, marginal association of a single SNP with a single outcome, joint association of multiple outcomes with a SNP, joint association of a SNP-set with a single outcome, joint association of a SNP-set with multiple outcomes. Bonferroni correction applied to any of these set of tests controls for the family-wise error rate or the probability of rejecting at least one true null hypothesis when the global null hypotheses of no association holds. However, one has to remember that the family of tests being performed in each method and the inferential focus as stated in terms of the respective null hypotheses is different, that in turn, leads to different number of tests and different cost of multiple testing. For each analysis we carried out meta-analysis to combine the ethnicity specific results. Among the six stress response genes investigated, the new gene-level multiple-phenotype analysis identified a significant association between one gene (ADRB2) and cortisol features (meta-analysis p-value = 0.0025), while none of the meta-analysis results of the other three strategies are significant after Bonferroni correction.

Materials and Methods

Study Participants

We analyzed data from the MESA Stress Study, which is an ancillary study to the original MESA study, a longitudinal cohort study focused on investigating the early stages of atherosclerosis. The MESA Stress Study was approved by the Institutional Review Boards at UCLA and Columbia University. This research was also approved by the Health Sciences and Behavioral Sciences Institutional Review Board at the University of Michigan. The data was anonymized and de-identified prior to analysis. All procedures were carried out with the adequate understanding and written consent of the subjects. Eligible participants were 45–84 years of age and free from history of cardiovascular disease at the baseline examination (2000–2002) [23]. The MESA Stress Study took place in the context of MESA examinations 3 and 4 conducted between 2004 and 2006, and obtained detailed stress hormone data on a subsample of 1002 MESA participants recruited from the New York and Los Angeles sites. Participants for the MESA Stress Study were African Americans, European Americans, and Hispanic Americans and were enrolled as they presented for follow-up, until approximately 500 participants were recruited from each location. Of the 1002 MESA Stress Study participants, after exclusions for 1) no consent for use of genetic information, and 2) quality control processes, our resultant sample size was 950 individuals. The ethnicity specific distribution of this sample is as follows: 181 European Americans, 254 African Americans, and 515 Hispanic Americans. Table 1 presents more descriptive statistics of the study.

[Figure omitted. See PDF.]

Table 1. Characteristics of MESA Stress Study participants.

https://doi.org/10.1371/journal.pone.0126637.t001

Cortisol Features

Each MESA Stress Study participant was asked to collect six saliva samples per day at pre-specified times over three consecutive weekdays, for a maximum of 18 samples, using Salivette collection tubes. The samples were collected using the following schedule: Sample 1) upon waking and before getting out of bed; 2) 30 minutes later; 3) Around 10:00am; 4) Around 12:00 noon or before lunch, whichever came first; 5) Around 6:00pm or before dinner, whichever came first; 6) just before bed. Because earlier work has shown that the use of a time tracking device improves sample collection compliance [24], each collection tube was equipped with a time tracking device, which recorded the time when the swabs were removed for sample collection.

We explored seven features of the diurnal cortisol cycle averaged across the three days for each individual (Fig 1, Table A in S1 Table): 1. Cortisol concentration at wakeup (Wakeup); 2. Difference in cortisol concentrations between the peak and wakeup measurements (CAR); 3. The slope from 0.5 hours and 2 hours since wakeup (EDSlope); 4. The slope from 2 hours to 16 hours since wakeup (LDSlope); 5. The overall decline slope ignoring the peak value from wakeup to bedtime pooled (ODSlope); 6. The cortisol concentration at bedtime (Bedtime); 7. Standardized area under the curve for the interval 0hr-16hr since wakeup averaged across all days for an individual (AUC). Features were selected for investigation because prior work has hypothesized or demonstrated their associations with health risk factors or health outcomes [6, 7, 25, 26]. Features were modeled using all available salivary cortisol data (up to six samples per day collected over three days). Raw cortisol concentrations, measured in nmol/L, were log-transformed to more closely approximate a normal distribution [25]. Descriptive statistics can be found in Table 1. Correlation across multiple cortisol features were calculated and summarized in Table B (in S1 Table).

[Figure omitted. See PDF.]

Fig 1. Representation of the diurnal cortisol curve describing our summary features of interest.

In our study we specifically used Wakeup, Bedtime, Cortisol awakening response (CAR), Area under the curve (AUC) from 0–16 hours, Early Decline Slope, Late Decline Slope, and Overall Decline Slope.

https://doi.org/10.1371/journal.pone.0126637.g001

Genetic Information

Genotyping data included both measured and imputed SNPs available through participation in the MESA SHARe (SNP Health Association Resource) project. Under the SHARe project, genome-wide genotyping was obtained using the Affymetrix Genome-Wide Human SNP Array 6.0 platform. Imputation to HapMap was completed at the MESA Genetics Centers using the IMPUTE2 [27] program with the following reference panels: the HapMap Phase I and II, the human genome reference sequence (NCBI Build 36) [28]. Imputation for African Americans and Hispanic Americans was performed using the CEU+YRI+CHB+JPT reference panels (release #22). Imputation for European Americans was performed using only the CEU reference panel (release #24). Based on the imputed allele probabilities (AA, AB, BB), dosage data was created and used for the analysis.

We defined the six stress response gene regions as the entire gene, plus a window ±5kb around each gene. Base pair start and end positions for each gene were assigned based on annotation from the UCSC Genome Browser [29], and the analysis was restricted to common SNPs with minor allele frequency (MAF) >5%. Due to the small ethnic group specific sample sizes, a threshold for the MAF of 5% was chosen to limit the influence of unstable frequency estimates being driven by small sample sizes. Details on the chromosomal locations of each of the six genes, overall size, and the number of SNPs in each gene region by ethnic group can be found in Table 2.

[Figure omitted. See PDF.]

Table 2. Characteristics of the stress response genes.

https://doi.org/10.1371/journal.pone.0126637.t002

Statistical Analysis

We conducted ethnicity specific analysis suggested by the evaluation of inter-ethnic variability of cortisol features using ANOVA by which race has an association with cortisol features (Table C in S1 Table). We adjusted for age and sex as precision variables related to the outcome (cortisol features), consistent with previous MESA publications [30, 31, 32]. Since socio-economic status was previously shown to be associated with cortisol levels in MESA [33], we also examined study site and education as possible precision variables (these can be thought of as a proxy for socio-economic status), but found no evidence of a strong association with cortisol features. Thus these two predictors were not included in the model. We included top five ethnicity specific genetic principal components as potential confounders to account for racial differences not captured by self-reported race (in the full sample) or genetic admixture (within each ethnic group). Bonferroni correction was applied for multiple comparison adjustment in all analytic strategies introduced below.

Method 1. Single SNP Analysis for single cortisol feature.

For each ethnic group, we applied multivariable linear regression to investigate associations between each cortisol feature and each SNP separately. Given the seven cortisol features and approximately 500 SNPs in total (460 in European Americans, 526 in African Americans and 463 in Hispanic Americans), this strategy results in 7×500 = 3500 tests. The specific model is as follows:where Yik is the k-th cortisol feature (k = 1,…,K in general with K = 7 for this example) corresponding to subject i; α0k is an intercept term; Xi is a vector of non-genetic covariates; gi is a minor allele count of single SNP in the regions (gi = 0,1,2) and measurement error εik follows any distribution with mean zero and variance . αk is a vector of regression coefficients for the covariates, and βk is a regression coefficient for the SNP regressed on the k-th cortisol feature. Standard Wald test is used to test the null hypothesis H0:βk = 0.

Method 2. Single SNP analysis with joint test for multiple cortisol features.

We then applied the MultiPhen method to test for the joint association of the multiple cortisol phenotypes with a given SNP [18]. This strategy results in approximately 500 tests (approximate number of SNPs). MultiPhen inverts the regression problems so that the SNP genotypes become the outcome variable, and the phenotypes become the covariates in a proportional odds logistic regression. The model is given by:where m (m = 0,1,2) is a minor allele count of SNPgi; Yik is the k-th cortisol feature corresponding to subject i; α0gm is an intercept term; Xi is a vector of non-genetic covariates. αgm is a vector of regression coefficients for the covariates, and βgk is a regression coefficient for the k-th cortisol feature. Then it uses a likelihood ratio test with K degrees of freedom to test the global null hypothesis H0:βg1 = … = βgk = 0. Since MultiPhen cannot be directly applied to dosage data, we use best call genotype data for this analysis.

Note that though MultiPhen can aggregate the information across multiple phenotypes, the statistical justification for covariate adjustment in such a reverse regression setting with the proportional odds link function remains unclear. The method only tests for conditional association of genes and outcomes (conditional on covariates) and is not designed to unbiasedly estimate the beta coefficients or arrive at a causal interpretation. We also note that MESA is a cohort study and modeling the genes as outcomes is not interpretable and the reverse regression setting is a mere statistical artifact to construct the multivariate tests for association. One cannot interpret the beta-coefficients obtained from this model of G|Y,X the method is only used to derive joint association test for the global null hypotheses that none of the multiple outcomes are related to the single genetic marker.

Method 3. Gene-level analysis for each cortisol feature separately.

To assess the impact of the six stress response genes on cortisol features, for each ethnic group, we performed gene-level analysis using the sequence kernel association test (SKAT) separately for each pair of cortisol feature and gene region [22]. This strategy results in 7×6 = 42 tests. The model for testing genetic main effects is as follows:where Yik is the k-th cortisol feature corresponding to subject i; α0k is an intercept term; Xi is a vector of non-genetic covariates; Gi a vector of genotypes and measurement error εik follows any distribution with mean zero and variance . αk is a vector of regression coefficients for the covariates, and βk = (β1k,…,βpk)' is a vector of regression coefficients for the genotypes corresponding to the p markers on a gene. In SKAT, one assumes that each of the βjk, j = 1,…,p, independently follows an arbitrary distribution with mean zero and variance , where wjk is a weighting parameter to up-weight rarer variants. Testing is equivalent to testing H0:βk = 0, i.e., H0:βjk = 0 for j = 1,…,p. We used wjk = 1 for all SNPs since the weighting scheme is developed for rare variant association tests and hence not necessary for testing for common variant associations. The asymptotic null distribution of the resultant score statistic is a mixture of chi-squares [22].

Method 4. Gene-level multiple-phenotype analysis for jointly testing 7 cortisol features.

We proposed and applied a new multiple phenotype analysis to integrate the information contained in the seven cortisol features but still using a gene as the unit of analysis. This strategy results in only 6 tests (number of genes). The method appropriately combines p-values of the gene-level tests corresponding to each individual cortisol feature obtained in Method 3 and provides a combined p-value using Fisher’s combined probability test [34]. In order to account for the correlation among cortisol features, and thus the gene-level statistics, a permutation procedure was then followed to adjust the combined p-value across multiple features, described as follows:

Step 1: Randomly permute the subject ID linked to the genotype data to reflect the null hypothesis of no association between the cortisol features and gene regions.

Step 2: Apply SKAT (Method 3) separately for each pair of cortisol feature and gene region, adjusting for the confounders. The model with permuted genotype data is given by:where Gi* is now a vector of the genotypes corresponding to another randomly chosen subject i*.

Step 3: Use Fisher’s combined probability test to combine the test for each cortisol feature and compute a combined p-value.

Step 4: Repeat Step 1–Step 3 for 10000 times and get the distribution of the combined p-values under the null hypothesis.

Step 5: Adjust the observed p-value and determine significance using the null distribution of the combined p-value generated by the permutation procedure.

Because the permutation only destroys the association between cortisol features and gene regions but retains the correlation structure among cortisol features, this procedure can account for the correlation structure of multiple phenotypes to take advantage of the rich information contained in the seven cortisol features and leading to a single multivariate association test P-value capturing gene-level association.

A limitation of the above permutation approach is that it cannot account for the possible correlation between genotypes and covariates. Genotype-covariate association may be expected, for example when principal components are used as covariates to adjust for population stratification. As a result, such a permutation procedure may not provide entirely accurate Type I error rates in some scenarios [35].

Meta Analysis.

Given that we have a multi-ethnic sample and the association may vary across groups, we employed meta-analysis for both gene-level analysis and single SNP analysis using three different tests: 1. Weighted Z-score test [36]; 2. Fisher’s combined probability test; 3. Weighted Fisher’s combined probability test [37]. The tests allow for a meta-analysis of genetic effects across ethnic groups only using the ethnicity specific p-values.

Results

Basic demographic and health status information on the 950 Stress Study participants at baseline is provided in Table 1. Hispanic Americans represented the largest proportion of participants (54.2%), relative to the African Americans (26.7%) and European Americans (19.1%). The gender distribution was fairly equal (52.0% female). Overall, cortisol feature means varied across ethnic groups. There was a statistically significant difference in means across the ethnic groups for all cortisol features except CAR (Table C in S1 Table). We ran ANOVA to assess the impact of the covariates on cortisol features. Age and gender were significant for multiple cortisol features across ethnic groups; Site was only significant for Wakeup and ODSlope in Hispanic Americans; Education did not show significant association (Table C in S1 Table). The correlation structure of seven cortisol features can be found in Table B (in S1 Table) and we can observe some strong correlations (maximum of 0.77).

The most significant SNP in each gene identified by the single-SNP analysis (Methods 1 and 2) are summarized in Table D (in S1 Table). We note that the single-SNP meta-analysis exhibits some small p-values, but none of these p-values is significant after Bonferroni correction (Bonferroni threshold = 1.4×10−5 for single cortisol analysis, 0.0001 for MultiPhen) due to the burden of multiple comparisons and the limited sample size (950 individuals after exclusions). These results are used as supportive information to further evaluate the association via other methods.

The results for the gene-level assessment of gene-by-cortisol associations by Method 3 are presented in Table 3. Three genes, ADRA2A, ADRB2 and NR3CA, exhibited small p-values (less than 0.05) with more than one cortisol features in at least one ethnic group, but none of the single cortisol by gene analysis p-values are significant after Bonferroni correction (Bonferroni threshold = 0.0012). We then further combined the p-values obtained by Method 3 for separate cortisol features and applied the permutation procedure in Method 4 to produce a single p-value combining all seven features. After Bonferroni correction for testing of the six stress response genes (Bonferroni threshold = 0.0083), only ADRB2 exhibited significant association with cortisol features by the multi-phenotype meta-analysis using Fisher’s combined probability test (p-value = 0.0025). The association was most significant in European Americans (p-value = 0.0019) less significant in African Americans (p-value = 0.0342), and not significant in Hispanic Americans (p-value = 0.6073). We note that the p-values of the weighted tests are larger because the weighted tests assign a higher weight on Hispanic Americans. There were four instances where the association p-value was less than 0.05 among the cortisol features by meta-analysis: CAR (p-value = 0.0176), Bedtime (p-value = 0.0261), AUC (p-value = 0.0443) and ODSlope (p-value = 0.0262). The association pattern is not uniformly consistent among ethnic groups and cortisol features: ADRB2 is associated with CAR, Bedtime and LDSlope in European Americans, but only associated with AUC and ODSlope in African Americans and not significant in Hispanic Americans.

[Figure omitted. See PDF.]

Table 3. Gene based analysis testing the association between the cortisol features and common variants (minor allele frequency > 0.05) in the stress response genes using SKAT (Method 3).

https://doi.org/10.1371/journal.pone.0126637.t003

To further investigate the association between ADRB2 and multiple cortisol features at a SNP level, we examined in more detail the individual SNP-by-cortisol joint association in this region as produced by Method 2. Results for this association using meta-analysis of MultiPhen p-values are shown using LocusZoom in Fig 2 [38]. The linkage disequilibrium structure was based on European Americans which exhibited the most significant signal. SNP level results stratified by ethnic groups can be found in S1–S3 Figs. The most strongly associated SNP is rs6580583 (meta-analysis p-value using Fisher’s method = 0.0051) which had the small p-values in European Americans (p-value = 0.0121) and Hispanic Americans (p-value = 0.0140) but was not significant in African Americans (Table 4). The other two suggestive SNP findings in ADRB2 are rs17778257 and rs12654778 (meta-analysis p-value using Fisher’s method = 0.0091 and 0.0072 respectively). From the present analysis we can observe that gene-level analysis is more consistent than single SNP analysis across ethnic groups and leveraging correlation among multiple cortisol features, in turn, reducing the number of tests can aid with enhancing the power of such an analysis.

[Figure omitted. See PDF.]

Fig 2. LocusZoom plot of the association between all SNPs in the ADRB2 gene region and cortisol features.

Each SNP in the ADRB2 gene region was analyzed by MultiPhen (O’Reilly et al., 2012) stratified by ethnicity. Then the meta-analysis p-values were calculated using Fisher’s probability test and plotted using LocusZoom (Pruim et al., 2010). The linkage disequilibrium is based on the European Americans. This is labeled as Method 2 in the paper.

https://doi.org/10.1371/journal.pone.0126637.g002

[Figure omitted. See PDF.]

Table 4. P-values for the top SNP in ADRB2 gene associated with the seven cortisol features according to multi-phenotype meta analysis (Method 2).

https://doi.org/10.1371/journal.pone.0126637.t004

Discussion

This candidate gene association study investigated the associations between a’priori selected stress response genes and cortisol features in a multi-ethnic population by utilizing a gene-level multi-phenotype analysis approach. We found statistical evidence that variation in established stress response gene regions is related to features of the diurnal cortisol curve. One of the six stress response gene regions, ADRB2 had a significant association in a joint meta-analysis of multiple cortisol features. There were significant associations in European Americans, suggestive associations in African Americans, but no associations in Hispanic Americans. This gene encodes beta-2-adrenergic receptor which is in the G protein-coupled receptor superfamily [39]. It plays a role in blood pressure regulation, and influences both resting and stress-related blood pressure [15].

The SKAT methodology used for this gene-level analysis has several advantages as a set-based variance component test. First, SKAT is a powerful method that can integrate all SNPs in the region and can reduce the test degree of freedom when the SNPs are correlated due to linkage disequilibrium, which is of particular importance given the small ethnic group specific sample sizes for these analyses. Second, SKAT allows for the individual variant effects to vary from mean zero in either direction, and does not lose power when variants have different direction of effect unlike the Cohort Allelic Sum Test (CAST) [40] or Weighted Sum Statistic (WSS) [41] which have optimal power properties when effects across the SNPs are in the same direction and are of similar magnitude. Thirdly, it allows for the adjustment of covariates and is computationally efficient. The multi-phenotype analysis proposed in this paper can further integrate the rich information contained in the seven cortisol features to give a more precise description of the diurnal cortisol curve and improve the power of identifying statistical associations. We note that this gene-level multi-phenotype analysis is especially powerful when the SNPs in the gene are correlated and there exists consistent associations between more than one phenotype and the gene.

The Fisher method is a simple and powerful approach to combine p-values to increase association signals. We used the Fisher method to combine gene based association test p-values across seven cortisol features and obtained more significant results than the individual phenotype tests. Permutation was used to obtain p-values of the Fisher method, which makes it computationally expensive. An alternative approach is to extend the SKAT method to multiple phenotypes tests [42]. We considered this multi-phenotype SKAT, however, we have found that power of this approach was very sensitive to underlying assumptions regarding homogeneity or heterogeneity of genetic effects across phenotypes. In contrast, the Fisher method provided robust power regardless of underlying genetic effects model.

The permutation procedure we used for the Fisher method can preserve the correlation among cortisol features and the correlation between cortisol features and confounders. This strategy has been used in genetics and other areas of biological sciences and provided a reasonable performance for type I error correction [43]. However, in presence of genotype-covariate association that is ignored in our permutation approach, we may incur slightly inflated Type 1 error as mentioned in Remark 2.

There are three main limitations to this work. The first is a design limitation due the use of HapMap imputed variants, which are not functional SNPs. However, as the HapMap tagging SNPs may be in linkage disequilibrium with causal SNPs they are still useful for identifying genomic regions of potential interest. Second, compliance with cortisol sampling protocols is necessary for estimating reliable cortisol features [44, 45]. Compliance with taking samples within 10 minutes the requested times was greatest for wakeup (68%) and bedtime (75%) collections, and poorest during the middle of the day, ranging from 43%-57%. Stability of the cortisol features is of particular importance for genetic analyses, compared to other MESA cortisol work, as the effect estimates of individual variants are expected to be modest and large variation in features estimates could mask true associations. A third limitation is that, in absence of other information about genes most relevant to HPA activity, we selected a set of candidate genes involved in stress responsivity generally. Several of these (including ADRB2) are more involved in sympathetic activity than HPA activity per se; however both systems are interrelated. In this sense our analyses serve to illustrate a useful methodologic approach to examine genetic predictors of cortisol than to identify the most relevant genes. Future analyses should include a broader set of potentially relevant genes.

Despite the limitations, this work is novel in the ability to integrate the information of multiple SNPs in the region, multiple cortisol features and multiple ethnic groups, which was possible through the use of the innovative SKAT and multi-phenotype methodologies as well as the unique, highly detailed cortisol phenotype information. The gene-level multi-phenotype approach allows us to address the concern that individual SNPs may not replicate across ethnic groups as well as cortisol features due to differences in underlying patterns of linkage disequilibrium, differences in allele frequencies [46, 47], and heterogeneous effects on different cortisol features. The multiple phenotype gene-based analyses presented here provide new insight into the relationship between stress response genes and cortisol features and illustrate the advantage of integrating information in both the phenotype and genotype space by considering different hypotheses regarding marginal or joint association between a set of cortisol features and a set of SNPs, rather than analyzing each SNP and each cortisol feature in isolation. The multivariate association analysis at a gene level can enhance power of gene discovery tests by considering the correlation structure of the SNPs and the correlation structure between the multiple cortisol features.

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. LocusZoom plot of the association between all SNPs in the ADRB2 gene region and cortisol features in European Americans.

Each SNP in the ADRB2 gene region was analyzed by MultiPhen (O’Reilly et al., 2012). The p-values were plotted using LocusZoom (Pruim et al., 2010).

https://doi.org/10.1371/journal.pone.0126637.s001

(TIFF)

S2 Fig. LocusZoom plot of the association between all SNPs in the ADRB2 gene region and cortisol features in African Americans.

Each SNP in the ADRB2 gene region was analyzed by MultiPhen (O’Reilly et al., 2012). The p-values were plotted using LocusZoom (Pruim et al., 2010).

https://doi.org/10.1371/journal.pone.0126637.s002

(TIFF)

S3 Fig. LocusZoom plot of the association between all SNPs in the ADRB2 gene region and cortisol features in Hispanic Americans.

Each SNP in the ADRB2 gene region was analyzed by MultiPhen (O’Reilly et al., 2012). The p-values were plotted using LocusZoom (Pruim et al., 2010).

https://doi.org/10.1371/journal.pone.0126637.s003

(TIFF)

S1 Tables. Supplementary Tables A—D.

https://doi.org/10.1371/journal.pone.0126637.s004

(DOCX)

Acknowledgments

MESA and the MESA SHARe project are conducted and supported by the National Heart, Lung, and Blood Institute (NHLBI) in collaboration with MESA investigators. Genotyping was performed at Affymetrix (Santa Clara, California, USA) and the Broad Institute of Harvard and MIT (Boston, Massachusetts, USA) using the Affymetrix Genome-Wide Human SNP Array 6.0.

Author Contributions

Conceived and designed the experiments: AVD SLRK BM. Analyzed the data: ZH EKP. Contributed reagents/materials/analysis tools: SL EBW JAS BNS TES. Wrote the paper: ZH EKP BM SL.

Citation: He Z, Payne EK, Mukherjee B, Lee S, Smith JA, Ware EB, et al. (2015) Association between Stress Response Genes and Features of Diurnal Cortisol Curves in the Multi-Ethnic Study of Atherosclerosis: A New Multi-Phenotype Approach for Gene-Based Association Tests. PLoS ONE 10(5): e0126637. https://doi.org/10.1371/journal.pone.0126637

References

1. Kirschbaum C, Hellhammer DH (1989) Salivary cortisol in psychobiological research: an overview. Neuropsychobiology 22(3):150–169. pmid:2485862

2. Young E, Abelson J, Lightman S (2004) Cortisol pulsatility and its role in stress regulation and health. Front.Neuroendocrinol. 25(2):69–76. pmid:15571755

3. Bjorntorp P, Holm G, Rosmond R (1999) Hypothalamic arousal, insulin resistance and Type 2 diabetes mellitus. Diabetic medicine 16(5):373–383. pmid:10342336

4. Rosmond R, Chagnon YC, Holm G, Chagnon M, Pérusse L, Lindell K, et al. (2000) A glucocorticoid receptor gene marker is associated with abdominal obesity, leptin, and dysregulation of the hypothalamic-pituitary-adrenal axis. Obes. Res. 8(3):211–218. pmid:10832763

5. Matthews K, Schwartz J, Cohen S, Seeman T (2006) Diurnal cortisol decline is related to coronary calcification: CARDIA study. Psychosomatic medicine 68(5):657–661. pmid:17012518

6. Champaneri S, Xu X, Carnethon MR, Bertoni AG, Seeman T, Diez Roux AV, et al. (2011) Diurnal salivary cortisol and urinary catecholamines are associated with diabetes mellitus: the Multi-Ethnic Study of Atherosclerosis. Metabolism.

7. Desantis AS, Diez roux AV, Hajat A, Aiello AE, Golden SH, Jenny NS, et al. (2011) Associations of salivary cortisol levels with inflammatory markers: The Multi-Ethnic Study of Atherosclerosis. Psychoneuroendocrinology.

8. Joels M, Karst H, DeRijk R, de Kloet ER (2008) The coming out of the brain mineralocorticoid receptor. Trends Neurosci 31(1):1–7.

9. DeRijk R, de Kloet ER (2005) Corticosteroid receptor genetic polymorphisms and stress responsivity. Endocrine 28(3):263–270. pmid:16388115

10. Hyman SE (2009) How adversity gets under the skin. Nat.Neurosci 12(3):241–243. pmid:19238182

11. DeRijk R and de Kloet ER (2005) Corticosteroid receptor genetic polymorphisms and stress responsivity. Endocrine 28(3):263–270. pmid:16388115

12. Rao F, Zhang L, Wessel J, Zhang K, Wen G, Kennedy BP, et al. (2007) Tyrosine hydroxylase, the rate-limiting enzyme in catecholamine biosynthesis: discovery of common human genetic variants governing transcription, autonomic activity, and blood pressure in vivo. Circulation 116(9):993–1006. pmid:17698732

13. Barbeau P, Litaker MS, Jackson RW, Treiber FA (2003) A tyrosine hydroxylase microsatellite and hemodynamic response to stress in a multi-ethnic sample of youth. Ethn Dis. 13(2):186–192. pmid:12785414

14. Finley JC Jr, O'Leary M, Wester D, MacKenzie S, Shepard N, Farrow S, et al. (2004) A genetic polymorphism of the alpha2-adrenergic receptor increases autonomic responses to stress. J Appl Physiol. 96(6):2231–2239. pmid:14742450

15. Li GH, Faulhaber HD, Rosenthal M, Schuster H, Jordan J, Timmermann B, et al. (2001) Beta-2 adrenergic receptor gene variations and blood pressure under stress in normal twins. Psychophysiology 38(3):485–489. pmid:11352137

16. Wüst S, Kumsta R, Treutlein J, Frank J, Entringer S, Schulze TG, et al. (2009) Sex-specific association between the 5-HTT gene-linked polymorphic region and basal cortisol secretion. Psychoneuroendocrinology.

17. Bolton JL, Hayward C, Direk N, Lewis JG, Hammond GL, Hill LA, et al. (2014) Genome Wide Association Identifies Common Variants at the SERPINA6/SERPINA1 Locus Influencing Plasma Cortisol and Corticosteroid Binding Globulin. PLoS genetics 10(7), e1004474. pmid:25010111

18. O’Reilly PF, Hoggart CJ, Pomyen Y, Calboli FC, Elliott P, Jarvelin MR, et al. (2012) MultiPhen: joint model of multiple phenotypes can increase discovery in GWAS. PLoS One 7(5), e34861. pmid:22567092

19. Stephens M (2013) A unified framework for association analysis with multiple related phenotypes. PLoS one 8(7), e65245. pmid:23861737

20. Bhattacharjee S, Rajaraman P, Jacobs KB, Wheeler WA, Melin BS, Hartge P, et al. (2012) A subset-based approach improves power and interpretation for the combined analysis of genetic association studies of heterogeneous traits. The American Journal of Human Genetics 90(5), 821–835. pmid:22560090

21. Wu MC, Kraft P, Epstein MP, Taylor DM, Chanock SJ, Hunter DJ, et al. (2010) Powerful SNP-set analysis for case-control genome-wide association studies. The American Journal of Human Genetics 86(6), 929–942. pmid:20560208

22. Wu MC, Lee S, Cai T, Li Y, Boehnke M, Lin X (2011) Rare-variant association testing for sequencing data with the sequence kernel association test. The American Journal of Human Genetics 89(1), 82–93. pmid:21737059

23. Bild DE, Bluemke DA, Burke GL, Detrano R, Diez Roux AV, Folsom AR, et al. (2002) Multi-ethnic study of atherosclerosis: objectives and design. Am.J.Epidemiol. 156(9):871–881. pmid:12397006

24. Kudielka B, Broderick J, Kirschbaum C (2003) Compliance with saliva sampling protocols: electronic monitoring reveals invalid cortisol daytime profiles in noncompliant subjects. Psychosom.Med. 65(2):313–319. pmid:12652000

25. DeSantis AS, Diez Roux AV, Hajat A, Golden SH, Jenny NS, Sanchez BN, et al. (2011) Associations of salivary cortisol levels with metabolic syndrome and its components: the multi-ethnic study of atherosclerosis. J.Clin.Endocrinol.Metab 96(11):3483–3492. pmid:21880797

26. Champaneri S, Xu X, Carnethon MR, Bertoni AG, Seeman T, DeSantis AS, et al. (2012) Diurnal salivary cortisol is associated with body mass index and waist circumference: The multi-ethnic study of atherosclerosis. Obesity (Silver Spring)

27. IMPUTE2. Available: http://mathgen.stats.ox.ac.uk/impute/impute_v2.html#home.

28. Gibbs RA, Belmont JW, Hardenbol P, Willis TD, Yu F, Yang H, et al. (2003) The international HapMap project. Nature 426(6968), 789–796. pmid:14685227

29. Kent WJ, Sugnet CW, Furey TS, Roskin KM, Pringle TH, Zahler AM, et al. (2002) The human genome browser at UCSC. Genome Res. 12(6):996–1006. pmid:12045153

30. Do DP, Diez Roux AV, Hajat A, Auchincloss AH, Merkin SS, Ranjit N, et al. (2011) Circadian rhythm of cortisol and neighborhood characteristics in a population-based sample: the Multi-Ethnic Study of Atherosclerosis. Health & place, 17(2), 625–632.

31. DeSantis AS, Diez Roux AV, Hajat A, Golden SH, Jenny NS, Sanchez BN, et al. (2011). Associations of salivary cortisol levels with metabolic syndrome and its components: the multi-ethnic study of atherosclerosis. The Journal of Clinical Endocrinology & Metabolism, 96(11), 3483–3492.

32. DeSantis AS, Diez Roux AV, Hajat A, Aiello AE, Golden SH, Jenny NS, et al. (2012). Associations of salivary cortisol levels with inflammatory markers: The Multi-Ethnic Study of Atherosclerosis. Psychoneuroendocrinology, 37(7), 1009–1018. pmid:22178583

33. Hajat A, Diez Roux AV, Franklin TG, Seeman T, Shrager S, Ranjit N, et al. (2010). Socioeconomic and race/ethnic differences in daily salivary cortisol profiles: the multi-ethnic study of atherosclerosis. Psychoneuroendocrinology, 35(6), 932–943. pmid:20116177

34. Fisher RA (1925) Statistical methods for research workers. Genesis Publishing Pvt Ltd..

35. Bůžková P, Lumley T, Rice K (2011) Permutation and parametric bootstrap tests for gene—gene and gene—environment interactions. Annals of human genetics 75(1), 36–45. pmid:20384625

36. Liu L, Sabo A, Neale BM, Nagaswamy U, Stevens C, Lim E, et al. (2013) Analysis of rare, exonic variation amongst subjects with autism spectrum disorders and population controls. PLoS genetics 9(4), e1003443. pmid:23593035

37. Liu DJ, Peloso GM, Zhan X, Holmen OL, Zawistowski M, Feng S, et al. (2014) Meta-analysis of gene-level tests for rare variant association. Nature genetics 46(2), 200–204. pmid:24336170

38. Pruim RJ, Welch RP, Sanna S, Teslovich TM, Chines PS, Gliedt TP, et al. (2010) LocusZoom: regional visualization of genome-wide association scan results. Bioinformatics 26(18), 2336–2337. pmid:20634204

39. Nygaard R, Zou Y, Dror RO, Mildorf TJ, Arlow DH, Manglik A, et al. (2013) The dynamic process of β2-adrenergic receptor activation. Cell 152(3), 532–542. pmid:23374348

40. Morgenthaler S, Thilly WG (2007) A strategy to discover genes that carry multi-allelic or mono-allelic risk for common diseases: a cohort allelic sums test (CAST). Mutat.Res. 615(1–2), 28–56.

41. Madsen BE, Browning SR (2009) A groupwise association test for rare mutations using a weighted sum statistic. PLoS Genet. 5(2), e1000384. pmid:19214210

42. Maity A, Sullivan PF, Tzeng JI (2012) Multivariate Phenotype Association Analysis by Marker-Set Kernel Machine Regression. Genetic epidemiology 36(7), 686–695. pmid:22899176

43. Bush WS, Moore JH (2012) Genome-wide association studies. PLoS computational biology 8(12), e1002822. pmid:23300413

44. Kudielka B, Broderick J, Kirschbaum C (2003) Compliance with saliva sampling protocols: electronic monitoring reveals invalid cortisol daytime profiles in noncompliant subjects. Psychosom. Med. 65(2):313–319. pmid:12652000

45. Broderick JE, Arnold D, Kudielka BM, Kirschbaum C (2004) Salivary cortisol sampling compliance: comparison of patients and healthy volunteers. Psychoneuroendocrinology 29(5):636–650. pmid:15041086

46. Lander ES, Schork NJ (1994) Genetic dissection of complex traits. Science 265(5181):2037–2048. pmid:8091226

47. Bryc K, Velez C, Karafet T, Moreno-Estrada A, Reynolds A, Auton A, et al. (2010) Colloquium paper: genome-wide patterns of population structure and admixture among Hispanic/Latino populations. Proc.Natl.Acad.Sci. U.S.A. 107 Suppl 2:8954–8961. pmid:20445096