Participants

The primary sample used involved participants from UK Biobank, an open-access resource established to examine the determinants of disease in middle-aged and older adults living in the United Kingdom55. Recruitment to UK Biobank occurred between 2006 and 2010, targeting community-dwelling individuals from both urban and rural environments across a broad range of socioeconomic circumstances. A total of 502,655 participants were assessed at baseline on a range of cognitive and other psychological measures, physical and mental health, and their SEP. They donated a number of biological samples, including DNA for genotyping. In order to reduce the effects of population stratification, only participants from a single-ancestry group, those of White British ancestry, were included in the analysis. High-quality genotyping was performed by UK Biobank on 332,050 participants. Ethical approval for UK Biobank was received from the Research Ethics Committee (REC reference 11/NW/0382). This work was conducted under UK Biobank application 10279.

Phenotype description

A total of 332,050 participants had genotype data and data on their level of household income. Self-reported household income was collected using a 5-point scale corresponding to the total household income before tax, 1 being less than £18,000, 2 being £18,000–£29,999, 3 being £30,000–£51,999, 4 being £52,000–£100,000 and 5 being greater than £100,000. This 5-point scale was analyzed by treating the categories of income as a continuous variable. Participants were removed from the analysis if they answered “do not know” (n = 12,721), or “prefer not to answer” (n = 31,947). This left a total number of 286,301 participants (138,425 male) aged 39–73 years (mean = 56.5, SD = 8.0 years) with genotype data who had reported, between 1 and 5, their level of household income.

UK Biobank genotyping

Full details of the UK Biobank genotyping procedure have been made available56. In brief, two custom genotyping arrays were used to genotype 49,950 participants (UK BiLEVE Axiom Array) and 438,427 participants (UK Biobank Axiom Array)56,57. Genotype data on 805,426 markers were available for 488,377 of the individuals in UK Biobank. Imputation to the Haplotype Reference Consortium (HRC) reference panel leads to 39,131,578 autosomal SNPs being available for 487,442 participants56. Allele frequency checks58 against the HRC59 and 1000G60 site lists were performed, and variants with minor allele frequencies (MAFs) differing more than ±0.2 from the reference sets were removed.

Additional quality control steps were conducted and described previously8,34. These included the removal of those with non-British ancestry based on self-report and a principal components analysis, as well as those with extreme scores based on heterozygosity and missingness. Individuals with neither XX nor XY chromosomes, along with those individuals whose reported sex was inconsistent with genetically inferred sex, were also removed. Finally, individuals with more than ten putative third-degree relatives (identified by Bycroft et al.56 by estimating the kinship coefficients for all pairs of samples using the software KING61) were also removed. Following these exclusions, a sample of 408,095 individuals remained. Using GCTA–GREML on 131,790 reportedly related participants62, related individuals were removed based on a genetic relationship threshold of 0.025. Following this quality control, household income data and genetic data were available on 286,301 participants. Following association analysis, SNPs with an MAF < 0.0005, and an imputation quality score <0.1 were removed. Finally, only biallelic SNPs were retained, resulting in 18,485,882 autosomal SNPs.

GWAS in the UK Biobank sample

The level of household income as measured on the 5-point scale was subjected to a regression using income as the outcome, as has been conducted previously7, and 40 genetic principal components (to control for population stratification), genotyping array, batch, age and sex as predictors. The residuals from this model were then used in a GWAS assuming an additive genetic model as implemented in BGENIE56.

Functional annotation and loci discovery

Genomic risk loci were derived using the summary data from the data set of household income derived in UK Biobank, using FUnctional Mapping and Annotation of genetic associations (FUMA)63. First, independent significant SNPs were defined using a P value cutoff of genome-wide significance (P < 5 × 10−8), as well as being independent from each other (r2 < 0.6) within a 1-mb window. Second, SNPs that were in LD with any independent SNP (r2 ≥ 0.6) and within a 1-mb window in addition to being in the HRC genomes reference panel with an MAF >0.001, were included for further annotation. Third, lead SNPs were identified using the independent significant SNPs as defined above. Lead SNPs were a subset of the independent significant SNPs that were in LD with each other at r2 <0.1, with a 1-mb window. Fourth, genomic risk loci were created by merging lead SNPs if they were closer than 250 kb apart. This means that a genomic risk locus could contain multiple independent significant SNPs and multiple lead SNPs. Finally, all SNPs in LD of r2 ≥0.6 with one of the independent significant SNPs formed the border, or edge, of the genomic risk loci.

The lead SNPs and those in LD with the lead SNPs were then mapped to genes based on their functional consequences, as described using ANNOVAR64 and the Ensemble genes build 85. Intergenic SNPs were annotated as the two closest flanking genes, which can result in them being assigned to multiple genes.

Gene mapping

Three strategies were used to link the income-associated independent genomic loci to genes. First, positional mapping65 was used to map SNPs to genes based on physical distance. SNPs were mapped to genes if they were within 10 kb of a known protein gene found in the human reference assembly (hg19).

Second, eQTL mapping was carried out by mapping SNPs to genes if allelic variation at the SNP is associated with expression levels of a gene. For eQTL mapping, information on 45 tissue types from three databases (GTEx v7, Blood eQTL browser and BIOS QTL browser) based on cis-QTLs was used and SNPs were mapped to genes up to 1 Mb away. A false discovery rate (FDR) of 0.05 was used as a cutoff to define significant eQTL associations.

Finally, chromatin interaction mapping was carried out to map SNPs to genes when there is a three-dimensional DNA–DNA interaction between the SNP and gene. No distance boundary was applied as chromatin interactions can be long-ranging and span multiple genes. Hi-C data of 14 tissue types were used for chromatin interaction mapping66. In order to reduce the total number of genes mapped using chromatin interactions and to increase the likelihood that those mapped are biologically relevant, an additional filter was added. We only retained interaction-mapped genes if one region involved with the interaction overlapped with a predicted enhancer region in any of the 111 tissue/cell types found in the Roadmap Epigenomics Project67, and the other region was located in a gene promoter region (i.e. 250 bp upstream and 500 bp downstream of the transcription start site and also predicted to be a promoter region by the Roadmap Epigenomics Project67). An FDR of 1 × 10−5 was used to define a significant interaction.

Gene-based GWAS

Gene-based analyses have been shown to increase the power to detect association due to the multiple testing burden being reduced, in addition to the effects of multiple SNPs being combined68. Gene-based GWAS was conducted using MAGMA13. All SNPs located within protein-coding genes were used to derive a P value describing the association found with household income. The NCBI build 37 was used to determine the location and boundaries of 18,782 autosomal genes, and linkage disequilibrium within and between genes was gauged using the HRC panel. In order to control for multiple testing, a Bonferroni correction was applied using each gene as an independent statistical unit (0.05/18,782 = 2.66 × 10−6). The gene-based statistics derived using MAGMA were then used to conduct the gene-set analysis, the gene-property analyses and the cell-type enrichment analysis.

Gene-set analysis

In order to understand the biological systems vulnerable to perturbation by common genetic variation, a competitive gene-set analysis was performed. Competitive testing, conducted in MAGMA13, examines if genes within the gene set are more strongly associated with the trait of interest than other genes, and differs from self-contained testing by controlling for type 1 error rate as well as being able to examine the biological relevance of the gene set under investigation69.

A total of 10,891 gene sets (sourced from Gene Ontology70, Reactome71 and MSigDB72) were examined for enrichment of household income. A Bonferroni correction was applied to control for the multiple tests performed on the 10,891 gene sets available for analysis.

Gene-property analysis

In order to identify the relative importance of particular tissue types, which may indicate the intermediary biological phenotypes that might act between genetic variation and SEP outcomes, a gene-property analysis was conducted using MAGMA. The goal of this analysis was to determine if, in 30 broad tissue types and 53 specific tissues, tissue-specific differential expression levels were predictive of the association of a gene with household income. Tissue types were taken from the GTEx v7 RNA-Seq (RNA-sequencing) database73 with expression values being log 2 transformed with a pseudocount of 1 after Winsorising at 50 with the average expression value being taken from each tissue. Multiple testing was accounted for using Bonferroni correction.

An additional gene-property analysis was conducted to determine if transcription in the brain at 11 developmental stages14, or across 29 different age groupings14, was associated with a gene’s link to household income. This RNA-Seq GEncode v10 summarized to gene data was accessed using the following link: http://www.brainspan.org/api/v2/well_known_file_download/267666525. The detailed descriptions of the normalization processes used can be found in the technical white paper at http://help.brain-map.org/download/attachments/3506181/Transcriptome_Profiling.pdf?version=1&modificationDate=1382036562736&api=v2, where a total of 524 samples were available for analysis. The developmental stages were assigned to each of the two groups (11 developmental stages and 29 age groupings) based on the age of the sample. The groupings of 25 post-conception weeks (pcw) and 35 pcw were excluded from the age groups as they contained fewer than three samples. Next, the 52,376 annotated genes were filtered so that the average reads per kilobase (RPKM) are >1. This was performed in the developmental group and in the age group separately. This resulted in 19,601 genes for the developmental- stage groupings and 21,001 genes for the age groupings. RPKM was then winsorized at 50 (RPKM >50 was replaced with 50). Then, the average of log-transformed RPKM with a pseudocount 1 (log 2 (RPKM + 1)) per group (for either 11 developmental stages or 29 age groups) was used as a covariate conditioning on the average across all the labels. To control for multiple tests, a Bonferroni correction was used to control for 11 and 29 tests separately.

Cell-type enrichment

As previous studies had indicated the importance of cortical tissues to differences in SEP7,10, a gene-property analysis was also conducted to examine a broad array of brain-specific cell types. Enrichment of heritability was tested against 173 types of brain cells (24 broad categories of cell types), which were calculated following the method described in Skene et al.31. Briefly, brain cell-type expression data were drawn from single-cell RNA-Seq data from mouse brains. For each gene, a specificity value for each cell type was calculated by dividing the mean Unique Molecular Identifier (UMI) counts for the given cell type by the summed mean UMI counts across all cell types. MAGMA13 was used to calculate cell-type enrichments where specificity values were then divided into 40 equal-sized bins for each cell type for the MAGMA analysis. A linear model was fitted over the 40 specificity bins (with the least specific bin indexed as 1 and the most specific as 40). This was done by passing the bin values for each gene using the ‘--gene-covar onesided’ argument.

Univariate LDSC

Univariate LDSC regression was performed on the summary statistics from the GWAS on household income in order to quantify the degree to which population stratification may have influenced these results.

For the GWAS on household income, LD score regression was carried out by regressing the GWA test statistics (χ2) from each GWAS onto the LD score (the sum of squared correlations between the MAF count of a SNP with the MAF count of every other SNP) of each SNP. This regression allows for the estimation of heritability from the slope, and a means to detect residual confounders using the intercept.

LD scores and weights were downloaded from http://www.broadinstitute.org/~bulik/eur_ldscores/ for use with European populations. SNPs were included if they had an MAF of >0.01 and an imputation quality score of > 0.9. Following this, SNPs were retained if they were found in HapMap 3 with MAF > 0.05 in the 1000 Genomes EUR reference sample. Following this, indels and structural variants were removed along with strand-ambiguous variants. SNPs whose alleles did not match those in the 1000 Genomes were also removed. The presence of outliers can increase the standard error in LDSC regression, and so SNPs where χ2 > 80 were also removed.

Partitioned heritability

Partitioned heritability was performed using stratified LDSC regression74,75. Stratified LD scores were calculated from the European-ancestry samples in the 1000 Genomes project and only included the HapMap 3 SNPs with an MAF of >0.05. The model was constructed using 60 overlapping, functional categories. In addition, ten MAF bins and six continuous annotations were included to control for LD-related bias in the partitioned heritability analysis by modelling regional LD, as well as MAF. Correction for multiple testing was performed using a Bonferroni test on the 60 functional categories (α = 0.00083). The continuous annotations were also analyzed by examining the enrichment of each quintile for the six continuous categories of predicted allele age, background selection, recombination rate, nucleotide diversity, low levels of linkage disequilibrium in African populations and CpG content. Here, control for multiple testing was performed using a Bonferroni correction within each of the six annotations (α = 0.05/5 = 0.01).

Cell-type analysis was conducted using the method of Finucane et al.76. Here, four data sets were used and examined for enrichment of household income. The first data set (gene expression) contained gene expression data from across 205 tissue or cell types taken from the GTEx73 database and from Franke lab data set77,78 from Finucane et al.76. The second data set (chromatin) contained data on 489 tissue and cell types taken from Roadmap Epigenomics consortium67 and from EN-TEx, a subgroup of ENCODE76,79. Data pertaining to expression in 13 regions, the brain was taken from GTEx73 and gene expression specific to the neuron, the astrocytes and the oligodendrocytes were taken from mouse data from the work of Cahoy et al.80.

Multiple testing for the partitioning of the heritability by cell types was conducted using a Bonferroni correction across the 13 brain regions (α = 0.05/13 = 0.004) and across the three types of neurons (α = 0.05/3 = 0.017). For the gene expression and chromatin groupings, an FDR81 was applied to the 205 tests performed to look at enrichment using gene expression (α = 0.006) and to the 489 tests examining chromatin-based annotations (α = 0.003).

Mendelian randomization

The causal effects of intelligence (termed the exposure in an MR analysis) on income (termed the outcome in an MR analysis) were investigated using univariate MR analysis. Here, the total causal effect of intelligence on income was examined by combining summary GWAS test statistics for intelligence and for income using an inverse-variance-weighted (IVW) regression model82. This is equivalent to a weighted regression of the SNP-outcome coefficients on the SNP-exposure coefficients, with the intercept constrained to zero (i.e. assuming no unbalanced horizontal pleiotropy).

The results of the IVW regression model were compared with the results obtained using MR–Egger regression83. MR analyses, which use multiple SNPs, are more likely to include invalid SNPs with horizontally pleiotropic effects84. By not constraining the intercept to zero (as done using IVW regression) MR–Egger relaxes the assumption that the effects of genetic variants on the outcome act solely through the exposure (in this case intelligence). The intercept parameter of the MR–Egger regression indicates the average directional pleiotropic effects of the SNPs on the outcome. As such, the direct pleiotropic effect that the SNPs have on the outcome, independent of the exposure, can be quantified, where a nonzero intercept provides evidence for bias due to directional pleiotropy and a violation of the MR IVW estimator assumptions. Of note is that the MR–Egger regression estimates only remain consistent if the magnitude of the gene-exposure associations, across all variants, are independent of their horizontally pleiotropic effects on the outcome (i.e. the InSIDE assumption holds)83. In addition, power is almost always lower for MR–Egger and it requires variation in the size of effect of the SNPs on the exposure (i.e. if all SNPs have similar sized effects on the exposure, then MR–Egger will have very low power).

For use with MR, two independent groups (n = 95,521 for intelligence and n = 271,732 for income) were created, whereby the GWAS on income was rerun using only those participants whose data were not included in the interim release of the UK Biobank genotype data. A GWAS data set on intelligence was created by meta-analyzing publicly available data on intelligence with a GWAS (conducted for this study) using data from the INTERVAL study16,17 (Supplementary Methods) where 19 SNPs were identified as being genome-wide significant and independent. These 19 SNPs were used as instrumental variables for intelligence in the MR analysis.

Genetic correlations

Genetic correlations were derived using bivariate LDSC regression. A total of 27 GWAS data sets on health, anthropometric, psychiatric and metabolic traits were examined for a genetic correlation with income (Supplementary Table 16). Genetic correlations were also derived between household income with education and intelligence. There were three objectives to our analysis examining genetic correlations using household income. First, we sought to replicate the results of Hill et al.7, who found genetic correlations between household income and other variables in a smaller data subset from the UK Biobank sample used here. Second, SEP is multi-dimensional in nature: it is composed of multiple measures, each of which are correlated imperfectly with the others. Because of this, different measures of SEP may have genetic variance that is both unique to them, and differentiates them from the others in the way it associates with health. To examine this, we compare how genetic correlations with household income and 27 health, anthropometric, psychiatric, cognitive and metabolic traits differed compared with the genetic correlations derived using a different, individual-level measure of SEP, that is, educational attainment as measured by the number of years one has spent in education11. Third, Hill et al.7 also found that the phenotypes with the strongest genetic correlations with income are those that are cognitive (verbal numerical reasoning, childhood IQ and years of education) in nature7. The magnitude of these genetic correlations might indicate the phenotypes that occur as potential mediators between molecular genetic inheritance and household income.

In addition, intelligence is known to be genetically correlated with many physical and mental health traits18,21,85. The role that intelligence might play in accounting for some of the genetic links between household income and 27 health and well-being, anthropometric, mental health and metabolic traits was examined using genetic correlations. Here, the GWAS of income was conditioned on a GWAS on intelligence using multi-trait-based conditional and joint analysis (mtCOJO). mtCOJO is used to perform conditional GWAS whereby the genetic effects from one GWAS are controlled for in another GWAS. Importantly, the mtCOJO method avoids well-known issues of collider bias that can occur by including heritable covariates86. In the current study, the GWAS on income was conditioned on a GWAS on intelligence (and the intelligence GWAS was conditioned on the income GWAS) before the genetic correlations between income (and intelligence) and 27 variables mentioned above were reran.

Genetic prediction

The summary statistics from our GWAS of household income PGRSs were derived using PRSice-287 and the Generation Scotland:Scottish Family Health Study (GS:SFHS) cohort. The recruitment protocol and sample characteristics of GS:SFHS are described in full elsewhere88,89. In brief, 23,690 participants were recruited through their GP from across Scotland. Participants were all aged 18 years and over and were not ascertained based on the presence of any specific disease. Following the removal of individuals who preferred not to answer, income was assessed in GS:SFHS by 5-point scale (1 less than £10,000, 2 between £10,000 and £30,000, 3 between £30,000 and £50,000, 4 between £50,000 and £70,000 and 5 more than £70,000). Individuals who preferred not to answer were excluded from the analysis. Individuals who had taken part in UK Biobank were also removed from the GS:SFHS data set (n = 174). SNPs were included in the data if they had an MAF of ≥0.01 and Hardy–Weinberg P value >0.000001. Finally, one from every pair of related individuals were removed from the data set by creating a genetic relationship matrix using GCTA90 and removing individuals who are related at ≥0.025. This yielded a final sample size of 6680 participants who had genotype data and income data.

The participant’s level of income was then used as a predictor in a regression analysis with age, sex and 20 principal components included to control for population stratification. The standardized residuals from this model were then used as each participant’s income phenotype. PGRSs were created using the income phenotype derived using UK Biobank.

In each instance PRSice-2 was used to create five PGRSs corresponding to one of five P value cutoffs (P ≤ 0.01, P ≤ 0.05, P ≤ 0.1, P ≤ 0.5 and P ≤ 1) applied to the association statistics from the summary data. The polygenic risk scores for each threshold were then standardized and used in a regression model to predict the income phenotype in GS:SFHS.

Multi-trait analysis of GWAS