Gene Expression Profiles in Peripheral Lymphocytes by Arsenic Exposure and Skin Lesion Status in a Bangladeshi Population

Arsenic contaminates the drinking water supply at a concentration greater than 50 μg/L for millions of people worldwide, with ∼35 to 50 million of these individuals in Bangladesh (1, 2). Arsenic is a well-established human carcinogen (3); however, the exact mechanism by which it causes cancer has not been established. This is partly because there is no good animal model to study the toxicity of arsenic. Many of the human health effects of arsenic have been established based on epidemiologic studies, which have shown a significant association between the consumption of arsenic through drinking water and cancers of the skin, lung, bladder, liver, and kidney (4-8), neurologic disease (9), cardiovascular disease (10), and other nonmalignant diseases (11, 12). Premalignant skin lesions (melanosis, leucomelanosis, and keratosis) are an early manifestation and hallmark of arsenic toxicity and may indicate increased future risk of arsenic-related cancer (13). However, the molecular basis of arsenic-induced skin lesions and its progression to cancer is poorly understood.

Several microarray-based gene expression studies have been conducted to investigate the mechanism of arsenic toxicity (14-39). Only two of these studies were conducted in human populations (27, 35), whereas the rest were in vitro or animal studies. Lu et al. (27) examined aberrant gene expression associated with arsenic-induced liver disease using cDNA microarrays that contained ∼600 genes and found ∼10% differentially expressed. The study used RNA that was extracted from liver biopsies of six exposed and six unexposed individuals. Wu et al. (35) examined differential gene expression associated with intermediate and high arsenic exposure (determined by blood arsenic concentration, 4.64-46.5 μg/L) compared with low arsenic exposure (blood arsenic concentration, 0-4.32 μg/L) using cDNA arrays that contained 708 known gene transcripts and found ∼9% differentially expressed. RNA was extracted from peripheral blood lymphocytes of 24 participants. Both of these studies explored gene expression of a limited subset of genes, with <1,000 transcripts interrogated.

We conducted a microarray-based gene expression study among individuals chronically exposed to arsenic to assess whether arsenical skin lesion status and arsenic exposure level are associated with differential gene expression patterns. This study is the first, to our knowledge, to examine the effect of arsenic on genome-wide expression in humans. Through identifying differential patterns of expression, we hope to characterize the effect of arsenic and its outcomes, to better understand arsenic carcinogenesis, and to identify potential molecular targets for chemoprevention.

Forty individuals were selected from the Health Effects of Arsenic Longitudinal Study (HEALS) cohort for examination of gene expression profiles (40). The HEALS is an ongoing, population-based study examining both short-term and long-term health effects of arsenic exposure from drinking water in Bangladesh. In 2000, the study was launched in Araihazar, Bangladesh, among a population of individuals chronically exposed to arsenic through groundwater consumption. The 40 individuals selected for this pilot study included 15 individuals with arsenical skin lesions and 25 individuals without such lesions. At the baseline interview of HEALS cohort participants, a clinical examination was conducted by study physicians to diagnose and identify skin lesions. Skin lesions included melanosis, leucomelanosis, and keratosis. Melanosis was characterized by the hyperpigmentation of the skin over wide body surface areas. Leucomelanosis was characterized by both hyperpigmentation and hypopigmentation of the skin over wide body surface areas. Keratosis was characterized by the general thickening of the skin of the palms and soles (41). The 25 individuals in this pilot study with no manifest arsenical skin lesions were sampled from the cohort by matching on gender and age (±5 years) to the individuals with skin lesions.

A blood sample was collected from each participant at the time of recruitment into the parent HEALS cohort between October 2000 and May 2002; however, because these blood samples were not processed and stored in a manner suitable for isolating RNA, a separate blood sample was collected from each of these 40 individuals between February 2003 and March 2003, specifically for the examination of gene expression in this pilot study. Individual-level arsenic exposure was assessed based on urine samples collected as part of the parent HEALS follow-up between December 2002 and March 2004 and used to examine correlation with gene expression patterns in this study. Urinary total arsenic concentration was analyzed by graphite furnace atomic absorption spectrometry according to the method of Nixon et al. (42). Urinary creatinine was measured by a colorimetric Sigma Diagnostics kit (Sigma, St. Louis, MO), and urinary total arsenic concentration was subsequently expressed as micrograms per gram creatinine. Among the 40 participants, 28 subjects had both blood and urine samples used for this analysis collected within a 3-month period of each other (range, <1-11 months). From questionnaire data available at the time of urine sample collection, 28 participants (3 arsenical skin lesion cases and 25 nonlesion participants) indicated they drank from the same well as reported at enrollment into the parent cohort study and 6 arsenical skin lesion participants reported having switched their well before the time of both blood and urine sample collections used in these analyses. For these 34 individuals, the urine sample should reasonably reflect arsenic exposure at the time of blood collection. The remaining 6 arsenical skin lesion cases switched their well in the period of time between the blood and urine sample collections, for whom correlation between blood and urine samples may be problematic.

RNA was isolated from peripheral blood lymphocytes, processed, and analyzed for gene expression of each individual; the 40 samples were individually run requiring 40 microarrays. In these analyses, gene expression comparisons were only evaluated across arrays with a similar 3′/5′ ratio for both glyceraldehyde-3-phosphate dehydrogenase and β-actin genes on the microarray. A 3′/5′ ratio measure greater than 3 for these internal control genes indicates degraded RNA or incomplete in vitro transcription. Among the 40 microarrays, a single batch of 20 samples processed using the same in vitro transcription kit all had 3′/5′ ratios greater than 3 for these internal control genes. All samples in the batch were from individuals with no manifest arsenical skin lesions and therefore were only used in analyses examining exposure level and not used in the comparison by skin lesion status. Among the 20 arrays run in the other in vitro transcription batch, all arrays had 3′/5′ ratio measures less than 3 for the internal control genes; however, 4 arrays (samples from individuals with skin lesions) were excluded because they contained a background signal greater than 100. Arrays with average background signals greater than 100 were taken to indicate a high level of nonspecific hybridization on those chips and excluded to preserve data quality. The remaining 16 arrays (11 skin lesion and 5 nonlesion samples) were comparable on other microarray quality metrics (43, 44), including noise (RawQ) <5, background signal <100, consistent scale factors, consistent detection of BioB and BioC spike controls, and consistent number of genes detected as present.

In subgroup exploratory analyses, subjects with no manifest arsenical skin lesions were subsequently categorized based on their creatinine-adjusted urinary total arsenic concentration. This study objective was explored separately in the two sets of nonlesion arrays (n = 5 and 20) to ensure the arrays were comparable in microarray quality variables as described previously. In one subgroup analysis (n = 5), 2 arrays were assigned low arsenic exposure based on a creatinine-adjusted urinary total arsenic concentration <70 μg/g and 3 arrays assigned high arsenic exposure based on a creatinine-adjusted urinary total arsenic concentration >70 μg/g. This cut point was selected based on the distribution of creatinine-adjusted urinary total arsenic concentration in this nonlesion study sample (n = 5). The lowest urinary arsenic tertile (corresponding to the concentration range 48-67 μg/g) was categorized as low arsenic exposure. The highest two urinary arsenic tertiles (corresponding to the concentration range 90-215 μg/g) were categorized as high arsenic exposure. In the second subgroup analysis (n = 20) exploring this objective, 6 arrays were assigned low arsenic exposure based on a creatinine-adjusted urinary total arsenic concentration <162 μg/g and 14 arrays assigned high arsenic exposure based on a creatinine-adjusted urinary total arsenic concentration >162 μg/g. This cut point was defined based on the distribution of creatinine-adjusted urinary total arsenic concentration in these nonlesion samples (n = 20). The lowest urinary arsenic tertile (corresponding to the concentration range 40-161 μg/g) was categorized as low arsenic exposure. The highest two urinary arsenic tertiles (corresponding to the concentration range 166-880 μg/g) were categorized as high arsenic exposure.

Whole venous blood samples were collected in 10 mL Vacutainer tubes containing 0.117 mL 15% (K3)EDTA and temporarily stored in a cooler box for transportation to the Dhaka laboratory. All blood samples were processed the day they were collected; processing took place within 2 to 8 hours after collection. Blood (8 mL) was transferred to a DNase- and RNase-free 15 mL tube and centrifuged at 2,400 rpm for 10 minutes at 25°C. Supernatant plasma was removed and preserved. Sterile 1× PBS stored at room temperature was added up to 15 mL and the blood was gently resuspended. The blood was then very slowly layered on approximately half-volume Ficoll-Paque (Amersham, Piscataway, NJ) and centrifuged at 2,700 rpm for 25 minutes. The opaque interphase layer, containing mononuclear cells, was transferred to another 15 mL DNase- and RNase-free tube and 1× PBS was added up to 15 mL. The tube was centrifuged at 1,000 rpm for 12 minutes and PBS was removed. Either 2 mL TRIzol reagent (Invitrogen, Carlsbad, CA) or 3.6 mL RLT buffer (Qiagen, Valencia, CA) was added to the mononuclear cells. After collection using syringe for homogenization, the mononuclear cells in TRIzol/RLT buffer were transferred to nonstick DNase- and RNase-free tubes for storage at −80°C. Within 2 to 3 months, the samples were transported on dry ice by air to the United States where they were stored at −80°C until RNA extraction.

RNA was extracted using the QIAamp RNA Blood Mini kit (Qiagen). The initial few steps were different for TRIzol and RLT preserved samples; however, the remainder was similar. TRIzol samples were incubated at 30°C for 5 minutes and chloroform (200 μL) was added per 1 mL TRIzol reagent. Samples were shaken vigorously by hand for 15 seconds and then incubated at 30°C for 3 minutes. The samples were then centrifuged at 10,000 rpm for 15 minutes at 4°C. The aqueous phase containing RNA was carefully removed to a new RNase-free tube and 700 μL RLT was added per 200 μL aqueous phase and mixed well. Then, 500 μL molecular-grade absolute alcohol was added and mixed by pipetting. The sample was put on a QIAamp spin column sitting in a 2 mL collection tube. The rest of the procedure was conducted as per the Qiagen recommendation in the QIAamp RNA Blood Mini kit for RLT samples (described below). The RLT preserved samples were incubated at 37°C for 10 minutes. The remainder of the procedure was carried out at room temperature. For every 600 μL RLT sample, one QIAshredder column was used for centrifugation at 13,000 rpm for 2 minutes. The filter was discarded and 70% molecular grade alcohol (600 μL) was added per 600 μL filtered sample in RLT buffer and mixed well. The sample was then put on QIAamp spin column sitting in a 2 mL collection tube. The subsequent steps were similar for both TRIzol and RLT preserved samples. In short, the sample-filled QIAamp spin column was centrifuged and the filtrate was discarded. RNA was bound to the filter membrane. The membrane was washed with RW1 buffer, treated with DNase, and washed again with RW1 buffer followed by RPE buffer. RNA was eluted in RNase-free water.

Using a 1 μL sample, RNA quality and quantity were measured using the RNA 6000 Nano Chips in the Agilent 2100 Bioanalyzer (Palo Alto, CA). The A₂₆₀/A₂₈₀ ratio was checked using 1 μL RNA on the ND-1000 spectrophotometer. All samples had A₂₆₀/A₂₈₀ ratio > 1.9, indicating high-quality RNA. RNeasy MinElute Cleanup kit (Qiagen) was used to further concentrate the RNA. Integrity of RNA was checked again on the Agilent 2100 Bioanalyzer (an example is shown in Fig. 1). As seen in Fig. 1A, the two most distinct and intense peaks on the electropherogram correspond to the 18S and 28S rRNA. A 28S/18S ratio > 1.4 and the sum of both rRNA peak areas accounting for >33% of all RNA were taken as indicators of good-quality RNA.

The Affymetrix (Santa Clara, CA) standard protocol has been described extensively elsewhere (45). Briefly, starting with 5 to 8 μg total RNA, we generated double-stranded cDNA using T7-(dT)₂₄ primer 5′-GGCCAGTGAATTGTAATACGACTCACTATAGGGAGGCGG-(dT)₂₄-3′ and SuperScript Choice System for cDNA synthesis (Invitrogen). The cDNA samples were purified using GeneChip Sample Cleanup Module (Affymetrix). Then, the in vitro transcription was done using the Enzo BioArray High-Yield RNA Transcription Labeling kit (Affymetrix) to produce biotin-labeled cRNA. Purification of cRNA was done using the GeneChip Sample Cleanup Module (Affymetrix). Quality and quantity of purified cRNA samples were checked by running a 1 μL sample on the ND-1000 spectrophotometer as well as on the Bioanalyzer. A A₂₆₀/A₂₈₀ ratio between 1.9 and 2.2 was accepted. The adjusted cRNA yield was calculated as below:

where RNAm is the amount of purified cRNA measured (μg), total RNAi is the starting amount of RNA (μg), and y is the fraction of cDNA used in the in vitro transcription reaction.

We fragmented 15 μg adjusted cRNA product using the 5× fragmentation buffer supplied in the GeneChip Sample Cleanup Module. Fragmented cRNA was run on the Bioanalyzer to ensure the fragmented size. Typical electropherograms of unfragmented and fragmented cRNA are shown in Fig. 1B and C, respectively. Biotin-labeled fragmented cRNA was then supplied to the core facility for hybridization onto the GeneChip HG-U133A array. After 16 to 18 hours of incubation, each microarray was washed, stained with streptavidin-phycoerythrin, and scanned by the high-resolution Affymetrix GeneChip Scanner 3000 according to GeneChip Expression Analysis Technical Manual procedure (Affymetrix). After scanning, the raw intensities of each probe were stored in electronic files (.DAT and .CEL formats) by GCOS software (Affymetrix).

The robust multiarray average method was applied to extract probe set signal intensities and to normalize all arrays (46). Before statistical analyses, genes were filtered by applying two exclusion criteria: (a) genes showing minimal variation (below the 80th percentile) across the set of arrays and (b) genes with >50% missing data were excluded from statistical analysis. Gene filtering was conducted separately for each distinct aim of the study and yielded separate lists of 4,456 genes each, which were used for the statistical analyses. Differentially expressed genes were identified using the significance analysis of microarrays statistic with a false discovery rate of 1% (47). The false discovery rate is the proportion of genes claimed to be differentially expressed that are false positives. Gene Ontology (GO) and pathway comparisons were also conducted based on the lists of 4,456 filtered genes using functional class scoring analysis as described by Pavlidis et al. (44). For each gene in a GO category, a P for differential expression comparing the two conditions of interest was computed. The set of Ps for a GO category were summarized by two statistics: (a) the Fisher statistic P and (b) the Kolmogorov-Smirnov statistic P. A GO category was considered significantly differentially expressed if either significance level was <0.005. All GO categories with between 5 and 100 genes represented on the array were considered. Functional class scoring is a more powerful method of identifying differentially expressed gene classes than overrepresentation analysis or annotation of gene lists based on individually analyzed genes (44). Analysis of pathway comparisons were similar to the GO analysis; however, genes were grouped by KEGG or BioCarta pathways instead of GO categories (48). All analyses were done using BRB ArrayTools version 3.3.0 Beta 2 (49).

There is currently no predominate method to handle confounding variables in gene expression microarray studies. An effective technique for eliminating confounding in traditional epidemiologic studies is to match on potential confounding variables in the study design (50). For the current study, gender and age were considered potential confounders of gene expression in relation to arsenic exposure and skin lesion status. Their effects were considered a priori in the design stage of the study and matched on in the study participant selection. Additionally, analyses were restricted to never-smokers; smoking status was not included as a matching factor because there was no a priori reason to believe that this variable was a potential confounder because it is not an independent risk factor for arsenical skin lesions nor is it associated with arsenic exposure status.

Gene expression profiles, using RNA derived from peripheral blood lymphocytes, were analyzed using Affymetrix HG-U133A GeneChip arrays. This microarray platform contains 22,283 probe sets for interrogation, including known genes and expressed sequence tags (51). To examine differential expression associated with arsenic-induced skin lesion status, the blood lymphocyte gene expression profiles of the subjects with no manifest skin lesions (n = 5) were compared with subjects with arsenical skin lesions (n = 11). The characteristics of the study participants are shown in Table 1. The mean (SD) well water arsenic concentration was 342.7 (258.1) μg/L for individuals with skin lesions and 39.6 (49.5) μg/L for individuals without skin lesions. The number of differentially expressed genes generated from the significance analysis of microarrays statistic was 468 (467 down-regulated and 1 up-regulated genes in the skin lesion group). Restricting the study sample to females only resulted in 330 differentially expressed genes (all down-regulated in the skin lesion group) between the skin lesion (n = 9) and nonlesion (n = 5) groups. Further restricting the study sample to female never-smokers resulted in 312 differentially expressed genes (all down-regulated in the skin lesion group) between the skin lesion (n = 8) and nonlesion (n = 5) groups. The findings of the female never-smoker subgroup analysis are reported here, because we consider this gene list to have the least influence of gene expression variability due to extraneous factors. However, there was considerable overlap across the three gene lists generated; all 312 differentially expressed genes were present in the gene list from the overall analysis and 310 were present in the gene list from the gender-restricted analysis. The top 50 down-regulated genes are listed in Table 2.

GO analysis, based on functional class scoring, revealed 41 significant GO categories, which are listed in Table 3. From the pathway comparison, one significant pathway was identified that could discriminate between the skin lesion (n = 11) and the nonlesion (n = 5) classes (Table 4). Additional analyses within the skin lesion group only (n = 11) to assess the influence of arsenic exposure based on well water and urinary arsenic concentrations did not reveal any differentially expressed genes. Analyses by skin lesion severity also did not reveal any differentially expressed genes.

To explore differentially expressed genes associated with arsenic exposure, the gene expression profiles from blood lymphocytes of individuals with low arsenic exposure were compared with those from individuals with high arsenic exposure, indicated by creatinine-adjusted urinary total arsenic concentration. Arsenic exposure was defined by tertile of creatinine-adjusted urinary total arsenic concentration for the arrays. These analyses were conducted separately for the two nonlesion sample populations, with arrays restricted to in vitro transcription batch (n = 5 and 20) to ensure comparability in microarray quality variables. The low and high arsenic exposure groups were similar in age, gender, and smoking status (data not shown). Within the 5-array nonlesion sample, the mean (SD) creatinine-adjusted urinary total arsenic concentration was 144.7 (63.9) μg/g for the high arsenic exposure group and 57.5 (13.4) μg/g for the low arsenic exposure group. Within the 20-array nonlesion sample, the mean (SD) creatinine-adjusted urinary total arsenic concentration was 346.6 (250.5) μg/g for the high arsenic exposure group and 123.5 (44.2) μg/g for the low arsenic exposure group. Statistical analysis did not reveal any differentially expressed genes by urinary arsenic exposure status in either sample population. Recategorization and analysis by well water arsenic concentration status in both sample populations yielded no differentially expressed genes.

To our knowledge, this is the first molecular epidemiologic study to explore alterations in genome-wide expression due to arsenic exposure and arsenic-induced toxicity in a currently exposed human population. Only two prior in vivo studies had been conducted, which explored a subset of genes differentially expressed in relation to arsenic exposure or arsenic-induced disease (27, 35). Using the Affymetrix HG-U133A GeneChip, we have identified differentially expressed genes in peripheral blood lymphocytes of arsenic-exposed individuals and individuals with arsenical skin lesions.

In the present study, the effect of arsenic exposure was assessed among a population of individuals who have been chronically exposed to arsenic through groundwater consumption and were currently exposed at the time of blood sample collection for this pilot study. The HEALS cohort provided the sampling frame to examine individuals exposed to arsenic with no clinical evidence of premalignant skin lesions compared with individuals with manifest arsenical skin lesions. Thus, the effect of both arsenic exposure and exposure-related disease process could be examined. Study participants were matched on age and gender to minimize variability in gene expression that may be due to these confounding factors. Despite the restriction of samples in these analyses to in vitro transcription batch, the comparison groups still remained comparable on gender and age. Furthermore, the samples were collected from a population of individuals who are relatively homogenous in ethnicity and dietary intake, which again served to minimize variability due to extraneous factors. Although smoking status was not included as a matching factor, the effect of smoking was removed by restricting analyses to never-smokers.

In this study, we evaluated the joint effect of arsenic exposure and arsenical skin lesion status on gene expression. GO categories found to be differentially expressed between individuals with arsenical skin lesions compared with exposed individuals without such lesions broadly included RNA metabolism, hydrolase activity, ribonucleoprotein complex, translation, cellular protein catabolism, amino acid activation, transport and transporter activity, and glycoprotein metabolism. Additionally, signal transduction through the interleukin (IL)-1 receptor was identified as a significant pathway of differentially expressed genes between the arsenical skin lesion and nonlesion groups. Research has shown that this pathway plays an important role in skin (52, 53), including the skin repair process (54). The IL-1 receptor pathway may lead to the activation of nuclear factor-κB and p38 pathways (54); thus, the down-regulation of which may be indicative of deficiency in wound healing among the arsenical skin lesion cases.

Most genes differentially expressed in analyses between arsenical skin lesion cases and noncases were down-regulated. However, in the analysis between arsenical skin lesion cases (n = 11) and noncases (n = 5), we found 3.3-fold up-regulation of thymopoietin, a gene that has been implicated in the initiation of DNA replication (55). Furthermore, an in vitro study found that thymopoietin was overexpressed in cancer cell lines (56).

We have identified several genes that have been implicated in carcinogenesis. The superoxide dismutase 2 gene, located at chromosomal region 6q25.3, encodes a free radical–scavenging enzyme that is the first line of defense against oxidative stress. Arsenic has been implicated as a source of reactive oxygen species (57-59). Additionally, polymorphisms in this gene have been associated with cancer (60). Some evidence from in vitro studies suggests that superoxide dismutase 2 may act as a tumor suppression gene (61-65), with increased expression of superoxide dismutase 2 having been shown to suppress the cancer phenotype in a large number of mammalian tumor cells, including human melanoma cells (61, 62). We observed down-regulated expression of superoxide dismutase 2 in individuals with skin lesions, which may indicate an increased vulnerability to reactive oxygen species generated by arsenic in individuals with manifest arsenic toxicity. A second gene of interest is tumor necrosis factor (TNF) that is located at chromosomal region 6p21.3. This gene encodes a cytokine that is involved in cell proliferation, differentiation, apoptosis, lipid metabolism, and coagulation. Furthermore, studies have shown arsenic to be an inducer of TNF (66, 67). TNF is thought to trigger a chemokine signal response, including chemokine factor CCL20 (68). In this study, we have observed down-regulated TNF expression among individuals with arsenical skin lesions and down-regulated expression of CCL20, providing credence to the suppression of a chemokine response pathway in arsenic toxicity. There is evidence to indicate that TNF is released early in response to skin wounds or inflammation (69), which may suggest that the down-regulated expression of TNF observed in arsenical skin lesion cases is also indicative of deficient wound healing.

In future analyses, we plan to examine the correspondence of the current study findings using peripheral blood lymphocytes to a study examining gene expression from the target tissue of interest (e.g., skin lesion tissues). It is currently not well known whether gene expression in blood lymphocytes is a valid proxy for gene expression in the target tissue of interest and needs to be explored in future research.

A limitation of the present study was the extent of incomplete in vitro transcription among the nonlesion samples, which limited our ability to explore differential gene expression by arsenic exposure status only. After RNA extraction, sample quantity and quality was assessed, which indicated good RNA quality for all samples; therefore, the poor 3′/5′ ratio in the glyceraldehyde-3-phosphate dehydrogenase and β-actin genes are indicative of incomplete in vitro transcription occurring in an entire batch as opposed to sample degradation in this study. The current analysis was limited to separate comparisons of 5 and 20 samples from individuals with relatively varied arsenic exposure. The lack of differentially expressed genes from these comparisons was most likely due to small sample size and heterogeneous exposure categories as opposed to the absence of a biological effect by arsenic on gene expression. In future studies, we hope to evaluate the effect of arsenic exposure in a dose-dependent manner on gene expression.

Although the Affymetrix HG-U133A GeneChip array measured genome-wide expression, the use of this platform assessed mRNA from coding regions of the genome only. Recent evidence has suggested the importance of RNA from noncoding regions that may play key roles in carcinogenesis and other processes (70-72). A limitation of the current microarray technology used in this experiment (or any other experiments using generic reverse transcription-PCR) is that only RNA with a polyadenylic acid tail is amplified and reverse transcribed. A recent study reported that a substantial proportion of transcribed messages in the genome may include RNA without polyadenylic acid tail (73) and thus may not be detected in the currently available microarray-based gene expression analyses.

In conclusion, our findings show that microarray-based gene expression analysis is a powerful method to characterize the molecular profile of arsenic exposure and arsenic-induced diseases. Genes identified from this analysis may provide insights into the underlying disease mechanisms of arsenic and represent targets for disease prevention.