Gene-Covariate Interaction between Dysplastic Nevi and the CDKN2A Gene in American Melanoma-prone Families¹

CMM³is a potentially fatal form of skin cancer, the etiology of which is heterogeneous and complex. Approximately 10% of malignant melanomas develop in individuals with a familial predisposition (1)and often in association with clinically dysplastic or atypical nevi (2).

To date, two genes have been implicated in melanoma pathogenesis. The first, CDKN2A, encodes a low molecular weight protein, p16,that inhibits the activity of the cyclin D1-CDK4 complex (3). This complex phosphorylates the retinoblastoma protein, allowing the cell to progress through the G₁ cell cycle checkpoint. Thus, p16 acts as a tumor suppressor and negatively regulates cell growth by arresting cells at G₁. Germ-line CDKN2Amutations have been detected in 10–25% of melanoma-prone families that have been examined from North America, Europe, and Australia. In contrast, the second melanoma gene CDK4 acts as an oncogene (4), and germ-line mutations have been detected in only three melanoma-prone families worldwide (5, 6). Other genetic factors remain to be identified.

Although CDKN2A and CDK4 confer substantial risk for melanoma, other environmental and genetic factors are likely involved in the pathogenesis of familial melanoma. Epidemiological studies of melanoma suggest that exposure to sunlight is the major environmental risk factor, although the exposure response relationship appears complex with intermittent sun exposure likely to be more important for risk than total lifetime exposure (reviewed in Ref. 7). The major host factors associated with melanoma are increased numbers of melanocytic nevi, both clinically banal and atypical (dysplastic; Refs. 7 and 8). Previous examination of 13 melanoma-prone families with CDKN2A mutations, a subset of the families in the current study, showed results consistent with that seen in epidemiological studies (9), i.e., the same types of clinical and environmental risk factors were observed.

Specifically, examination of the relationship between clinical,environmental, and genetic factors and the risk of melanoma in these 13 families with CDKN2A mutations revealed that several variables, in addition to CDKN2A mutation status, appeared to enhance disease risk (9). The exposures that most consistently increased melanoma risk were TN, the presence of DN, and evidence of SI, as measured on the backs of shoulders (9, 10). To search for interactions between the CDKN2Agene and these clinicoepidemiological risk factors, we conducted combined segregation/linkage analysis using the regressive models (11, 12, 13) in 20 American melanoma-prone families.

The study subjects were drawn from 20 families in which there was a history of invasive melanoma in at least two first-degree relatives. The families were referred to the National Cancer Institute through health-care professionals or self-referrals. All of the families were Caucasian and resided in various regions of the United States. The families have been followed prospectively for 4–22 years.

Written informed consent was obtained prior to participation under an Institutional Review Board-approved protocol. All family members willing to participate in the study were clinically evaluated. Clinical evaluation of family members and spouses included complete skin examination, routine medical history, and phlebotomy to obtain lymphocytes. Variables recorded during the clinical examination included the type, distribution, and TN, and evidence for SI (defined as clinical assessment of absence/presence of solar keratosis, thinning of epithelium, idiopathic guttate hypomelanosis, wrinkling, and telangiectasia on backs of shoulders). For subjects who had multiple clinical evaluations, variables for the current study were taken from the most recent examination during which covariate data were recorded. All diagnoses of melanoma were confirmed using histological review of pathological material, pathology reports, or death certificates. Invasive and in situ cutaneous melanomas were classified as melanoma. Only confirmed melanomas were included in the analyses.

Table 1 shows the numbers of CMM cases, median age at CMM diagnosis, total family members, and specific alteration by family. The families ranged in size from 6 to 68 members with 2–12 CMM cases. The 20 families included 101 invasive and 10 in situ CMM cases and their 396 relatives (275 blood relatives and 121 spouses). Thirteen of the families had germ-line CDKN2A mutations that cosegregated with melanoma (9, 14). The mutations included one nonsense (Arg58Ter), one splice-donor-site (IVS2 + 1), five missense(Val126Asp, Gly101Trp, Arg87Pro, Asn71Ser, and Ala148Thr) mutations,one insertion (23ins24) and two deletions (234del14, 167del31). Melanoma susceptibility gene(s) have not yet been identified in the other seven families; these seven families were negative for both CDKN2A and CDK4 mutations. In addition, these seven families did not show statistically significant evidence for linkage to chromosome 9p21, the location of the CDKN2A gene.

Combined segregation/linkage analysis of CMM was conducted by use of the class D regressive logistic model (11) extended to take into account variable age at diagnosis of disease (15) and linked marker loci (12). The regressive models are constructed by specifying a regression relationship between each person’s phenotype (affected/unaffected with CMM) and a set of explanatory variables, including the person’s major genotype, the phenotype of older relatives (to take into account family dependence of unspecified origin, genetic, and/or environmental), and measured covariates. The cotransmission of CMM and the CDKN2A locus was analyzed as follows. When considering the cosegregation of a trait (Y) and a marker (M),the likelihood of a family can be written as:

where g_Y is the vector of underlying genotypes at the disease/trait locus, g_M, the vector of genotypes at the observed marker locus, and X is a vector of measured covariates. P(g_Y, g_M) is the joint probability of the genotypes at the disease and marker loci and P(Y| g_Y, X) is the penetrance function. Note that P(M |g_M) = 1 for a codominant marker. For individuals with no ancestors in the pedigree (founders and spouses), P(g_Y, g_M) is a function of the haplotype frequencies, p(d-1), p(d-2), p(D-1), and p(D-2),where d and D are the alleles at the disease locus and 1 and 2 are the marker alleles; the sum of these haplotype frequencies is constrained to be one. For individuals with ancestors in the pedigree, P(g_Y, g_M) is a function of the Mendelian probabilities and the recombination fraction θ between the two loci. Because the dominantly inherited CDKN2A gene is here assumed to be the gene causing CMM in our family sample, there is complete identity between the disease and marker loci, which leads to p(d-2) = p(D-1) = 0 and θ = 0. Thus, in this situation, when g_M =g_Y, the likelihood (A) becomes:

where the summation is over all possible genotypes(g_{Y_unobs}) in subjects untyped for the marker conditional on the observed genotypes(g_{Y_obs}) in their family members. Thus, only one allele frequency needs to be estimated, which will be designated CDKN2A− for the allele not carrying the CDKN2Amutation, whereas CDKN2A+ is the mutated allele. The likelihood was computed with the REGRESS computer program, which uses Eq. A for general situations. As stated above, Eq. A is equivalent to Eq. B when there is complete confounding between the putative disease gene and the marker, as is the situation here.

The penetrance function P(Y |g_Y, X) over nindividuals in the family is decomposed in a product of penetrance functions for each individual i: P(Y| g_Y, X) =∏ _i=1ⁿ(P_i|g_Yi, Y_ji, X_i), where g_Yi is the ith person’s genotype at the CDKN2A locus, Y_ji are the phenotypes of jantecedents of i, and X_i,the vector of covariates of i. Survival analysis concepts are introduced to take into account the variable age at diagnosis of CMM. The period of follow-up (taken here from 5 years of age because the risk of CMM is negligible before age 5), to age at diagnosis for affecteds, age at examination for unaffecteds (or affecteds with unknown age at diagnosis), or age at death for deceased subjects is partitioned into K mutually exclusive intervals. In each interval, we compute the hazard function, λ(k), the probability of being affected in the kth interval given not being affected before. The penetrance function is then derived fromλ(k). It is a density function at a given age of diagnosis included in the kth interval for affecteds, f(k) = λ(k)∏ _h=1^k−1 [1 −λ(h)] and a cumulative function for unaffecteds with age at examination in the kth interval, S(k) =∏ _h=1^k [1 −λ(h)]. For affected individuals with an unknown age at diagnosis and age at examination in the kth interval, the penetrance becomes F(k) = 1 −S(k). The penetrance function is one for individuals with unknown disease status. The hazard function λ(k) for the ith person in the kth interval is a logistic function λ(k) = exp[θ_i(k)]/(1 +exp[θ_i(k)]) whereθ _i(k), the logit of the hazard function, is:

where α_g is the genotype-specific baseline parameter, Γ_ji is a vector of regression coefficients on j antecedents’ phenotypes of the ith person, and β_g is a vector of genotype-specific regression coefficients for covariates. We assumed a dominant mode of inheritance for the CDKN2A locus. The class D model specifies four types of family dependences of the ith person on his/her antecedents: spouse(Γ_Si), father (Γ_Fi),mother (Γ_Mi), and preceding siblings(Γ_Ci). The covariates included age and the three clinicoepidemiological covariates (TN, DN, and SI). Age was the only time-dependent covariate, whereas the effects of the others were assumed to be constant over time. Age was treated as a continuous covariate in the model, and different functions of age were considered,polynomial in age and natural logarithm of age. The logarithm function was found to fit the data better and was subsequently used in all analyses. In addition, no analyses showed evidence for an age by CDKN2A mutation interaction; therefore, all models presented include only one regression coefficient for the logarithm of age.

Parameter estimation and tests of hypotheses were carried out using maximum-likelihood methods as implemented in the program REGRESS (12, 13), which incorporates the regressive approach in the ILINK program of the LINKAGE package (16) and uses the GEMINI optimization routine (17). The joint likelihood of the CMM phenotypes and CDKN2A genotypes was maximized under different models, always including an autosomal dominant major gene effect at the CDKN2A locus. Incorporation of familial dependencies into the models showed results similar to the analyses that ignored familial dependencies but with less power to differentiate between models because of the increased number of parameters (data not shown). Therefore, our results are presented without residual familial dependencies. However, additional analyses were conducted by introducing the parental affection status as a covariate in the model to allow for aggregation of disease in the seven families without cosegregating CDKN2A mutations (as explained in the“Discussion”). Gene-covariate interactions were tested by comparison of submodels in which the βgs were set equal (no interaction) versus models in which two βgs (under a dominant model) were estimated(interaction). Effects and interactions of the three risk factors (TN,DN, and SI) with the major gene were tested separately and jointly.

Two types of analyses were performed to deal with missing/unknown covariates (18). The first approach (MI method) created two dummy variables for each covariate: a missing-value indicator set equal to 1 for missing/unknown and 0 for known; and a second variable set equal to 1 for exposed and 0 for unexposed or unknown. Thus, tests for gene-covariate interactions required comparison of four models: 1,β(MI) and β(E) for variables missing value indicator (MI) and exposure (E), both independent of genotype; 2, βg(MI) dependent on genotype but β(E) independent of genotype; 3, βg(E) dependent on genotype but β(MI) not dependent; and 4, βg(MI) and βg(E), both dependent on genotypes. If a test of model 1 versus model 2 was rejected, i.e., a test of whether the effect of the missing covariate on disease risk was independent of genotype, the test of interaction for the exposure variable would compare model 2 versus model 4; otherwise, the test of interaction would compare model 1 versus model 3. The second approach to deal with missing covariates coded subjects with missing covariates as unknown for the disease status (the CC method). Therefore, for the CC method, tests for gene-covariate interaction required comparison of only β(E) and βg(E).

To search for gene-covariate interactions, we dichotomized the three previously identified risk factors (9) into referent/unexposed and risk/exposed categories. The baseline and at-risk categories for the covariates were defined as follows: TN, if age <16 years or age >50 years, <50 nevi versus ≥50 nevi, and if age 16–50 years, <100 nevi versus ≥100 nevi;DN, absent versus present; and SI, absent versuspresent. Table 1 presents a distribution of the covariates in each family.

Table 2 shows the results of the combined segregation/linkage analyses for the MI method. Overall, there was a significant improvement in the fit of the model when TN (model 2 versus model 1:χ ²₂ = 25.7, P < 0.001), DN (model 3 versus model 1:χ ²₂ = 52.3, P < 0.001), or both nevi variables (model 5 versus model 1:χ ²₄ = 65.9, P < 0.001) were included. The effect of either one of these variables was significant when the other one was included in the model. In contrast to DN and TN, inclusion of SI did not significantly improve the likelihood (model 4 versus model 1:χ ²₂ = 3.4, P =0.19). TN showed no evidence for gene-covariate interaction either alone (model 2a versus model 2:χ ²₁ = 0.2) or with DN in the model (model 5a versus model 5:χ ²₁ = 0.6). In contrast, there was significant evidence for interaction between DN and CDKN2A when DN was the only covariate in the model (model 3a versus model 3: χ²₁= 4.3, P = 0.04). We also found that the regression coefficient for the unknown/known dummy variable was independent of genotype (data not shown). The DN-CDKN2A interaction was no longer significant when TN was added to the model (model 5b versus model 5: χ²₁= 2.2). Interestingly, the estimate of the regression coefficient for DN (model 3a) was greater in subjects without CDKN2Amutations (β = 2.8) versus those with mutations(β = 1.8). The OR for DN in subjects without mutations was 16.4(95% CI, 6.1–46.1) versus OR, 6.0 (95% CI, 2.3–15.5)among subjects with CDKN2A mutations.

The CC method (Table 3) showed results similar to the MI method with respect to addition of covariates to the model. Inclusion of DN, TN, or both significantly improved the likelihood (model 3 versus model 1:χ ²₁= 47.6, P < 0.001; model 2 versus model 1:χ ²₁ = 21.1, P < 0.001; model 5 versus model 1:χ ²₂= 52.1, P < 0.001, respectively), whereas adding SI did not significantly improve the fit of the model(χ²₁ = 0.6, P = 0.44). TN revealed no evidence for interaction with CDKN2A. DN showed significant evidence for interaction with CDKN2A when it was the only covariate (model 3a versus model 3: χ²₁= 4.0, P = 0.046) or when TN was also in the model(model 5b versus model 5:χ ²₁ = 4.2, P =0.04). This latter result contrasted with the MI analysis, where a significant DN-CDKN2A interaction was observed only when DN was the sole covariate in the model. For the CC method, the estimates of the regression coefficients for DN were similar to what was observed in the MI analysis; also, β(DN) was significantly greater in subjects without CDKN2A mutations (β − (DN) = 3.0) versus those with mutations (β + (DN) = 1.2). The OR for DN was 20.1 (95% CI, 4.8–92.8) in subjects without CDKN2A mutations compared with OR, 3.3 (95% CI, 1.1–10.0)among subjects with mutations.

Restriction of the analysis to the 13 families with cosegregating CDKN2A mutations produced very similar results. The same evidence for inclusion of covariates (addition of DN to the model, P < 0.001; TN: P < 0.001; DN and TN: P < 0.001; SI: P = 0.40) as well as DN-CDKN2A interaction [DN-CDKN2A interaction with DN as the only covariate:χ ²₁ = 3.9, P < 0.05, β − (DN) = 10.0, β + (DN) = 1.5;DN-CDKN2A interaction with DN and TN as covariates:χ ²₁ = 4.2, P < 0.05, β − (DN) = 10.2, β + (DN) = 1.2]. The estimate of the regression coefficient for DN in subjects with CDKN2A mutations (β+) was similar to what was observed in the 20 families. However, in subjects without CDKN2Amutations (β−), this estimate, although higher than in the whole sample, should be regarded with caution because of the small numbers of CMM subjects without mutations.

Combined segregation/linkage analysis of CMM in 20 American melanoma-prone families showed significant evidence for inclusion of the covariates DN and TN in the regressive models. A significant gene-covariate interaction between DN and CDKN2A was observed when DN was the only covariate in the model (MI and CC methods) or when both DN and TN were in the model (CC method only).

The CDKN2A gene has been implicated in melanoma pathogenesis from studies of melanoma-prone families, cultured and noncultured melanoma tumors, and functional studies of specific mutations. The pattern of CDKN2A mutations observed in cultured and noncultured melanoma tumors (19, 20, 21, 22)indicates a possible role for UV radiation in the tumorigenic process. One hypothesis suggests that early loss of the p16 pathway may result in a net increase in melanocyte proliferation (23). This proliferation may then increase the probability that a melanocyte may accumulate additional damage during regular or UV-induced cell division, leading to development of a melanoma tumor. The finding of a significant association between an increased number of melanocytic nevi, the presence of DN, and the risk of melanoma are consistent with this proposed tumorigenic process.

Epidemiological studies of melanoma have demonstrated the importance of TN and DN as major risk factors for melanoma (reviewed in Ref. 7), yet few studies have been able to separate the effects of TN and DN. When TN and DN have been separated, the effect of DN appeared to be much greater than for TN (8). For this study, TN exposure was defined as an increased number of typical and atypical (dysplastic) nevi. DN exposure was defined by the presence of one or more unequivocal DN. Thus, there is overlap in the phenotypic definition of these two covariates. However, TN and DN were each found to influence separately the risk of melanoma, with DN showing an interaction with CDKN2A mutations.

Previous examinations of melanoma-prone families showed that DN did not appear to cosegregate with the CDKN2Amutations found in many families (14, 24, 25). The current study confirmed and quantified this relationship and showed that DN was a risk factor for melanoma, separate from CDKN2A; the risk for CMM associated with DN was greater in subjects without CDKN2A mutations compared with subjects with mutations. That is, among CDKN2A mutation carriers, the presence of DN multiplied the hazard function by 3 compared with mutation carriers without DN. In contrast, among noncarriers, the presence of DN multiplied the hazard function by 20 versus subjects without DN. In addition, the estimate of the hazard function among noncarriers with DN was similar to the estimate of the hazard function among CDKN2A mutation carriers without DN. Restriction of the analysis to the 13 families with cosegregating CDKN2Amutations showed a very similar pattern of results. In addition,because seven families did not have CDKN2A mutations, we performed additional analyses that allowed for aggregation of disease separate from CDKN2A. That is, we added covariates for the seven families based on whether a parent had CMM, was unaffected, or phenotype unknown. The “parental” covariates for the 13 CDKN2A families were coded as unknown. This analysis yielded results that were very similar to the CC analysis. We observed the same evidence for inclusion of DN and TN in the model as well as DN-CDKN2A interaction with comparable estimates of the regression coefficients (data not shown). Thus, although DN was a major risk factor for CMM in subjects with CDKN2A mutations, it appeared to be a stronger risk factor in subjects without CDKN2A mutations. This relative strength of effect for DN may change when the susceptibility gene(s) in the families without CDKN2A mutations are identified.

One potential problem in genetic epidemiological studies of families is that covariate data may be missing because not all family members are available or willing to participate in a research study. To deal with missing or unknown covariates, we conducted two analyses. The basic idea of the MI method is to obtain regression information from those family members for whom such information is available. In the CC method, family members with missing covariates were coded as unknown for the affection status. We observed similar results for the MI and CC method, although the CC method showed additional evidence for the DN-CDKN2A interaction. One difference between the two methods is that the MI method requires estimation of twice as many parameters for covariates, one covariate for known/unknown and one covariate for exposed/unexposed, as the CC method. This requirement,plus the observation that the DN-CDKN2A interaction was no longer significant after TN was added to the model in the MI method,raises a concern about adequate power for detecting gene-covariate interactions. It also raises concerns about possible biases in the methods. For example, the CC method may produce biased estimates of the parameters if the data are not missing completely at random. Missing completely at random is equivalent to assuming that, for each variable,the observed values effectively constitute a simple random sample of the values for all study subjects (26). For DN, similar percentages of missing values were observed in melanoma patients (27%)and those who did not have melanoma (28%). The MI method has also been shown to yield biased results in some situations and particularly if the covariate of interest is an important confounder of the study effect (26). Although the problems associated with missing data have been extensively evaluated for logistic regression (26), the problems have been rarely evaluated for family studies, where the issues may be more complicated and difficult to resolve. The use of more complex methods to allow for missing data,such as Gibbs sampling, have been proposed (27, 28, 29) but these methods have not yet been studied extensively. These issues will be further investigated.

A correction for ascertainment for these families was not applied to the likelihood calculation because this requires the mode of ascertainment to be known precisely. When it is unknown, as is the situation here, the MOD score strategy, which considers the likelihood of disease status and marker conditionally on the disease status of all family members, would be appropriate and lead to unbiased parameter estimates when there is no association between disease and marker (30, 31). In the situation of a disease-marker association, the MOD score should include linkage disequilibrium between disease and marker alleles not only in the numerator but also in the denominator. However, such conditional likelihood may remove most of the information, especially when complete linkage disequilibrium is specified, and may not be efficient. Although use of joint likelihood may affect genetic parameter estimates when the gene is unknown, this is less true when the gene is known in most subjects,as is the situation here. Although some error in inferring the gene carrier status in those subjects with unknown CDKN2A status adds noise, it is expected to decrease power for detecting gene-covariate interactions rather than causing a false-positive result.

The 20 families in the current study were ascertained because they had at least two living first-degree relatives with melanoma. The average number of melanoma cases in the families was 5.5 (± 2.9). Although the families are not representative of all melanoma in the general population, they likely represent a rarer subset of melanoma-prone families with a strong genetic risk. These are the types of families that have been used to map melanoma and other disease genes(e.g., BRCA1/2). These types of families have yielded substantial information about the genetics of cancer and interactions with other genetic and environmental covariates. Thus, the results from these melanoma-prone families may be generalizable to that percentage of melanoma-prone families with strong genetic risks.

The etiology of melanoma is heterogeneous and complex. To date, two genes have been implicated in the inherited form(s) of melanoma. The major known gene CDKN2A confers substantial risk for melanoma. However, other genetic and environmental factors likely also contribute to disease expression. Interactions between a melanoma susceptibility gene and risk factors were suggested from multiple-case pedigrees from Utah where putative gene carriers, inferred from linkage analysis with 9p21 markers, had a higher number of nevi than noncarriers, and among gene carriers sun exposure was greater in subjects with than without melanoma (32). A significant interaction between a putative melanoma gene, detected by segregation analysis, and propensity to sunburn was reported in a large series of 295 melanoma French families ascertained through one melanoma proband (33). The current study of 20 American melanoma-prone families is the first to detect evidence for a significant interaction of a covariate with a known melanoma gene. The interaction between CDKN2A mutations and DN illustrate the complex etiology of melanoma. These results need to be confirmed in larger samples and in families from other geographic regions. In addition, the best approaches to deal with missing covariate data and the problems of adequate power to detect gene-environment interactions using the regressive models needs to be further assessed.

We are indebted to the participating families, whose generosity and cooperation have made this study possible. We acknowledge the Biomedical Supercomputing Center of the National Cancer Institute-Frederick Cancer Research and Development Center for allocation of computing time and staff support. We also acknowledge the contributions to this work that were made by Mary Fraser, Laura Fontaine, and Deborah Zametkin.