A genome-wide survey for SNPs altering microRNA seed sites identifies functional candidates in GWAS

Approximately 5% of validated MREs contain SNPs in their seed site

No population frequency data are available for 29% of the MRE SNPs (9 of 31), a value that will change as whole genome sequence data from more individuals surface and with completion of more encompassing GWAS studies. Nine SNPs have minor allele frequencies (MAF) < = 2%, and may be considered rare in the general population and therefore unlikely as common factors in complex disease. The 13 remaining SNPs have allele frequencies above 2% in at least one population listed in dbSNP (Table 1). Of note, 7 of the 31 (22%) SNPS identified have shown association with disease traits, emphasizing the potential importance of MRESS SNPs as modulators of disease risk.

Previous studies have estimated MRESS SNP density to be lower than that observed in regions outside the MRESS, suggesting a higher rate of negative selection on MRESS [3, 10]. In light of the updated account of variation in the human genome available in dbSNP build132, we estimated the frequency of SNPs falling within MRESS and those falling outside of MRESSs using the 606 sequences from validated target sites identified from the above analysis. We performed a sliding window search of 6 bases, (the size of seed positions 2-7) starting 18 bases upstream and continuing to 24 bases downstream of each validated MRE site, sliding at a 1 base step. The 0 mark of the x-axis in Figure 1 demarcates the second position of the MRESS (or first position of the 6mer seed site). Our data indicate that the MRESS contains the lowest amount of variation across the region, an observation in agreement with prior analyses [3, 10]. Although there is less variation across the MRESSs there appears a considerable density of SNPs (5.5/kb) falling within what is thought to be the most important sequence for miR-mRNA interactions.

Genome wide survey of predicted MRESS and CNM SNPs

Although many variants have been associated with the modulation of phenotypes relevant to disease in GWAS studies, the challenge of determining which of them may be casual remains [20]. To identify potential causal variants mediated by MRESS creation or disruption we first performed a genome-wide survey for SNPs falling within computationally predicted conserved MRESSs. To do this, we utilized the http://microRNA.org portal to access a collection of predicted miR-mRNA interactions. These predictions were derived using an algorithm that incorporates an array of the most recent miR prediction guidelines, such as seed-site pairing, site context, free-energy, and target conservation across multiple vertebrates [21]. We focused our analysis on MRESSs conserved across mammals as these are more likely to be of functional significance (see methods) [16]. By comparing the gDNA coordinates of each predicted MRESS against the gDNA coordinates of dbSNP132 SNPs, we identified 2723 MRESS SNPs interrupting 5797 conserved predicted interactions. To further prioritize these hits, we classified them by the type of seed match the MRESS SNP was predicted to interrupt; 8mer (2245), 7mer-8m (3251), 7-mer-A1 (180) or 6mer (121). Although there is overlap in the degree of efficiency of repression by these different seed type classes (likely dependent of site sequence context), there remains a hierarchy with 8mer sites being most efficient [1]. Interestingly 38% (2245) of MRESS SNPs fall within predicted 8mer MRESSs. It has been estimated that ~50% of predicted MREs are potentially functional and it is likely that a portion of the SNPs identified here fall within bona fide MREs [22]. Overall, we estimate that 3% of high confidence predicted conserved MRESSs contain SNPs.

In addition to SNPs that may interrupt MRESSs, SNP alleles may also create MRESSs. To identify potential CNM SNPs we performed a genome-wide computational survey for predicted MRESSs that are created when the mRNA sequence contains the non-reference allele of hg19. Using the Ensembl variation API tools we retrieved the flanking 22 nt sequence from both upstream and downstream of the non-reference allele of every 3'UTR dbSNP132 SNP. Each sequence containing the non-reference allele was analyzed for potential miR-target sites using the miRanda software [23]. This analysis provided us with 22295 CNM SNP creating 49047 miR-mRNA predictions which where also categorized by the seed-type they created; 8mer (10333), 7mer-8m (36188), and 7-mer-A1 (526) or 6mer (2000). It should be noted that there are many more predictions for CNM than MRESS SNPs. This is due to the fact that no conservation constraint was imposed on the CNM SNP predictions (See methods). Considering that many CNM SNPs presumably arise to create new regulatory sites, filtering our hits on conservation status would be counter-intuitive. Of note, we found that approximately 28% (6946/25018) of predicted MRESS and CNM SNPs identified here were first identified by the 1000 Genomes project.

Some SNPs identified in GWAS are in LD with predicted conserved MRESS and CNM SNPs

GWAS have been a powerful approach to identify genetic variants that contribute to disease risk. However, a functional role for many of the SNPs identified has not been elucidated. It is likely some of these SNPs are in strong LD with unknown functional ones, some of which could be among the predicted conserved MRESS and CNM SNPs. To investigate this possibility, we searched the resulting MRESS and CNM SNP data for variants in LD with SNPs showing association, of GWAS significance, with disease traits and related phenotypes. To do so, we retrieved a dataset of 4817 reported associations, collected from GWAS studies, between 3943 unique SNPs and disease traits shown to have P-values meeting a threshold of < 1.0 × 10-5 [24]. These data were processed through SNAP http://www.broadinstitute.org/mpg/snap/ldsearch.php which yielded a list of SNPs (including those from the pilot 1000 Genomes Project data) in LD with those reported in the GWAS. We limited our search to an r² of > 0.8 for the CEU population. The results of this query were searched against both the MRESS SNP and the CNM SNP predictions. This query identified 35 instances of an MRESS SNP (some MRESS SNPs are in LD with more than one reported GWAS SNP associating with multiple phenotypes), in LD (r² > 0.8) with a least one reported GWAS SNP or an original GWAS SNP, associating with disease phenotypes (Additional file 1). In total there were 14 MRESS SNPs in 11 genes associating with 16 traits. We also identified 124 instances of a GWAS SNP that associates with disease traits and is in LD (r² > 0.8) with CNM SNPs (Additional file 2). There were 73 CNM SNPs in 73 genes associating with 52 traits.

In total we identified 87 SNPs (14 MRESS and 73 CNM SNPs) in very strong LD with SNPs reported as associating with disease related phenotypes. These 87 SNPs represent 2.22% of the total unique SNPs reported in the GWAS data. Using the SNAP pairwise LD tool we determined that 6 of the 3943 GWAS SNPs are in LD with each other, giving us 3940 SNPs or regions associating with disease traits. To determine the possibility of this number occurring by chance we first filtered dbSNP 132 for SNPs having a minor allele frequency (MAF) > = 1%, which was the lowest MAF reported in the GWAS data. We next selected randomly 3940 SNPs and, as we did with the GWAS data, ran them through SNAP to determine the SNPs in LD. From this list we determined the number of SNPs found in our MRESS and CNM SNP data, and this analysis was repeated 1000 times. From these 1000 simulations the mean number of SNPs found in the MRESS and CNM SNP data was 24.99 and the standard deviation 6.19. The probability of finding 87 SNPs by chance based on this distribution was calculated to be 1.08 × 10^-23. These SNPs may be considered the likely putative functional variants which represent proxy SNPs identified by GWAS.

Co-expression data identify functional candidates

eQTL data support several MRE target predictions when a MRESS or CNM SNP is present

Variation in gene transcript levels is thought to be an important modulator of disease risk in humans and SNPs that may mediate this variation are thought to be of great functional significance [26]. To investigate the contribution of SNPs associating with disease traits to transcript level variation, a number of expression Quantitative Trait Loci (eQTL) studies have been performed [27–29]. These results have demonstrated a number of SNPs associating significantly with expression differences across collected tissue samples. Importantly, these studies have noted differences in the amount of transcript variation across tissue samples, suggesting SNPs may modulate regulatory mechanisms, in some cases, in a tissue specific manner [30]. Interestingly, a recent study has estimated that > 80% of miRs act to lower mRNA levels demonstrating mRNA destabilization is the primary mode of action of miRs on target mRNAs [31].

To determine if the 39 miR predictions with co-expression data identified in the previous section are supported by eQTL data we utilized the Genevar eQTL database web tool [32]. Genevar allows for querying and visualization of eQTL data for loci of interest using data from various studies. We utilized the results from the recent MuTHER study which reports eQTL data from twin pairs in 3 tissue types; 78 twin-pair lymphoblastoid cell line (LCLs) biopsies, 80 twin-pair skin cell biopsies and 83 twin-pair fat cell biopsies [29]. Searching for eQTL data on each of the 39 MRESS and CNM SNPs we found 11 of the 39 had genotype specific transcript level data in at least one of the 3 tissue samples investigated in the MuTHER study for which there was also evidence of co-expression in mimiR for this tissue. Four of these 11 SNPs showed marginally significant trends in the differences in transcript levels across genotypes (Figure 2).

SYS1 transcript levels were shown to be significantly different among rs2245717 genotypes from LCL in one of the two twin study groups (Figure 2). Although, the second twin group failed to achieve significance, the direction of the effect was in agreement with the first group. Furthermore, the lower transcript levels associate with the allele predicted to create the miR-150 binding site. Neither of the two twin adipose tissue sample groups for GFOD2 transcript levels showed a significant difference among rs12449157 genotypes. Interestingly, both p-values are of nominal significance and the direction of the effect supports the predicted miR interaction and subsequent effect of the CNM SNP (Figure 2). IKZF3 transcript levels measured in LCL cells showed significant differences among the rs907091 genotypes in both twin groups. Lower IKZF3 levels were observed in carriers of the G allele which is predicted to create a miR-326 MRESS. However, a second probe found on the Illumina whole genome expression array, used in the study, shows conflicting data where there is no difference between transcript levels in either group. The rs3810291 SNP in the ZC3H4 3'UTR shows no significant difference in transcript levels between alleles in adipose samples used in the MuTHER study. However, literature mining for "eQTL, " and the corresponding gene and phenotype identified an additional study showing eQTL data that supports an allelic difference, in the correct direction, and in adipose tissue for rs3810291 [33]. Taken together, this information suggests these four SNPs may have functional significance.

MRESS and CNM SNPs in Positive Selection

Genetic variants that have been subject to selection are most likely the functional variants [34, 35]. The fixation index (F_ST) statistic measures population differentiation and provides a test for the influence of selective pressures, where higher F_ST values indicate local positive adaptation and lower values negative or neutral selection [36]. As adaptive genetic variants have been driven to higher frequencies by environmental factors (i.e, positive selection), SNPs showing high F_ST values may be considered high priority candidates for association studies for gene by environment studies. Such variants also may play a role in the observed variation and potentially influence disease prevalence across populations [37, 38]. To determine if the identified 3'UTR SNPs that create or disrupt predicted MRESSs may be under positive selection we first downloaded genome wide F_ST calculations for HapMap Phase 3 data [39]. We found a significant difference (P = 0.0004) between the mean transformed F_ST values of combined MRESS and CNM SNPs (n = 2448) and the remaining (i.e. non-MRESS or non-CNM) 3'UTR SNPs (n = 19906) for which F_ST data were available. Figure 3 shows the number of F_ST values between MRESS and CNM SNPs and non-MRE SNPs across 10 F_ST bins. As F_ST values increase there is a clear increase in MRESS and CNM SNPs compared to the remainder. This observation further supports that MRESS and CNM SNPs are likely functional variants.

Retrieval and use of dbSNP information

We retrieved all dbSNP build 132 SNP (as of 11-31-10) information by downloading the vcf file available through the 1000 genomes home page [18]. To ensure we only surveyed variation in the form of SNPs and not indels and/or copy number variants we removed all SNPs not reported as bi-allelic. A subset of data containing all 3'UTR SNPs (n = 210042) was extracted using Perl. This dataset was used for all subsequent analyses. To determine the percentage of SNPs that were submitted to dbSNP by the 1000 Genomes Project we used UCSC Genome Browser to identify all dbSNP build 132 SNPs where the submitter status handle was equal to only the 1000GENOMES tag. This data set was then searched against our MRESS and CNM data sets to identify those SNPs contributed by the 1000 Genomes Project. To retrieve allele frequency data on SNPs reported in MRESSs and CNMs, we utilized the Perl API variation tools accessing the latest human genome variation data, build 61. All data analysis was performed on the NUGO information network [51].

Identifying MRE SNPs in validated targets

miRecords hosts a collection of validated miR-mRNA interactions built from an exhaustive literature search and the database of records was download in a tab-delimited format [19]. We next annotated each hit for target site functionality, by checking the literature source for evidence of a loss of function experiment, which provided us with 606 validated MRE targets. To identify SNPs falling within these 606 validated targets sites, a Perl script was written to search each SNP gDNA coordinate against the gDNA coordinates of each target transcript MRESS, retrieved from Ensembl.

SNP Density determination

To calculate SNP density, a Perl script was written to perform a sliding window search (W = 6 bases as this corresponds in size to bases 2-7 of the MRESS) of the 606 validated MREs for SNPs, starting 18 bases upstream of the first base of the MRESS and ending 18 bases downstream of the last base of the MRESS. We report the number of SNPs for each position of the window across this sequence. Values are reported as SNPs/kb.

Identification of conserved MRESS and CNM SNPs

To identify SNPs falling in conserved MRESSs, we downloaded the "good mirsrv_score conserved miRNA" datafile from the http://microrna.org website. This file contains all predicted mRNA target motifs for targeting microRNAs which belong to conserved microRNA families. Conservation signal is used to predict functional MREs, however it has been determined that a conservation signal above background for MREs of the most recent mammalian specific miRNA families (non-conserved) was unlikely due to the relatively short time between the emergence of these miRs and the occurrence of new MREs within 3 'UTRs [16]. Therefore, to eliminate false positives that would arise from this form of analysis we utilize predictions for only conserved miR families - which are contained in the "good mirsrv_score conserved miRNA" datafile. To add an additional measure of conservation, we implemented a conservation score cutoff for predicted miR targets using a Phastcon score of > = 0.57, which authors of the work estimate corresponds to conservation across the mammalian genomes used in their study [21, 52, 53]. The Phastcon scores are provided in the predictions data file from http://microRNA.org. Additionally, predictions generated by the http://microRNA.org tool provide a ranking score (mirSVR score) which is calibrated to correlate linearly with the extent of down regulation of a miR on its target. Importantly, these scores may be interpreted as an empirical probability of down regulation. From these data we selected a mirsvr score cutoff < = -0.6, representing the top 12% of all predictions. A Perl script was then used to compare the gDNA coordinates of each predicted MRESS (n = 197287) against the gDNA coordinates of each 3'UTR SNP in dbSNP132.

To identify CNM SNPs, we utilized the Ensembl variation Perl API tools (Build 61) to retrieve the 22 bases flanking the 5' and 3' regions of every 3'UTR SNP in dbSNP 132. We generated the reverse complement for those mRNAs transcribed from the negative strand. These data were then run locally through the miRanda target prediction algorithm. To limit identification of potential false positives we implemented an arbitrary paring score cutoff of > = 150 and an energy cutoff of < = -20. We identified all predicted MRESS created by CNM SNPs by filtering hits on the position of their target prediction on the mRNA, where each SNP is located at position 23 of 45.

Retrieval of GWAS results and LD calculations

We first downloaded a catalog of human variation associating with disease phenotypes (1-25-11) [24]. The list was then submitted to the SNAP tool http://www.broadinstitute.org/mpg/snap/ldsearch.php using r² of > 0.8 for the CEU population. The resulting list was then searched against each of the CNM and MRESS SNP lists to identify 87 MRESS and CNM SNPs. To determine the probability of observing this number by chance, 3940 SNPs were randomly selected from the dbSNP build 132 data set filtered for SNPs with a MAF > = 1% and run through the SNAP tool to identify all SNPS in LD, using r² of > 0.8 for the CEU population. This procedure was repeated 1000 times. The probability of observing 87 MRESS or CNM SNPs randomly from the genome was determined based on the normal distribution generated from the 1000 simulations.

F_STcalculations

Genotype characteristics of 11 HapMap Phase 3 populations were split into 4 groups of similar ancestry; Asian, African, European and American. F_ST values where calculated for each HapMap Phase 3 SNP between these 4 groups and reported in a downloadable file [39]. Using a Perl script we extracted all MRESS and CNM SNPs with F_ST values from this dataset. We used the statistical analysis software (SAS) boxcoxar macro to transform the F_ST data to fit a normal distribution. We then performed an unpaired Student's t-test using the transformed values for these two groups to determine if they were significantly different. To identify MRESS and CNM F_ST outliers we selected those SNPs with F_ST > 2 SDs from the mean.

Authors have no competing interests to declare.