Heritability and GWAS Analyses of Acne in Australian Adolescent Twins

Abstract

Acne vulgaris is a skin disease with a multifactorial and complex pathology. While several twin studies have estimated that acne has a heritability of up to 80%, the genomic elements responsible for the origin and pathology of acne are still undiscovered. Here we performed a twin-based structural equation model, using available data on acne severity for an Australian sample of 4,491 twins and their siblings aged from 10 to 24. This study extends by a factor of 3 an earlier analysis of the genetic factors of acne. Acne severity was rated by nurses on a 4-point scale (1 = absent to 4 = severe) on up to three body sites (face, back, chest) and on up to three occasions (age 12, 14, and 16). The phenotype that we analyzed was the most severe rating at any site or age. The polychoric correlation for monozygotic twins was higher (r_MZ = 0.86, 95% CI [0.81, 0.90]) than for dizygotic twins (r_DZ = 0.42, 95% CI [0.35, 0.47]). A model that includes additive genetic effects and unique environmental effects was the most parsimonious model to explain the genetic variance of acne severity, and the estimated heritability was 0.85 (95% CI [0.82, 0.87]). We then conducted a genome-wide analysis including an additional 271 siblings — for a total of 4,762 individuals. A genome-wide association study (GWAS) scan did not detect loci associated with the severity of acne at the threshold of 5E-08 but suggestive association was found for three SNPs: rs10515088 locus 5q13.1 (p = 3.9E-07), rs12738078 locus 1p35.5 (p = 6.7E-07), and rs117943429 locus 18q21.2 (p = 9.1E-07). The 5q13.1 locus is close to PIK3R1, a gene that has a potential regulatory effect on sebocyte differentiation.

Acne vulgaris (acne) is an inflammatory and chronic skin disease that primarily affects adolescents and young adults (ICD-10; WHO, 2015). Epidemiological data indicates that the distribution of acne is almost universal, and the incidence in young adults is on the rise (Lynn et al., 2016; Tan & Bhate, 2015). The disease is caused by an obstruction of the hair follicle, involving hyperactivity of the sebaceous glands, alterations of the composition of the skin bacteria, inflammation, and keratosis. It can develop in any area of the body, but lesions are generally present on the face, neck, chest, and upper back (Williams et al., 2012). The severity of the inflammation ranges from mildly inflamed lesions (comodones) to severely inflamed lesions (pustules, nodules, and papules). In extreme cases the lesions develop into nodulocystic acne. The scarring generated in extreme cases is associated with psychological distress and social inability (Barnes et al., 2012; Dunn et al., 2011; Halvorsen et al., 2011).

Despite the pervasiveness of the disease, little is known about the causes or origin of acne. Factors, such as diet, seasonal changes, hormonal changes, and ethnicity have been suggested as risk factors for the occurrence and severity of acne. However, there are no conclusive findings that show any of these factors as causative of acne (Bataille et al., 2002; Cho et al., 2014; Magin et al., 2005). Multiple twin studies (Bataille et al., 2012; Cho et al., 2014; Friedman, 1984; Walton et al., 1988) have reported that genetic factors influence the development of acne, including an earlier analysis of our cohort (Evans et al., 2005). There are also indications that the severity of the disease and time of onset are predisposed by genetic factors (Evans et al., 2005; Goulden et al., 1999). Several genes within the pathways of the cytochromes P450 (CYPs) family have been suggested as candidates in the pathogenesis of acne (Paraskevaidis et al., 1998). Insulin, insulin receptors (Arora et al., 2011; Melnik, 2010), insulin like growth-factor-I (Melnik & Schmitz, 2009; Tasli et al., 2013), and androgen receptors (Sawaya & Shalita, 1998; Yang et al., 2009) have also been reported as having an effect on the development of acne. Recently, Lichtenberger et al. (2017) reviewed 21 studies of candidate genes. A common observation around the reviewed studies was that further analyses are required to corroborate or refute the findings from candidate genes studies, which were often inconclusive due to the small number of individuals in the studies, lack of statistical power, and/or failure to replicate.

While GWAS has proved to be a powerful tool to detect genomic variants influencing complex traits, few GWAS of acne have been conducted. A GWAS study and meta-analysis in a Han Chinese cohort (He et al., 2014) reported suggestive association at two loci (1q24.2 and 11p11) located within the regions of DDB2 and SELL, genes that are implicated in the process of androgen metabolism, inflammation, and scarification. A follow-up of this study (Wang et al., 2015) confirmed the potential involvement of these loci in the pathogenesis of acne. Another GWAS study in European Americans (Zhang et al., 2014) did not detect SNPs with genome-wide association. A SNP at the proximity of the gene MYC (locus 8q24) had the highest association (1.7E-06), but the study lacked statistical power as it had only 81 cases and 847 controls and failed to replicate. A GWAS in a UK cohort (Navarini et al., 2014) reported genome-wide association at 1q41, 5q11.2, and 11q13. These loci are within the regions of genes that have functions related to skin homeostasis and are also associated with the TGFβ pathway, a growth factor with a putative role in the pathogenesis of acne. The UK study did not detect association at 8q24 or any of the loci reported for the Han Chinese cohort. Overall, the findings of these studies suggest that acne predisposition may have a complex genetic architecture of small genetic effects.

Here we report a study of heritability of acne for an Australian cohort of adolescent twins, extending a previous analysis by Evans et al. (2005) on the heritability of acne severity. The phenotype in our study differs from the Evans study as here we consider the highest score of acne among three longitudinal measures; we also include data from siblings, which gives more statistical power to the analysis. We apply twin modeling to re-estimate the heritability of acne using a sample three times larger (1,906 twin pairs plus 671 siblings) than that used by Evans et al. (2005), and then we conducted a GWAS of these data.

Materials and Methods

Participants

Participants were part of the Brisbane Longitudinal Twin study (BLTS) conducted at the QIMR Berghofer Medical Research Institute. Since 1992, the BLTS study has collected data from twins, their non-twin (singleton) siblings, and their parents. Families were first recruited into this study when the twins were 12-year-old, with additional follow-up data collections at ages 14, 16, 19, and 25. Within the BLTS study, data is collected on a wide range of biological, psychological, and social traits, as well as environmental exposures. Further details of the recruitment process, data collection, and determination of zygosity can be found in Wright and Martin (2004). All participants in the BLTS gave informed written consent; ethical approval was obtained from the Human Research Ethics Committee.

The present study used longitudinal data on the severity of acne for 3,817 twins and 674 siblings; the twin data includes 752 monozygotic (MZ) complete pairs, 1,154 dizygotic (DZ) complete pairs, and five unpaired individuals; 2,342 (52% of the total sample) were females and 2,149 (48% of the total sample) were males. Totals of pairs by zygosity and acne scores are presented in Table 1.

TABLE 1 Description of the Sample: Distribution of Acne Scores for Each Zygosity Group and Twin Polychoric Correlations of the Liability to Acne Severity

Note: The score for acne severity among twins and siblings corresponds to the highest rating for acne severity by site and location as explained in the methods. The table presents the number of individuals and complete pairs per zygosity group (MZM = Monozygotic males; MZF = monozygotic females; DZM = dizygotic males; DZF = dizygotic females; DZOS = dizygotic opposite-sex). r = polychoric correlations of the liability to acne severity for co-twin in each zygosity group and for twins-sibling. Polychoric correlations were estimated by SEM with full information maximum likelihood estimation and are standardized for sex, age, age², age × sex, and age² × sex. CI = 95% confidence interval of the polychoric correlation.

Measures

The severity of acne was rated by a nurse using a 4-point scale (1 = absent; 2 = mild; 3 = moderate; 4 = severe) on back, face, and chest. The measures were taken in a longitudinal manner at ages 12 and 14 years and from the face only at age 16 years. Details of the rating and validation of the scale are given in Evans et al. (2005). For this analysis we chose the score that corresponds to the most severe across all sites at all ages, to minimize missing data and to allow for individual differences in the time onset of acne, since boys and girls typically develop acne at different ages. Figure 1 presents the age distribution of the most severe acne score among females and males; the mean age of the twins was 14.6 (range 12–23, SD 1.4) and the mean age of the siblings was 14.7 (range 10–23, SD 2.4).

Note: Individuals are grouped by age in three categories (9–12, 13–15, and 16 or more years old). Categories of acne severity: 1 = absent; 2 = mild; 3 = moderate; 4 = severe, correspond to the most severe rating that each individual has across up to three body sites (face, chest, back) and up to three occasions (12, 14, and 16 years, age of the twin pair). The percentages of individuals within each category are shown in the text boxes.

FIGURE 1 Summary of the number of females and males who were scored for acne severity.

Model-Fitting Analysis

We conducted structural equation modeling (SEM) with full information maximum likelihood estimation of parameters, which makes use of all available information. As acne was scored as ordered categories, we fitted a liability-threshold model (Rijsdijk & Sham, 2002), which assumes that there is an underlying liability toward acne that is normally distributed (with mean zero and variance one). The categories can be considered as a series of thresholds that divide the liability distribution into discrete classes.

We began with a fully saturated model in which thresholds and correlations were allowed to differ by sex, zygosity, and whether it was a twin or a twin or sibling. Age and sex are associated with the onset of acne; so these fixed effects were removed from the thresholds by regressing sex, age, age², age by sex, and age² by sex on the thresholds. Thus, this is under the expectation that the beta effects are comparable between zygosity groups and siblings.

The saturated model consisted of 36 parameters. Within same sex pairs we equated thresholds for first- and second-born twins so there were three thresholds for each twin group (MZF, MZM, DZF, DZM, OS-female, OS-male) and three thresholds for each sibling group (male siblings and female siblings), for a total of 24 thresholds, four betas (age, age², age by sex, and age² by sex) and eight correlations (MZ-females, MZ-males, DZ-females, DZ-males, DZ-OS pairs, female twins-sibling, male twins-sibling, and OS-twins-sibling). A parameters reduction was conducted by testing the significance of the parameters thorough several nested models; we tested for differences in thresholds within twin pairs, zygosity groups, and twins and siblings.

Model fit was assessed using maximum likelihood-ratio test, which is asymptotically distributed as a chi-squared (χ ²) with degrees of freedom equal to the difference in the number of parameters between models. A non-significant p value (p > .05) indicates that the model with fewer parameters can be retained without a significant loss of fit. We also evaluated Akaike's information criterion (AIC), calculated as Δχ ²–2Δdf (Akaike, 1987), which combines parsimony and explanatory power by penalizing improving fit for the addition of parameters.

Genetic Analysis

As explained elsewhere (Neale & Cardon, 2013), the classical twin method makes use of the differences in genetic relatedness between MZ and DZ twins to estimate the proportion of variance due to additive genetic (A), unique environmental (E), and common environmental (C) or non-additive genetic (D) effects. If twice the DZ correlation is greater than the MZ correlation (2r _DZ > r _MZ), this is indicative of shared environmental effects so an ACE model is used. On the other hand when the MZ correlation is larger than twice the DZ correlation (r _MZ > 2r _DZ) this indicates non-additive genetic effects and an ADE model is used.

For the liability of acne severity r _MZ > 2r _DZ (Table 1) an ADE model was used, and the sub-models AE and E were fitted to test whether a model with fewer parameters could satisfactorily explain the variance in liability to acne severity.

Data handling and descriptive statistics were performed using SPSS (IBM Corp. Version 22.0. Armonk, NY: IBM Corp). Polychoric correlations and SEM analyses were estimated using the R package OpenMx version 2.6.7 (Boker et al., 2011; Boker et al., 2016).

Genome-Wide and Gene-Based Association Analysis

The individuals in this study were genotyped as part of a larger GWAS study. Briefly, genotyping was performed on HumanCoreExome-12v1-0_C or IlluminaHuman610W-Quad bead chip. Details of the genotyping procedures and quality control are given elsewhere (Medland et al., 2009). Genotypes were imputed to phase 3 version 5 of the 1000 Genomes Build37 (hg19) (1000 Genomes Project Consortium et al., 2015). We performed GWAS using RAREMETALWORKER (Feng et al., 2014) to explicitly correct for relatedness.

We used the summary data from the GWAS to conduct a gene-based association analysis using VEGAS2 (Mishra et al., 2015). This analysis aims to determine whether genes harbored an excess of variants with small p values. The test uses the summary statistics from a GWAS analysis to determine the evidence of association at the level of the gene, taking into account the SNP density, gene size and LD between SNPs. These analyses can be more powerful than individual SNP association (Liu et al., 2010). In these analyses, we used a MAF >0.01 and SNP coordinates based on Build37 (hg19). Our test included 20,944 genes and we used a p value of 2.3E-06 to declare significance.

Genotyping was available for 4,762 individuals (including 271 second siblings that were not included in the twin modeling). Following imputation quality control, data for ~7.5 million SNPs was available for the GWAS.

Results

Twin Modeling and Variance Components for Acne

The majority (63%) of the individuals assessed in this analysis had mild or moderate acne, with only 7% of the cohort showing severe acne. Thresholds could be equated within twin pairs, across zygosity groups, and for twins and siblings, without any loss of fit (p > .05).

Polychoric co-twin correlations (r _MZ =0.86, 95% CI [0.81, 0.90]; r _DZ = 0.42, 95% CI [0.35, 0.47]), and twin-sibling correlations (r _sibling = 0.34, 95% CI [0.26, 0.52]), indicates high heritability and additive and/or dominant effect (Table 1).

There was no significant change of fit (p = .122) when D was removed from the ADE model. However, the fit changed significantly (p < .001) when both A and D where removed (E model vs. AE model), which indicates that there are highly significant genetic contributions to the liability to acne. Estimates from the reduced AE model indicate that 85% (95% CI [0.82, 0.87]) of the variance was explained by additive genetic components and 15% (95% CI [0.11, 0.16]) of the variation was explained by unique environmental components (Table 2).

TABLE 2 Decomposition of the Genetic Variance of Acne Severity

Note: Results from variance components analysis for liability to acne severity; the genetic variance was decomposed into additive (A), non-additive (D) and unique environmental (E) factors. The AE model (in bold) was the most parsimonious model to explain the liability to acne severity. Confidence intervals for the parameter estimates are in brackets. −2LL = Minus twice the log-likelihood of the model. df = degrees of freedom.

GWAS and Genetic Association

Figure 2 shows the Manhattan plot of the associated p values and the quantile–quantile plot of the observed and expected –log₁₀ (p value). There was no evidence of confounding effects or inflation of the data (λ = 1.02). In Table 3 we provide a list of the top 10 SNPs from the GWAS analysis. No genome-wide significant (5E-08) variants were identified. But, suggestive associations were found for three SNPs: rs10515088 locus 5q13.1 (p = 3.9E-07), rs12738078 locus 1p35.5 (p = 6.7E-07), and rs117943429 locus 18q21.2 (p = 9.1E-07).

Note: (A) Manhattan plot of the genome wide association analysis for acne. Showing the best associated SNPs (p <E −06) and nearest gene. The vertical axis shows the –log₁₀ p values across the genome, and the horizontal axis shows the chromosome numbers. (B) Quantile–quantile plot of the GWAS analysis showing the genomic inflation factor (λ) and the distribution of -log₁₀ p values of the observed association for acne against the expected distribution.

FIGURE 2 Plots for the GWAS analysis of acne severity.

TABLE 3 Top 10 SNPs With the Highest Association with Acne

Note: Top 10 SNPs (with the lowest p value) from GWAS analysis of acne severity. Reference SNP identification, p value for the significance of the association, chromosome (in bold), followed by the position of the SNP in the human genome reference (hg19) and the nearest gene or genes to the SNP.

In order to search for potentially causative SNPs or other variants in high LD with these three SNPs, we made regional association plots using LocusZoom (Pruim et al., 2010); see Figure 3. The pairwise search found a group of eight SNPs that were in moderate LD (R ² > 0.6 < 0.8) with the locus rs12738078 (Figure 3A), but there were no other SNPs in high LD levels within ~1.5Mb of rs10515088 (Figure 3B), and rs117943429 (Figure 3C). All of these variants were intergenic.

Note: (A) rs12738078 in chromosome 1, (B) rs10515088 in chromosome 5, (C) rs117943429 in chromosome 18. Showing location, recombination rates, LD with possible associated regions and near genes.

FIGURE 3 Regional association plots for the top three SNPs associated with acne severity.

Figure 4 shows the Manhattan plot for the gene-association analysis; no gene reached genome-wide significance. Three genes were located in chromosome 17 (DHRS11, KAT7, PIGW) and one gene in chromosome 9 (ASS1). Genes DHRS11 and PIGW were located at the same locus and it is not possible to distinguish which of the two is tagged by the SNPs signals. The strongest evidence for association was found for the gene OR10J5 (p = 4.4E-05) an olfactory receptor located in chromosome five. Gene positions and p values are presented in Table 4. Regional association plots for the top five genes identified by the gene-base analysis are presented in Figure 5.

Note: The analysis input was a GWAS for acne severity with p values from ~7.5 million autosomal SNPs. Each circle in the plot indicates the start (base-pair) position of a gene, Build37 (hg19) was used as reference. Locations for the top five associate genes are shown, DHRS11 and PIGW are located in the same locus. No gene reached genome-wide significance which was –log10 (p value) >2.3E-06.

FIGURE 4 Manhattan plot for gene-association analysis from VEGAS2.

TABLE 4 Gene Association Analysis for Acne Severity Using VEGAS2

Note: Presenting the top five genes from the VEGAS2 analysis, chromosome (in bold), number of SNPs associated with each gene, p-value of the gene association analysis, the SNP with the highest association with each gene and p value of the SNP in the GWAS analysis for acne. * DHRS11and PIGW are located in the same locus and cannot be distinguished.

Note: Presenting the gene and associated SNP for OR10J5, in chromosome 1 (A), ASS1, PIGW, and DHRS11 in chromosome 17 (B and C) and KAT7 in chromosome 9 (D).

FIGURE 5 Regional plots for the top five genes identified by the gene-association analysis with VEGAS2.

Discussion

In this study, we estimated the heritability of acne severity in an unselected sample of Australian teenage twins and their siblings and searched for molecular variants that contribute to this genetic variance. Although our data did not include clinical cases and the majority of individuals had mild to moderate levels of acne, our analysis showed a high heritability for acne severity (h ² = 0.85, 95% CI [0.82, 0.87]) that was best explained by additive genetic effects and non-shared environments. These results are in concordance with previous studies of heritability of acne by us (Evans et al., 2005) and others (Bataille et al., 2002; Liddell, 1976; Szabó et al., 2011; Walton et al., 1988). Particularly, our results agree with the study by Evans et al. (2005), which used an earlier subset of the QIMR cohort analyzed here and found the heritability for acne to up to 97% (95% CI [0.91, 1.0]) at the site on the back.

Our study found suggestive association (p < 1E-07) for three SNPs: rs10515088, rs117943429, rs12738078, none of which have previously been reported as associated with acne. Interestingly, the SNP rs10515088 is in close proximity to the gene PIK3R1 (Figure 3B), which encodes for a phosphatidylinositol 3-kinase protein member of the PI3-Ks lipid kinases family. PI3K is part of the metabolic pathway of insulin, and a recent study by Ju et al. (2017) has identified that it plays a role in the interplay between androgens, insulin, insulin-like growth factor, and acne. In acne patients, the stimulation of PI3K was associated with a significantly reduced proliferation but increased differentiation of sebocytes. We suggest that the 5q13.1 locus may be a genomic region of interest for acne given the proximity with the PIK3R1.

While our analysis did not find genome-wide associated SNPs, the sample size of our study is relatively low, limiting the power of our study.

Our results support the hypothesis that acne is a trait influenced by the action of multiple loci with small effect. At least three metabolic pathways have shown to be involved in acne susceptibility: TGFb pathways, PI3K, and DDB, and it is likely that there are others that will be identified in the future. Detecting genomic regions associated with acne requires a larger sample size. As such, a meta-analysis that summarizes the data from different cohorts could be a powerful analysis to find more variants with more robust associations with acne.