Introduction Systemic lupus erythematosus (SLE) is a chronic autoimmune disease with a complex etiology and diverse clinical manifestations [1]. The role of genetic factors in the SLE risk has long been established, and demonstrated in familial aggregations, twin studies, and sibling recurrence rates [2]. Recently, high-throughput technologies have facilitated genome-wide association studies (GWASs) across different populations. This approach, accompanied by large-scale replications, has not only confirmed the association of many established susceptibility genes, but has also presented convincing evidence of novel genetic loci involved in SLE [3]–[8]. As members of the Asian Lupus Genetics Consortium, we have also performed a GWAS in Asian populations and have identified variants in ETS1 and WDFY4 that are associated with SLE [9]. A combination of GWAS data from different ethnic groups will clearly provide new insights into the genetics of SLE and further our understanding of the pathogenesis of lupus [10], [11]. To use genomic tools to study the mechanisms of SLE, we and others have independently identified a gene expression signature for lupus patients using microarray profiling [12]–[14], which highlights the pathogenic role of the abnormal activation of the type I interferon (IFN) pathway in human lupus [15]–[17]. Intriguingly, recent investigations suggest a genetic contribution to the variability observed among individuals in the production and signaling of IFN [17], and advances in the genetics of SLE highlight the strong association between the risk of developing lupus and gene variants connected to the production and effects of type I IFN [11], [18]. We recently used a microRNA (miRNA) profiling assay to examine the involvement of miRNAs in SLE, because miRNAs are novel gene expression regulators [19] and important players in shaping the immune response [20]–[22]. This profiling identified a reduction in miR-146a expression in lupus patients, and we showed that the underexpression of miR-146a contributes to lupus pathogenesis by deregulating the activation of the IFN pathway [23]. However, why miR-146a levels are reduced in patients with SLE remains unresolved. miR-146a is encoded at 5q33.3. Interestingly, recent data from GWASs in both European and Asian populations have indicated that this region is a novel susceptibility locus for SLE [3], [7], [9], suggesting a plausible role for a genetic variant around miR-146a in modulating its expression and thus the disease risk. Several studies have demonstrated unambiguously that genetic variants in miRNA precursors (pre-miRNA) can affect miRNA expression levels by interfering with the miRNA maturation process and are thus associated with disease susceptibility [24]–[26]. We postulate that genetic variants in both the miRNA promoter and the precursor region may alter mature miRNA production. Given the critical regulatory role of miR-146a in the type I IFN pathway and the abovementioned genetic association between this pathway and SLE susceptibility, polymorphisms in the miR-146a gene could also potentially confer a disease risk. To assess whether genetic variants modulate miR-146a expression and thus contribute to the risk of developing SLE, we sequenced the promoter and key regulatory regions of the miR-146a precursor to identify potential functional variants that might be associated with SLE susceptibility. Our subsequent replication and functional studies provide evidence that single-nucleotide polymorphism (SNP) rs57095329 in the miR-146a promoter, which affects its mature level, can confer SLE susceptibility. Results Discovery of SLE-associated miR-146a promoter SNPs miR-146a is located at 5q33.3. The transcription start site (TSS) of its primary transcript (pri-miR-146a) has been identified [27]. To characterize the essential regulatory region for subsequent genetic analysis, we first cloned miR-146a upstream fragments with variable 5′ ends into the pGL3-basic reporter plasmid to analyze its promoter activity. We found that the inclusion of a fragment from nucleotide (nt) −1,091 to nt −611, which contains a known NF-κB-binding site characterized in THP-1 cells [27], was consistently robust to promote luciferase activity in HeLa cells (Figure S1). The inclusion of the more distal region (nt −1,998 to nt −1,091) enhanced neither the basal nor phorbol myristate acetate and ionomycin (hereafter referred to as “PMA+Iono”) -induced activity of the promoter (Figure S1B ). Therefore, to look for new genetic variants and to characterize their potential association with SLE, we designed four pairs of primers with which to sequence the upstream region that spans the 1,105-bp promoter (nt −1,091 to nt +14) and the consecutive first exon of pri-miR-146a (Figure S1A ), in 360 individuals (180 SLE patients and 180 controls), which served as the discovery panel. We also sequenced the 452-bp region centered on miR-146a precursor or exon 2 of pri-miR-146a, because it potentially affects mature miR-146a production [26], in the same discovery panel. A total of 12 variants were identified, with nine already reported in the dbSNP database Build 130 (Table S1). Five variants had a minor allele frequency (MAF) of >1% (rs17057381, rs73318382, rs57095329, rs6864584, and rs2910164; Figure S1A ). Therefore, we extended our sequencing analysis to examine these five SNPs in up to 816 patients and 1,080 controls, who were all Chinese Han individuals living in Shanghai. In this expanded study panel, only rs73318382 and rs57095329 showed an association with SLE (Table S2). These two SNPs are separated by 304 bp and are in strong linkage disequilibrium (LD; r2 = 0.81; Figure S2). When a Bonferroni correction was applied, the association of rs57095329 with SLE remained highly significant (P = 4×10−4). Given that rs57095329 is identified through our candidate region sequencing approach and not included in the HapMap database, it is not surprising that this SNP has not been included in commercial SNP arrays. Because published GWASs in SLE of both Asian and European populations detected association signals at rs2431697 and rs2431099 (Figure S3), 15 kb and 8 kb upstream from miR-146a TSS, respectively, we extended our genotyping of rs2431697 and rs2431099 using 1,896 Shanghai samples. Both SNPs showed significant association with the disease (Table S2), while rs57095329 produced the best association signal among the three SNPs in the same dataset (Table S2). Therefore, we focused on rs57095329 in the subsequent experiments. Replicated association of rs57095329 with SLE in independent cohorts We replicated the association between rs57095329 and SLE using a TaqMan genotyping assay in another two panels from Hong Kong, China, and Bangkok, Thailand. We also added 1,536 patients from the central China area to our mainland China cohort, and the newly added patients showed an allele frequency for rs57095329 very similar to that in the discovery panel (MAFs of 20.53% and 20.77%, respectively). This replication provided consistent evidence for the association, revealed by an allelic association analysis (Table 1). When all the samples were included (3,968 patients and 3,214 controls in total) to conduct a meta-analysis, there was strong evidence that the minor G allele of rs57095329 conferred a risk of SLE (P meta = 2.74×10−8, odds ratio [OR] = 1.29, 95% confidence interval [CI] = 1.18–1.40; Table 1). There was no significant difference among the ORs for the three independent cohorts (P = 0.33), when the Breslow–Day test installed in PLINK was used [28], although the SNP showed significant allele frequency differences in respective controls. Recessive mode of action seemed to be supported in the Chinese mainland cohort and the cohort from Thailand (OR = 2.47 and 2.11 for the two cohorts, respectively), compared with the allelic OR of 1.35–1.36 for the two cohorts. However, this was not supported by the result for the Hong Kong cohort, where the same OR was observed for both the recessive mode and the allelic test (OR = 1.18), reflecting certain variations among the different cohorts. 10.1371/journal.pgen.1002128.t001 Table 1 Association between the rs57095329 G allele and SLE in independent cohorts and in the combined sample. Cohorts case/control Frequency χ2 P value OR(95% CI) case Control Mainland China 2,352/1,080 0.21 0.16 17.55 2.79E-5 1.35(1.17–1.55) Hong Kong 1,152/1,152 0.23 0.20 4.88 0.027 1.18(1.02–1.36) Bangkok 464/982 0.29 0.23 10.99 9.16E-4 1.36 (1.13–1.63) Meta-analysis 2.74E-8 1.29(1.18–1.40) We also examined whether the genetic variant is specifically associated with disease risk in patients with lupus nephritis. Although only the discovery panel in the Chinese mainland cohort showed a significant association in a patient-only analysis, a similar trend was also observed in the Hong Kong and Bangkok cohorts, with a marginal P value of 0.093 and an OR of 1.105 when patients with nephritis were compared with patients without it (Table S3). Association between rs57095329 and miR-146a expression We explored the association between rs57095329 and miR-146a expression. Mature miR-146a levels were determined with a TaqMan microRNA assay in 86 healthy controls with known genotypes and available RNA samples. Compared with individuals with the AA genotype, individuals with heterozygous AG genotype for rs57095329 had lower levels of miR-146a (P = 0.0438; Figure 1), while individuals with GG genotype had the lowest miR-146a levels (P = 0.0197; Figure 1). This association indicates that rs57095329, located in the miR-146a promoter, may function by regulating the transcription activity and expression levels of miR-146a. 10.1371/journal.pgen.1002128.g001 Figure 1 Comparison of miR-146a expression levels between groups of healthy individuals with different rs57095329 genotypes. The horizontal line indicates the mean expression level within each group. * indicates P 0.01) in the controls of all three cohorts. The average call rate for all samples was 92%. Real-time PCR Total RNA was extracted from peripheral blood leukocytes or cultured cells using TRIzol (Invitrogen), followed by reverse transcription using a reverse transcriptase kit obtained from Takara. To determine the quantity of pri-miR-146a, the cDNA was amplified by real-time PCR with SYBR Green RT–PCR kit (Takara), and the expression of RPL13A was used as the internal control. The primers used are shown in Table S7. To determine the quantity of mature miR-146a, the specific TaqMan MicroRNA Assay kit (Applied Biosystems) was used, and the expression levels were normalized to snRNA U6. The assays were performed on a 7900HT real-time instrument (Applied Biosystems). Relative expression levels were calculated using the 2−ΔΔCt method. Constructs To create the miR-146a promoter–luciferase reporter constructs, three fragments of variable lengths, corresponding to the upstream region of the TSS of pri-miR-146a, were amplified and cloned into the pGL3-basic luciferase vector (Promega). To compare the activities of miR-146a promoters containing the different rs57095329 alleles, the full-length 1105-bp fragment was amplified from individual homozygous templates. The ETS1 overexpression vector was a kind gift from Dr Gang Pei, and the PBX1 overexpression plasmid was created by replacing the inserted ETS1 sequence with the PBX1 coding sequence. The primers used are shown in Table S7. All constructs were verified by sequencing. Cell culture, transfection, and stimulation Jurkat and Raji cells were grown in RPMI 1640 medium supplemented with 10% fetal bovine serum. These two cell lines were electroporated with 2 µg of the indicated luciferase reporter vector and 0.2 µg of a modified pRL-TK vector, using a nucleofector device (Amaxa). Alternatively, the reporter gene vectors were electroporated in combination with 1.5 µg of an ETS1- or PBX1-expressing vector. For the knockdown of ETS1, 3 µg of ETS1 siRNA or negative control oligonucleotides (all from GenePharma, Shanghai) were transfected. HeLa and 293T cells were grown in Dulbecco's modified Eagle's medium supplemented with 10% fetal bovine serum. These two cell lines were transfected using Lipofectamine 2000 (Invitrogen), with the ETS1- or PBX1-expressing vector or ETS1 siRNA alone, or in combination with 50 ng of the indicated luciferase reporter vector and 5 ng of a modified pRL-TK plasmid. Where indicated, an irrelevant “carrier” vector was added to ensure that equal total amounts of plasmid DNA were transfected among the groups. For cell activation, Jurkat and HeLa cells were stimulated with PMA (100 ng/mL; Sigma) and ionomycin (1 µM; Sigma) for the indicated times. Alternatively, Jurkat cells were activated with plate-bound anti-CD3 antibody (coating solution: 5 µg/mL; eBioscience) and soluble anti-CD28 antibody (2 µg/mL; eBioscience). Reporter gene assay Cells were cultured for 24 hours or 48 hours after transfection with the reporter gene vectors together with the ETS1 expression vector or siRNA, respectively. The cells were then maintained resting or activated for 6 hours and lysed. Their luciferase activity was measured on a luminometer (LB960; Berthold) using the Dual-Luciferase Reporter Assay System (Promega). The ratio of firefly luciferase to Renilla luciferase was calculated for each well. EMSA Jurkat and HeLa cells (1×107) were activated or left to rest for 2 hours, and then their nuclear proteins were extracted with a Nuclear Extract Kit (Active Motif), according to the manufacturer's protocol. The protein concentrations were determined with the DC Protein Assay Kit (Bio-Rad). Double-stranded allelic probes were synthesized and labeled with biotin by Takara (the sequence is shown in Figure S6). EMSA was carried out with a gel-shift kit purchased from Active Motif. The competition assay was performed by adding cognate unlabeled oligonucleotides. After incubation, the protein–DNA complexes were separated on a nondenaturing 6% polyacrylamide gel and then transferred to a nitrocellulose membrane (Millipore). The signals were detected using a luminoimage analyzer. Streptavidin–agarose pulldown and western blotting The pulldown assay was performed following a protocol described elsewhere [40], with slight modification. Biotin-labeled allelic probes were incubated with equal amounts of nuclear extract from activated Jurkat cells for 2 hours at room temperature, in the presence of streptavidin–agarose beads (GE Healthcare) and protein inhibitors. The precipitated protein–DNA complex was dissociated from the agarose beads by suspending the pellet in Laemmli sample buffer (Bio-Rad) and heating it. The supernatants were then subjected to SDS–PAGE. The proteins were transferred onto a PVDF membrane (Bio-Rad), blotted with an anti-Ets-1 antibody, and detected with ECL solution (Pierce). To evaluate the Ets-1 protein levels after the transfection of the overexpression vectors or siRNAs, the Jurkat and HeLa cells were lysed in RIPA buffer (Thermo Scientific), and the supernatants were similarly used for immunoblotting. Anti-ETS1, anti-GAPDH, and horseradish-peroxidase-conjugated secondary antibodies were all obtained from Santa Cruz Biotechnology. Data analysis For single SNP analysis, PLINK was used for the basic allelic test and other tests in the patients and the controls [28]. LD patterns were analyzed and displayed with HaploView [41]. Review manager was used to perform meta-analysis. IMPUTE version 2 was used to perform imputation. Other data were analyzed with GraphPad Prism 4 software, version 4.03. The nonparametric Mann–Whitney test was used to compare miR-146a expression between the genotype groups, and an unpaired t test was used to compare reporter gene activities. Two-tailed P values 0.01 are shown. TSS, transcription start site. (B) Schematic representation of reporter gene constructs driven by various miR-146a upstream fragments (left) and their corresponding relative luciferase activities in HeLa cells (right), under rested (medium) and PMA+Iono-activated conditions. The data shown are means ± SEM and are representative of three independent experiments performed in triplicate. (TIF) Click here for additional data file. Figure S2 Linkage disequilibrium of six common SNPs in or upstream of the miR-146a promoter. Data are based on 816 SLE patients and 1,080 controls from Shanghai and were analyzed with HaploView. (PNG) Click here for additional data file. Figure S3 Plot of −log10 P values for SNPs genotyped in the GWAS spanning 5q33.3 region. Data were from three GWAS on both Asian (A) and European population (B and C). The linkage disequilibrium in the region derived from the Asian GWAS data is shown with r2 values as indicated. (PNG) Click here for additional data file. Figure S4 Reporter gene activity of constructs containing either rs57095329 allele in different cell lines. Shown are the relative luciferase activities of the two constructs driven by the miR-146a promoter containing either rs57095329 allele (A or G) in rested (medium) and activated (with the addition of PMA+Iono for 6 hours) HeLa cells (A), steady-state Raji cells (B), and steady-state 293T cells (C). The data shown are means ± SEM and are representative of three independent experiments performed in triplicate or quadruplicate. * indicates P<0.05, ** P<0.01, *** P<0.0001. (TIF) Click here for additional data file. Figure S5 Gel-shift assay of allelic probes with nuclear proteins from different cell lines. Shown is a comparison of the binding affinities of different rs57095329 alleles for the nuclear extracts (N.E.) from rested (medium) or PMA+Iono-activated HeLa cells (left), and from PMA+Iono-activated Jurkat cells (right). Also shown are the results of a competition assay, which was performed with the addition of 50- to 200-fold unlabeled cognate oligonucleotides. The assays were repeated at least three times. (TIF) Click here for additional data file. Figure S6 Predicted binding sites of Ets-1 on the miR-146a promoter. (A) The 200-nt sequence around rs57095329 (A/G) was used as the input for the Genomatix online tool, which predicted a nearby Ets-1-binding site (indicated by the blue box, with the red letters indicating the core sequence). Also shown is the probe sequence for the EMSA, indicated by the pink line below the sequence. (B) Conservation of rs57095329 residue and Ets-1 binding site. Shown is the UCSC Genome Bioinformatics search result by alignment of the sequence around rs5705329 (indicated by the red box) in 7 species. (TIF) Click here for additional data file. Figure S7 Analysis of the regulation of miR-146a expression by Ets-1 in HeLa cells. (A) Comparisons of the miR-146a promoter–reporter gene activity after cotransfection of an equal amount of an irrelevant carrier vector or an ETS1- or PBX1-expressing vector. The data shown are means ± SEM and are representative of three independent experiments performed in triplicate. *** indicates P<0.001. (B) Comparisons of the miR-146a promoter–reporter gene activity after the cotransfection of ETS1 siRNA (siETS1 #1 and siETS1 #2) or a negative control (NC). The data shown are means ± SEM and are representative of three independent experiments performed in triplicate. * indicates P<0.05. (C) Western blot analysis of Ets-1 levels after the transfection of the indicated expression vectors or siRNA. In the overexpression assay, the cells were collected 24 hours after transfection; in the siRNA-mediated knockdown assay, the cells were collected 48 hours after transfection. GAPDH was used as the loading control. (TIF) Click here for additional data file. Figure S8 Effect of ectopic Ets-1 expression on the activity of the allelic miR-146a promoter–reporter gene constructs. Reporter gene constructs containing the A or G miR-146a sequence were cotransfected into HeLa cells with different amounts of ETS1-expressing vector (0, 50, or 200 ng). For these three groups, 200 ng, 150 ng, or 0 ng of an irrelevant carrier vector was cotransfected, respectively, so that equal amounts of total plasmid DNA were used in all groups. The relative luciferase activity was measured 24 hours after transfection (upper). Cotransfection of a PBX1-expressing vector was used as the negative control. Also shown are the average G/A ratios of the luciferase activity of the allelic constructs (lower). (TIF) Click here for additional data file. Figure S9 Comparison of miR-146a expression levels in healthy individuals with different genotypes of rs2431697 or rs2431099. The horizontal line indicates the mean expression level within each group. (JPG) Click here for additional data file. Table S1 A list of the variants identified by the initial sequencing of the miR-146a region. (DOC) Click here for additional data file. Table S2 Association between the seven common SNPs around miR-146a and SLE. (DOC) Click here for additional data file. Table S3 Association between the rs57095329 G allele and lupus nephritis. (DOC) Click here for additional data file. Table S4 Analysis of OR in case-control groups carrying different numbers of risk alleles of either miR-146a or ETS1 SNP. (DOC) Click here for additional data file. Table S5 Conditional analysis of three SNPs in 5q33.3 in SLE cases and controls. (DOC) Click here for additional data file. Table S6 Haplotypic association of three SNPs in 5q33.3 with SLE. (DOC) Click here for additional data file. Table S7 A list of the primers used for the various assays. (DOC) Click here for additional data file.
