signature=a702a44d1f46acac87955d98d4ef1534,Complex signatures of selection for the melanogenic loci

Cells and cell cultures

All melanocytes used were obtained from human skin. Melanocyte line M0202 was obtained from a lightly pigmented donor by one of us (NS) in the lab. Melanocytes M1 ("Caucasian", 1.5 years old, male donor) and M4 ("Negroid", 4 years old, male donor) were purchased from Promocell (Heidelberg, Germany). The 19-HEM cell line (lightly pigmented) was purchased from Gentaur (Belgium), and NHEM cells (two lightly pigmented and three darkly pigmented) were purchased from Cascade Biologics (Nottinghamshire, UK).

Cell cultures were maintained in an incubator at 37°C and 5% CO2. Melanocyte culture conditions were as indicated by the suppliers. Briefly: a) M0202 cells were initially grown in FETI medium: H-10, 2% FCS, 1% Ultroser G (BioSepra, Ciphergen, France), 4 ng/ml bFGF, 2 ng/ml endothelin-1, 5.3 nm TPA, 0.05 mM IBMX, but were later grown in H-10 supplemented with HEPES 6 mM, 5% FBS and MelanoMax (Gentaur), which presumably contains TPA, CT and contains BPE at 40 μg/ml (C. Stefanidis, Gentaur; personal communication); b) M1 and M4 melanocytes were cultured in melanocyte growth medium M2 (Promocell). M2 is a serum-free medium, without PMA (TPA) or other tumor promoting or toxic agents, consisting of a basal medium plus a Supplement Mix. The detailed formulation has not been disclosed by the supplier, but it seems that no BPE is present in this growth medium (Ute Liegibel, PromoCell, personal communication). c) 19-HEM melanocytes were grown in H-10 supplemented with HEPES 6 mM, 5% FBS and MelanoMax (Gentaur); d) NHEM melanocytes culture medium was Cascade Biologics Medium 254 supplemented with 1% HMGS (containing BPE, FBS, bovine insuline, bovine transferrin, bFGF, hydrocortisone, heparin and PMA) (Cascade Biologics).

The culture medium was changed every two days until the culture was approximately 80–90% confluent, and everyday thereafter.

M0202, HEM and NHEM melanocytes were passaged (split 1:3) routinely every 10 or 11 days, or harvested when they had reached confluency. The growth rate for M1 and M4 melanocytes was slow, and passaging (1:2) or harvesting at confluency was performed routinely after 2 to 3 weeks. All melanocytes used were from passage 5 to 15, and all were from normal non-transformed primary cell culture isolates.

Irradiation

Subconfluence cultures were irradiated once a day for 5 consecutive days with UVA+B (50 mJ/cm2:25 mJ/cm2) light in an ICH2 photoreactor (LuzChem, Canada) at 37°C. These doses assumed an absorbance of the plastic flask of approximately 5% for UVA and 11% for UVB (flasks were opaque to UVC), estimated from spectrophotometer absorbance readings at 255 nm, 305 nm and 360 nm of cuvette-size splinters of the flasks. By trials with murine melanoma cells (B16F10), this dose was shown to have no effect on cell viability. In order to prevent the generation of toxic metabolites, the culture medium was replaced by PBS with magnesium and calcium immediately before irradiation. After irradiation, PBS was replaced again by the culture medium. Irradiation control cultures were subject to the same procedure, except that they were covered by aluminum foil during irradiation. Cultures were harvested 24 h after the last irradiation dose.

Microarray gene expression

Immediately after harvesting, cells were resuspended in lysis buffer. Total RNA was obtained following the supplier's protocol (Ambion's total RNA extraction kit), including DNAse treatment. cDNA was synthesized from 2 μgr of total RNA using the Affymetrix One-Cycle cDNA Synthesis Kit and following the Affymetrix Expression Analysis Technical Manual. From this cDNA, cRNA was synthesized using the Affymetrix IVT Labeling Kit, which was then purified with the Affymetrix GeneChip Sample Cleanup Module. The purified cRNA (15 μgr) was fragmented and hybridized to Affymetrix U133A 2.0 arrays using standard Affymetrix protocols. A total of 18 microarrays from 9 irradiated cell lines (5 from lightly pigmented donors and 4 from darkly pigmented donors) plus their corresponding unirradiated control cultures were analyzed.

Microarray Data Analysis

Raw data were log2 transformed and quantile normalized using DNAMR v1.1 [46] for R (2.4.1). The Pattern Discovery option of SAM software v3.0 [47] was used to analyze the normalized data. Gene Ontology analysis was done using FatiGO [48].

Quantitative PCR

To demonstrate that the media composition, in particular the presence/absence of BPE, affects the expression levels of TYR, TYRP1 and DCT we purchased a new melanocyte cell line from Cascade Biologics (lightly pigmented) and grew this cell line in Cascade growth medium 254 supplemented with 1% HGMS (Cascade + for short). We sub-cultured and propagated the cells for approximately 3 weeks (three passages). Three days after the third passage we generated 4 subcultures of the same cell line. These 4 subcultures were from now on grown in 6 different media: Medium 1: Cascade +; Medium 2: Cascade medium, PromoCell Supplement (no BPE); Medium 3: PromoCell growth medium plus PromoCell Supplement. Medium 4: PromoCell growth medium plus PromoCell Supplement, plus 0.2% (v/v) of a 13 mg/ml solution of BPE. Cells were grown in these media for four days to allow some adaptation to the new media. Media were refreshed every two days.

After this time, we extracted total RNA from each subculture (Ambion). cDNA was synthesized using the Invitrogen SuperScript First-Strand synthesis system for RT-PCR kit, and then, we quantified the expression of TYR, TYRP1 and DCT for each of these subcultures using the BIO-RAD iQ SYBR green Supermix system in combination with a BIO-RAD iCycler machine. For mRNA quantification the following primers were used: TYR: 5'-AGAATGCTCTGGCTGTTTTG-3' and 5'-TCCATCAGGTTCTTAGAGGAGACAC-3'. For TYRP1: 5'-CATGCAGGAAATGTTGCAAGAG-3' and 5'-AGTTTGGGCTTATTAGAGTGGAATC-3'; For DCT: 5'-TATTAGGACCAGGACGCCCC-3' and 5'-TGGTACCGGTGCCAGGTAAC-3'. For normalization GADPH was used; primer: 5'-CCTGTTCGACAGTCAGCC-3' and 5'-CGACCAAATCCGTTGACTCC-3'. In all cases, annealing temperatures were fixed at 56°C. For DCT, Mg2+ concentration was increased in 1 mM above the standard reaction conditions.

DNA samples and Resequencing

We have resequenced the following in 116 human chromosomes: a) 4.1 kb of 5' DCT, including 187 bp of the first intron, the first CDS plus the 5'-UTR and 3,204 bp of the upstream flanking sequence; b) approximately 4 kb of 5' TYR, including 17 bp of the first intron, the first CDS plus the 5'-UTR and 2,773 bp of the upstream flanking sequence; and c) approximately 5 kb of 5' TYRP1, including 433 bp of the first intron, 47 bp of the second intron, the first CDS plus the 5'-UTR and 4,200 bp of the upstream flanking sequence. Sample individuals come from diverse geographical origins and include: 20 chromosomes from Biaka and Mbuti Pygmies (DNA purchased from the European Collection of Cell Cultures, ECACC), 20 chromosomes from Senegalese individuals resident in Spain, 42 European (N. Spain) chromosomes (including 20 chromosomes from melanoma patients), 20 Asian chromosomes including Chinese samples from Coriell Cell Repositories and Chinese residents in Spain, and 14 chromosomes from Australian Aborigines (DNA purchased from the ECACC).

DNA was PCR amplified in overlapping ~1 kb long segments and these were resequenced using BigDye 3.1 chemistry and ABI PRISM 3730 and ABI 310 DNA analyzers. Sequences were edited with Genalys v3.3.45a (M. Takahashi). Primers and PCR conditions are available on request. Sequences have been deposited in GenBank, accession numbers DQ821585-DQ821701 for DCT, EF675246-EF675361 for TYR and EF675362-EF675477 for TYRP1.

Estimation of haplotypes

To solve the haplotypes phase, we first run PHASE [49]. For those pairs of SNPs that did not reach a PHASE probability greater than 0.95, we solved their phase experimentally by ARMS-PCR and/or cloning (using the TOPO-TA kit from Invitrogene) plus resequencing.

Population parameters and neutrality tests

After removing polymorphic positions in humans, divergence (K) between a chimp sequence and a human sequence was estimated using K-estimator 6.0 [50]. Initially, population diversity parameters and neutrality statistics like Tajima's D [29] were obtained by means of DnaSP 4.1 [51]. These tests were corrected for demography as specified below. HEW and DHEW tests [32] were carried out using software kindly provided by Kai Zeng.

Optimization of demographic parameters

To correct for demography in the coalescent neutral simulations of the neutrality tests (excluding HEW and DHEW), we optimized the fit between pairwise FST distributions obtained from real genomic data from the three major geographical human groups in the HapMap project, and those FST distributions from coalescent simulations obtained varying the demographic parameters. The optimization criterion was the p-value of the Kolmogorov-Smirnov D statistic between the real and simulated FST distributions. For the real neutral distribution of the FST statistic [52], we used that obtained previously in [15] in which we selected 43 regions distributed across the autosomal genome that belong to broader regions of low gene density, and which are at least 150 kb away from the closest exon. Each of these regions spanned an average of 1.96 Mb, and in total they account for 84.3 Mb. For each region we downloaded the SNP frequency information available from the HapMap browser (data Rel #20/phase II on NCBI B35 assembly, dbSNP b125) for the 3 major populations (Caucasians: 153,339 SNPS, Yorubas: 123,798 SNPs and Chinese: 33,190 SNPs). We further filtered the number of SNPs to include only those SNPs that: a) were at least 100 kb away from each other, b) have been genotyped in all three populations and c) at least one of the three major populations had a minor allele frequency (MAF) higher than 0.1. A final list of 546 SNPs satisfied these criteria. We used the ms program [53] for the simulations with 3 populations, with sample sizes equal to those in the HapMap population (120 chromosomes for the African population, 120 for Caucasians and 90 for Asians). As starting points in our simulations to optimize demographic parameters, we used those values described in [54]. Simulations were fixed on one segregating site. These SNPs were matched for a MAF of 0.1 in at least one population.

The mean and 95% upper limits (between brackets) of the observed FST distributions were [15]: Caucasians-Chinese: 0.08 (0.33); Caucasians-Yorubas: 0.14 (0.47) and Chinese-Yorubas: 0.16 (0.45). Final optimized values were obtained in the simulations under the following conditions: we used an ancestral population size of 24,000 for the African population (population 1) and 7,700 for both Asians and Caucasians (populations 2 and 3, respectively), with a migration rate matrix Mij = {0, 0.05, 0.4, 0.1, 0 3, 0.8, 2.5, 0} for i and j values from 1 to 3. Looking back in time, we assumed two bottlenecks with instant population reduction, each followed by a population fusion: one approximately 40,000 years ago, in which the Chinese population reduced its size to approximately one sixth. About 2,000 years after this episode, the Chinese population fuses with the European population. Assuming a generation time of 20 years, this represents an F value of 0.04 for this bottleneck. A second population bottleneck takes place about 90,000 years ago. On this occasion, the Eurasian population suffers a reduction in size to one sixth of its previous size. About 10,000 years after this bottleneck (F = 0.21), the Eurasian population fuses with the African population. The mean and 95% upper limits (between brackets) of the simulated FST distributions were: Caucasians-Chinese: 0.08 (0.31); Caucasians-Yorubas: 0.15 (0.44) and Chinese-Yorubas: 0.15 (0.45). The optimized, simulated FST distribution and the real distribution recorded Kolmogorov-Smirnov D values of 0.0356 for Caucasians vs. Asians (p-value 0.891), 0.0748 for Caucasians vs. Africans (p-value 0.104), and 0.0441 for Asians vs. Africans (p-value 0.683). As both simulated and real data are not statistically different, we used those demographical parameters used in the simulations for the subsequent neutrality tests.

Correction for demography in the neutrality tests

These demographic and genetic parameters were subsequently used in further simulations for estimating the critical points in Tajima's D. These simulations also allowed us to obtain the distribution of Tajima's D under the genetic and demographic parameters described above.

Extended Haplotype Homozygosity (EHH) test

For the EHH test, we initially used Sweep 1.0 [55] to scan and select the core haplotypes. We initially focused on those haplotypes that showed both substantial frequency and high EHH values from the core SNPs, in both 5' and 3' directions. By means of a Perl script, we then calculated the EHH values as in [26] for a region extending about 50 kb from the closer core SNP. To test the significance of the test, we ran coalescent simulations using the population and demographical parameters described above for the FST distributions, but in this case several aspects are different from our previous approach [15]. First, in this case we included the variable recombination rate information obtained from the HapMap webpage for the particular region tested. In each case, we obtained a discrete number of recombination rate classes using the following approach: for each region we obtained the average and standard deviation (sd) of the distribution of recombination rates obtained from the HapMap link. All recombination values greater than the average plus 2 sd were considered outliers. If these outlier regions were consecutive, a local average was calculated; otherwise, a single outlier value was assigned for that region flanked between the previous and the next recombination values. We then excluded these outliers and repeated the process. In this second round, recombination rates that were higher than the new global average plus 1 sd were considered again as a class each, except if they were consecutive, in which case a local average was estimated. The average of the rest of recombination rates was considered as the background recombination rate, which was used as a reference to estimate the relative intensity of recombination for all other recombination classes.

For the coalescent simulations with heterogeneous recombination rates, we used msHOT [56], a modification of [53] coalescent-based program (ms) for simulating genetic variation data for a sample of chromosomes from a population. Second, in order to obtain a null distribution that reflects a neutral scenario for the chosen core haplotypes, we proceeded as follows: a) in the simulations, the "core" was the set of the first n SNPs, where n is the number of SNPs in the original core for each case. We then discarded those simulations that did not result in a number of distinct haplotypes (as defined from the simulated core set of SNPs) identical to that observed in the HapMap data using the corresponding observed core SNPs; b) one of the simulated haplotypes had to match in frequency (allowing for a 2% difference) our observed core haplotype being tested, and in addition, it had to show the same ancestral/derived states for the SNPs composing the core. For the latter, ancestrality was obtained by comparing the corresponding orthologous regions from the chimpanzee genome sequence and the Macaca mulatta genome sequence obtained from the UCSC genome browser [57] or the Ensembl genome browser [58]; c) for each set of simulations that fulfilled these conditions (we typically ran the program until we obtained about 500 simulations satisfying all the conditions), we obtained the 95% upper percentile of the expected EHH distribution for each of a series of consecutive 1 kb-long windows spanning the region. We finally declared a core haplotype as under selection if the distance for which the EHH values were equal to 1 for the SNPs in the HapMap data was longer than the distance observed from the 95% upper percentile distribution of the EHH values obtained in the simulations.

This approach offers several advantages over other methods currently being used to detect selection. For instance, it allows for a direct comparison with a "neutral" null distribution. In contrast, null distributions obtained from genomic regions are not representative of a homogeneous neutral scenario, but rather the result of heterogeneous evolutionary processes. In addition, it allows statistical inference even when no alternative haplotypes for the same core SNPs are available in order to obtain relative EHH values, as done in other approaches. Finally, ancestral-derived state information and information on the core haplotypes frequency can be incorporated, which helps fine-tune the nature of the elements being compared (observed and simulated data).

Multiple testing corrections

To control for multiple testing in the FST tests, we used the approach by [59], which sets the false discovery rate (FDR) at a level α by ranking the initial p values in ascending order P(1) ≤ P(2) ≤ ... ≤ P(m), with m being the number of tests, and then by specifying P(i) ≤ αi/m as the point below which there is no rejection at an FDR of α.

Genealogical relationships among haplotypes

Graphical representations of the genealogical relationships among haplotypes were estimated by the Median-Joining (MJ) algorithm implemented in Network 4.1.1.2 [60]. When feasible, the ancestral haplotype was inferred using parsimony by comparison to the chimp and rhesus macaque sequences.

Test for overdominance

To evaluate the possibility of overdominance (heterozygote advantage), we scored the ratio of "observed heterozygosity" to "expected heterozygosity" for single SNPs. Observed heterozygosity was estimated by counting heterozygote individuals in the HapMap data set for each SNP in question. The expected heterozygosity was estimated by calculating gene diversity from allele frequencies. Significance for ratio values greater than 1 was obtained by simulation using ms [53].