Identification of a new susceptibility variant for multiple sclerosis in OAS1 by population genetics analysis.

Contrasting results have been reported concerning the association of a splice-site polymorphism (rs10774671) in OAS1 with multiple sclerosis (MS). We analysed two OAS1 regions encompassing alternatively spliced exons. While the region carrying the splice-site variant is neutrally evolving, a signature of long-standing balancing selection was observed across an alternative exon 7. Analysis of variants in this exon identified an insertion/deletion polymorphism (rs11352835, A/-) that originates predicted products with distinct C termini. This variant is located along the major branch of the haplotype genealogy, suggesting that it may represent the selection target. A case/control study for MS indicated that rs11352835 is associated with disease susceptibility (for an allelic model with the deleted allele predisposing to MS, OR 1.27, 95% CI 1.072-1.513, p = 0.010). No association was found between rs10774671 and MS. As the two SNPs are in linkage disequilibrium in Europeans, the previously reported association between rs10774671 and MS susceptibility might be driven by rs11352835, possibly explaining the contrasting results previously observed for the splice-site polymorphism. Thus, we describe a novel susceptibility variant for MS in OAS1 and show that population genetic analyses can be instrumental to the identification of selection targets and, consequently, of functional polymorphisms with an effect on phenotypic traits.

Introduction

In humans four 2′,5′-oligoadenylate synthetase genes (OAS1, OAS2, OAS3 and OASL) are located on the long arm of chromosome 12 and play a central role in the innate immune response against viruses. These enzymes are activated either by the presence of double strand RNA or by single strand RNA with secondary structure, and catalyze the oligomerization of ATP into 2′,5′-linked oligoadenylates. These, in turn, bind to and activate the latent ribonuclease L (RNase L), which degrades viral and cellular RNA and blocks protein synthesis. In line with this function, human OAS1 and its murine ortholog, Oas1b, modulate the susceptibility to viral infections (Lim et al. 2009; He et al. 2006; Hamano et al. 2005; Mashimo et al. 2002; Perelygin et al. 2002). In particular, variants in human OAS1 have been associated with West Nile Virus and SARS infection (Lim et al. 2009; He et al. 2006; Hamano et al. 2005), and with liver fibrosis progression in HCV infected individuals (Li et al. 2009).

One of the most intensely studied variants in the human OAS1 gene is a G>A substitution (rs10774671) at the splice-acceptor site of exon 6. The A allele abrogates the production of OAS1 p46 isoform which displays high antiviral activity (Bonnevie-Nielsen et al. 2005) resulting in the production of the p48 and p52 isoforms (Fig. 1). Additionally, alternative splicing of OAS1 transcripts originates a p42 isoform which displays an alternative exon 5, and a p44 product deriving from skipping of exon 6 and inclusion of a downstream alternative exon (Fig. 1). A recent study indicated that the different protein products display variable activity in blocking Dengue virus replication via activation of RNase L (Lin et al. 2009).

Fig. 1

Schematic representation of the OAS1 gene region and alternative transcripts. Transcribed regions are shown in grey and the corresponding protein products are reported. The direction of transcription is indicated by the arrows. Alternatively spliced exons are shown in darker grey. The two regions we resequenced (OAS1 and OAS1 ) are indicated by the hatched lines. The asterisk shows the position of rs11352835; rs10774671 is circled. The LD (r 2) plot refers to CEU and data were derived from HapMap

In analogy to other genes involved in the antiviral response, polymorphisms and haplotypes in OAS1 have been identified as potential risk factors for autoimmune conditions, including type 1 diabetes (T1D) and multiple sclerosis (MS). Interestingly, viral agents have been proposed to trigger or exacerbate both diseases. The association between OAS1 polymorphisms and T1D is still controversial, with discordant results obtained in different studies (Qu et al. 2009; Smyth et al. 2006; Tessier et al. 2006; Field et al. 2005). In the case of MS, OAS1 SNPs and haplotypes have been associated with both disease susceptibility and severity (Fedetz et al. 2006; O’Brien et al. 2010). In particular, the A allele of rs10774671 (the splice-altering variant in exon 6), which originates protein isoforms with lower antiviral activity, was shown to confer both an increased risk to develop MS and a higher likelihood to have a more severe form of disease in a population from Ireland (O’Brien et al. 2010). Conversely, a previous study in a Spanish cohort found no association between this variant and MS susceptibility, although the authors found a rare haplotype to be more common in patients versus controls (Fedetz et al. 2006).

Antiviral response genes are an obvious target of natural selection as virus-borne infections have represented, and still represent, a major threat to human populations. Our previous data (Fumagalli et al. 2010) suggest that some overlap may exist among genes subjected to virus-driven selective pressure and loci involved in the pathogenesis of MS. Population genetic analyses of candidate gene regions can be regarded as instrumental to the identification of selection targets and, consequently, of functional polymorphisms with an effect on phenotypic traits. Here we analysed two gene regions in OAS1 and verified that while the region carrying the functional splice-site altering variant is neutrally evolving, long-standing balancing selection has shaped nucleotide diversity in a gene portion encompassing the alternative exon 7. A case/control association analysis indicated that variants in this downstream region, but not the splice-site altering variant, are associated with susceptibility to MS in an Italian cohort.

Materials and methods

HapMap samples and sequencing

Human genomic DNA from HapMap subjects (20 individuals for YRI, CEU and EAS) was obtained from the Coriell Institute for Medical Research. All analysed regions were PCR amplified and directly sequenced; primer sequences are available upon request. PCR products were treated with ExoSAP-IT (USB Corporation, Cleveland, OH, USA), directly sequenced on both strands with a Big Dye Terminator sequencing Kit (v3.1 Applied Biosystems) and run on an Applied Biosystems ABI 3130 XL Genetic Analyzer (Applied Biosystems). Sequences were assembled using AutoAssembler version 1.4.0 (Applied Biosystems), and inspected manually by two distinct operators. The genomic DNA of 3 Pan troglodytes was obtained from the Gene Bank of Primates, Primate Genetics, Germany (http://dpz.eu/index.php).

Patients, controls, and genotyping

For the case/control association study, 660 MS subjects (66% females) and 503 age- and sex-matched healthy controls (65% females) were recruited. All subjects were Italian of Caucasian origin. Median age was 42.5 (SD 12.1) and 43.9 (SD 20.1) years for MS and controls, respectively. Patients and controls were recruited at the MS Centre of Don Gnocchi Foundation in Milan and at Department of Neurological Sciences, University of Milan. All subjects gave informed consent according to protocols approved by the local Ethics Committees. MS patients underwent a standard battery of examinations, including medical history, physical and neurological examination, screening laboratory test, brain magnetic resonance imaging (MRI). All patients with MS fulfilled the McDonald’s criteria (McDonald et al. 2001).

Genotyping of rs11352835 and rs10774671 was performed by direct resequencing, as described above, using genomic DNA extracted from peripheral blood.

The power to detect genetic association was estimated using the QUANTO software (version 1.2; http://hydra.usc.edu/GxE/) under a multiplicative mode of inheritance, assuming a disease prevalence in the Italian population (non-insular) of 0.0007 (Pugliatti et al. 2006), and a type I error rate of 0.05.

Data retrieval and haplotype construction

Genotype data for 5 kb regions from 238 resequenced human genes were derived from the NIEHS (National Institute of Environmental Health Sciences) SNPs Program web site (http://egp.gs.washington.edu). In particular, we selected genes that had been resequenced in populations of defined ethnicity including CEU, YRI and EAS (NIEHS panel 2).

Haplotypes were inferred using PHASE version 2.1 (Stephens et al. 2001; Stephens and Scheet 2005). Haplotypes for individuals resequenced in this study are available as online resource (Online Table 1).

Linkage disequilibrium analyses were performed using the Haploview (v. 4.1) (Barrett et al. 2005). Data for LD analysis were derived from HapMap.

Statistical analysis

Tajima’s D (Tajima 1989), Fu and Li’s D* and F* (Fu and Li 1993) statistics, as well as diversity parameters θ W (Watterson 1975) and π (Nei and Li 1979) were calculated using libsequence (Thornton 2003). Calibrated coalescent simulations were performed using the cosi package (Schaffner et al. 2005) and its best-fit parameters for YRI, CEU, and EAS populations with 10,000 iterations. Coalescent simulations were conditioned on mutation rate and recombination rate. The maximum-likelihood-ratio HKA test was performed using the MLHKA software (Wright and Charlesworth 2004), as previously proposed (Fumagalli et al. 2009). Briefly, 16 reference loci were randomly selected among NIEHS loci shorter than 20 kb that have been resequenced in the 3 populations; the only criterion was that Tajima’s D did not suggest the action of natural selection (i.e. Tajima’s D is higher than the 5th and lower than the 95th percentiles in the distribution of NIEHS genes). In all analyses, the chimpanzee sequence was used as the out-group. All calculations were carried out in the R environment (Team R Development Core 2008).

Haplotype analysis and TMRCA calculation

Median-joining networks to infer haplotype genealogy were constructed using NETWORK 4.5 (Bandelt et al. 1999). Estimate of the time to the most common ancestor (TMRCA) was obtained using a phylogeny-based approach implemented in NETWORK 4.5 using a mutation rate based on the number of fixed differences between chimpanzee and humans. An additional TMRCA estimate derived from application of a maximum-likelihood coalescent method implemented in GENETREE (Griffiths and Tavare 1995; Griffiths and Tavare 1994). The method assumes an infinite-site model without recombination; therefore, haplotypes and sites that violate these assumptions need to be removed: for the analysis of OAS1 we removed 1 variant. Again, the mutation rate μ was obtained on the basis of the divergence between human and chimpanzee and under the assumption both that the species separation occurred 6 MY ago (Glazko and Nei 2003) and of a generation time of 25 years. The migration matrix was derived from previous estimated migration rates (Schaffner et al. 2005). Using this μ and θ maximum likelihood (θ ML), we estimated the effective population size parameter (N e) which resulted equal to 22,838. With these assumptions, the coalescence time, scaled in 2N e units, was converted into years. For the coalescence process, 106 simulations were performed.

Real-time PCR for the OAS1 p44 isoform

Whole blood was collected from 21 Italian healthy controls by venipuncture in Vacutainer tubes containing EDTA (Becton-Dickinson, NJ, USA), and peripheral blood mononuclear cells (PBMC) were separated on lymphocyte separation medium (Organon Teknica, Malvern, PA). Total RNA was extracted PBMCs by using the acid guanidinium thiocyanate–phenol–chloroform method. The RNA was dissolved in RNase-free water, and purified from genomic DNA with RNase-free DNase (RQ1 DNase, Promega, Madison, WI, USA). One microgram of RNA was reverse transcribed into first-strand cDNA in a 20-μl final volume containing 1 μM random hexanucleotide primers, 1 μM oligo dT and 200 U Moloney murine leukemia virus reverse transcriptase (Clontech, Palo Alto, CA, USA). cDNA quantification for OAS1 p44 and GAPDH was performed by real-time PCR (DNA Engine Opticon 2; MJ Research, Ramsey, USA). Reactions were performed using a SYBR Green PCR mix (RealMasterMix SYBR ROX, 5 PRIME) and the results were expressed as ΔΔC (where “C ” is the cycle threshold) and presented as ratios between the target gene and the GAPDH housekeeping mRNA. Primers were designed to specifically amplify the p44 isoform of OAS1 and were located on exon 5 and 7 (Fig. 1) (primer sequences available upon request).

Results

Nucleotide diversity and neutrality tests

In order to study the evolutionary pattern of OAS1 in human populations we resequenced three HapMap population samples (Yoruba, YRI; Europeans, CEU; East Asians, EAS) across two gene regions. Region 1 (thereafter referred to as OAS1 ) covers exons 5 and 6, that undergo alternative splicing events to originate p46, p42, p48 and p52 (Fig. 1). Region 2 (OAS1 ) is downstream exon 6 and covers the alternative exon 7, whose inclusion originates p44 (Fig. 1). The two regions are in strong linkage disequilibrium (LD) in CEU and EAS (Fig. 1; Online Fig. S1).

A total of 9 and 41 variants were identified in OAS1 (2.5 kb) and OAS1 (5.5 kb), respectively. In the second region we analysed, rs11352835 is an A insertion/deletion (A/−) polymorphism falling within the coding sequence of the alternative exon 7. Therefore, due to a reading frame shift, the variant potentially originates two C termini with different length and amino acidic composition. Analysis of genomic alignments and BLAST searches against resequenced primate genomes indicated that the A allele is ancestral and that sequences with homology to exon 7 can be identified in hominoid primates (Pan troglodytes, Gorilla gorilla, Pongo pygmaeus), Old World monkeys (Macaca mulatta, Papio hamadryas) and New World monkeys (Callithrix jacchus), but not in lemurs and non-primate mammals.

For the OAS1 and OAS1 regions we calculated nucleotide diversity by means of two indexes: θ W (Watterson 1975), an estimate of the expected per site heterozygosity, and π (Nei and Li 1979), the average number of pairwise sequence nucleotide differences. In order to compare the values we obtained for the two OAS1 regions, we calculated θ W and π for 5 kb windows (thereafter referred to as reference windows) deriving from 238 genes resequenced by the NIEHS program in the same population samples; the percentile rank corresponding to OAS1 and OAS1 in the distribution of NIEHS gene values is reported in Table 1 and indicates that OAS1 displays extremely high nucleotide diversity in both CEU and EAS; conversely, no exceptional values are observed for OAS1 , with the exclusion of YRI that shows reduced genetic variability in the region.

Table 1

Summary statistics for the OAS1 and OAS1 gene regions

Reg.	L a	P b	S c	θ W d		πe		Tajima’s D			Fu and Li’s D*			Fu and Li’s F*‘
				Value	Rankf	Value	Rankf	Value	Rankf	p g	Value	Rankf	p g	Value	Rankf	p g
OAS1 r1	2.5	YRI	4	3.72	0.025	2.49	0.020	−0.79	0.32	0.61	−1.10	0.22	0.75	−1.17	0.22	0.74
		CEU	4	3.72	0.28	7.62	0.75	2.49	0.98	<0.01	1.03	0.86	0.10	1.71	0.97	0.023
		EAS	6	5.59	0.59	6.80	0.67	0.58	0.67	0.32	1.19	0.95	0.08	1.17	0.89	0.12
OAS1 r2	5.5	YRI	11	4.73	0.13	3.74	0.15	−0.63	0.39	0.53	−0.19	0.55	0.38	−0.39	0.51	0.42
		CEU	40	17.20	0.98	31.10	0.99	2.84	0.98	<0.01	0.90	0.83	0.14	1.86	0.97	0.013
		EAS	37	15.91	0.99	21.92	0.98	1.32	0.87	0.13	1.64	>0.99	0.01	1.81	0.99	0.02

aLength of analyzed sequenced region (kb)

bPopulation (40 chromosomes were analysed in each population)

cNumber of segregating sites

dWatterson’s θ estimation per site (×10−4)

eNucleotide diversity per site (×10−4)

fPercentile rank relative to a distribution of 238 5 kb segments from NIEHS genes

g p value applying demographic coalescent simulations

High levels of nucleotide diversity are consistent with the action of balancing selection as linked polymorphisms are maintained together with the selected variant(s). Another effect of balancing selection is a distortion of the site frequency spectrum (SFS) towards intermediate frequency alleles.

Common neutrality tests based on the SFS include Tajima’s D (D T) (Tajima 1989) and Fu and Li’s D* and F* (Fu and Li 1993). D T tests the departure from neutrality by comparing θ W and π and positive values indicate an excess of intermediate frequency variants. Fu and Li’s F* and D* are also based on SNP frequency spectra and differ from D T in that they also take into account whether mutations occur in external or internal branches of a genealogy. Since, population history, in addition to selective processes, is known to affect the SFS, the significance of neutrality tests was evaluated by performing coalescent simulations with population genetics models that incorporate demographic scenarios (see “Materials and methods”). As above, we also applied an empirical comparison by calculating the percentile rank of D T, F* and D* in the OAS1 regions relative to 5 kb reference windows. Neutrality tests for OAS1 region indicated departure from neutrality with significantly positive values for most statistics in CEU and EAS (Table 1). In line with these findings, D T, as well as Fu and Li’s F* and D* calculated for OAS1 rank above the 95th percentile of the distribution of 5 kb reference windows in these two populations. Conversely, no departure from neutrality was observed for OAS1 (Table 1) with the exception of CEU where D T and F* were significantly high. Yet these values are calculated over a very small number of segregating sites (n = 4) and are likely due to LD with OAS1 .

Overall, these results suggest that nucleotide diversity at OAS1 has been shaped by balancing selection in European and Asian populations, with the OAS1 region being neutrally evolving.

As mentioned above, a hallmark of balancing selection is an excess of polymorphism compared to neutral expectations. Indeed, our data (Table 1) indicate that nucleotide diversity indexes are extremely high for OAS1 in CEU and EAS. Yet, polymorphism level also depends on local mutation rates; therefore, under neutral evolution, the amount of within- and between-species diversity is expected to be similar at all loci in the genome (Kimura 1983). The multi-locus HKA test was developed to verify this expectation (Wright and Charlesworth 2004). We applied a multi-locus MLHKA (maximum-likelihood HKA) test by comparing polymorphism and divergence levels at the OAS1 r1 and OAS1 r2 genomic regions with 16 NIEHS genes resequenced in YRI, CEU and EAS. As shown in Table 2, a significant excess of polymorphism compared to divergence is observed for CEU and EAS in OAS1 , while no deviation from expectations is detectable at OAS1 .

Table 2

MLHKA test for the two gene regions we analysed

Region	Pop.a	MLHKA
		k b	p
OAS1 r1	YRI	0.46	0.48
	CEU	1.16	0.74
	EAS	1.28	0.31
OAS1 r2	YRI	1.08	0.18
	CEU	6.11	3.73 × 10−6
	EAS	7.18	5.42 × 10−6

aPopulation

bSelection parameter (k > 1 indicates an excess of polymorphism relative to divergence)

Haplotype analysis

Further insight into the evolutionary history of a gene region can be gained by inferring haplotype genealogies. This has both a descriptive purpose (i.e. showing the relationship among alleles and their distribution in human populations) and can be used to test for selection. In particular, balancing selection is expected to result in two or more major haplotype clades with a deep coalescence time. Here we constructed haplotype genealogies using two approaches: a neighbour-joining network and a maximum-likelihood coalescent method implemented in GENETREE.

Consistently with the extended LD pattern (Fig. 1; Online Fig. S1), the haplotype network for the OAS1 gene region presented no reticulations or recurrent mutations (Fig. 2). The haplotype genealogy is split into two major clades in CEU and EAS, while all African chromosomes cluster within clade B. In line with these results, calculation of population genetic differentiation (F ST) between CEU and YRI resulted in a value of 0.366, corresponding to a percentile rank of 0.96 in the distribution of F ST values calculated for 5 kb reference windows.

Fig. 2

Network analysis for OAS1 . Each node represents a different haplotype, with the size of the circle proportional to frequency. Nucleotide differences between haplotypes are indicated on the branches of the network. The position of rs11352835 (A insertion/deletion in the alternative exon 7) is shown. Circles are colour-coded according to population (green YRI, blue CEU, red EAS). The most recent common ancestor (MRCA) is also shown (black circle). The relative position of mutations along a branch is arbitrary

The absence of African chromosomes in clade A is not merely due to the relatively small sample of individuals we analysed, as HapMap data from 120 YRI chromosomes indicated that several SNPs along the major branches are monomorphic in this population.

In order to estimate the time to the most recent common ancestor (TMRCA) of the haplotype genealogy, we applied a phylogeny-based method (Bandelt et al. 1999). For this analysis only single base pair substitutions were included (i.e. rs11352835, A ins/del was excluded). Using a mutation rate based on the number of fixed differences with chimpanzee and a separation time of 6 million years (MY) (Glazko and Nei 2003), we estimated a TMRCA of 7.4 MY (SD 1.36 MY). Consistent with these results, GENETREE analysis (Fig. 3) yielded a TMRCA estimate of 5.61 MY (SD 0.88 MY). These coalesce times are much deeper than estimated under neutrality.

Fig. 3

GENETREE analysis for OAS1 . Mutations are represented as black dots and named for their physical position along the region. The absolute frequency of each haplotype is also reported. Note that mutation numbering does not correspond to that reported in Fig. 2

Given the deep coalescence time of OAS1 haplotypes, we resequenced corresponding region in eight chimpanzees (Pan troglodytes) but no variants shared with humans was identified.

Finally, we verified that SNPs defining the two major haplotype clades do not affect (at least overtly) the inclusion in the transcript of the alternative exon 7: a PCR amplification of cDNA derived from lymphoblastoid cell lines of four subject homozygous for A clade haplotypes and four homozygous for clade B revealed similar levels of exon 7 inclusion (not shown). In order to further address this issue, we performed real-time PCR experiments to analyse the level of p44 expression in PBMCs derived from 21 healthy volunteers. No difference in p44 expression levels was observed among individuals with different genotype at rs11352835 (Kruskal–Wallis test, 2df, p = 0.54) (Online Fig. S2).

OAS1 polymorphisms in MS

As mentioned above, OAS1 has been recently described as a multiple sclerosis susceptibility gene (Fedetz et al. 2006; O’Brien et al. 2010). We therefore analysed a population of 660 subjects suffering from MS and 503 sex-matched healthy controls (HC). All subjects were Italian of Caucasian origin and were genotyped for two SNPs: rs10774671 (located at the splice-acceptor site of exon 6) and rs11352835 (the A ins/del polymorphism in the alternative exon 7) (Fig. 1). This latter was selected because of its putative functional role and as a consequence of its location in the haplotype phylogeny. Indeed, as shown in Fig. 2, rs11352835 is located on the basal branch of the haplotype network, suggesting that it may represent (or be in full LD with) the selection target. The two SNPs we selected for genotype analysis are in relatively tight LD (r 2 = 0.87) and both complied to Hardy–Weinberg equilibrium in both cases and controls. We found the allele frequency of both SNPs to be significantly different in MS compared to HC (Table 3), but only rs11352835 withstood Bonferroni correction for multiple tests (Fisher’s exact test, Bonferroni corrected p = 0.010, OR 1.27, 95% CI 1.072–1.513). Comparison of genotype frequencies in the MS and HC samples indicated no significant difference for rs10774671 (Table 3). Conversely, the genotype distribution of rs11352835 was significantly different in the two cohorts (Bonferroni corrected p = 0.030, Table 3). The odds ratio for a recessive model with the deleted allele predisposing to MS was 1.41 (95% CI 1.103–1.811, p = 0.0051). Haplotype analysis identified no significant association.

Table 3

Association of OAS1 SNPs with susceptibility to MS

SNP	Alleles/genotypes	Counts		p a (corrected)
		MS (n = 660)	HC (n = 503)
Alleles
rs10774671	G/A	488/832	416/590	n.s.
rs11352835	A/−	469/851	415/591	0.010
Genotypes
rs10774671	GG/AG/AA	95/298/267	87/242/174	n.s.
rs11352835	AA/A−/−−	87/295/278	83/249/171	0.030

n.s. not significant

aBonferroni corrected p values (p values were corrected for two tested SNPs)

Discussion

Most human genes undergo extensive alternative splicing, which is generally regarded as a powerful mechanism to increase proteome diversity (Modrek and Lee 2002). This is likely even more true for immune response genes, which are more frequently engaged in alternative splicing events compared to other functional classes (Modrek et al. 2001). Several studies have addressed the evolutionary forces underlying the appearance and fate of alternatively spliced exons in inter-specific comparisons (Xing and Lee 2005; Lu et al. 2009a, b; Ke et al. 2008; Zhang and Chasin 2006; Plass and Eyras 2006; Chen and Chuang 2006; Ermakova et al. 2006). Conversely, fewer analyses have focused on the selective patterns shaping the intra-specific diversity of alternatively spliced genes. Here we selected for analysis the two gene regions in OAS1 that encompass alternatively spliced exons. As mentioned above, at least five transcripts originate from OAS1, one of these is genetically determined by the presence of a splice-acceptor variant and gives rise to protein products with different antiviral activity (Bonnevie-Nielsen et al. 2005; Lin et al. 2009). The splice-acceptor polymorphism has also been associated with altered susceptibility to viral infections (Lim et al. 2009), indicating this region as a possible selection target. Surprisingly, our resequencing analysis in three human populations showed that OAS1 is neutrally evolving, as all tests failed to reject the null hypothesis of selective neutrality. Low diversity indexes were observed for this region in YRI, but the MLHKA test revealed that polymorphism levels were not exceptionally low compared to inter-specific diversity; moreover, calculation of Fay and Wu’s H (Fay and Wu 2000), which detects an excess of high-frequency derived alleles, in both OAS1 and OAS1 yielded results consistent with selective neutrality in YRI (not shown), suggesting that weak negative selection (as opposed to positive selection) may be acting on OAS1 in African populations.

Conversely, data herein indicate that the region encompassing the alternatively spliced exon 7 has been a target of long-standing balancing selection in Europeans and Asians. Inclusion of this exon in the OAS1 transcript originates a p44 isoform (Fig. 1) whose antiviral activity has been poorly characterized. A recent report indicated that, among all OAS1 isoforms, only p42 and p46 display strong antiviral activity against Dengue virus (Lin et al. 2009). The reason for the differential activity of OAS1 proteins is presently unknown, as all isoforms display the catalytic OAS domain (which is located at the N terminus). Our analysis indicated that an A insertion/deletion polymorphism (A/-, rs11352835) is located on the major branch that separates the two haplotype clades, indicating that this variant, or one in linkage to it, may represent the selection target in OAS1 . The derived deleted allele determines a frame shift that is predicted to result in the inclusion of 22 additional amino acid residues compared to p44. Given the differential activity observed against Dengue virus for p42 and p44 (which only differ by 35 terminal amino acids) (Lin et al. 2009), it is tempting to speculate that similarly, such a minor change in sequence composition may alter the antiviral potential of OAS1. One interesting possibility is that the different composition of the C terminus affects the specificity of OAS1 against different viral species. An alternative explanation for our findings is that the selection target is accounted for by a regulatory variant/haplotype with a role in transcription modulation, as several expression QTLs for OAS1 have been identified in the gene region encompassing exon 7 (http://eqtl.uchicago.edu/).

Although further analyses will be required to identify the selection target in OAS1 , our data provide evidence that a functional polymorphism is located within this region, as signatures of long-standing balancing selection are expected to extend over relatively short genomic regions (Charlesworth 2006). Thus, we reasoned that a susceptibility variant for MS might be located in this region, as well, possibly explaining the previous contrasting results obtained for the splice-acceptor SNP in different MS cohorts. Indeed, while a strong relationship was recently described between the G allele of rs10774671 and protection from MS in a population from Ireland (O’Brien et al. 2010), a previous report in a Spanish sample revealed no such association, despite similar cohort sizes in the two studies (about 400 cases/controls) (Fedetz et al. 2006). Our results on an Italian sample indicated that the A allele of rs11352835, but not the splice-site variant, is significantly associated with protection from MS. Given that the two variants have very similar frequency in the Italian sample, the estimated power to detect association is also comparable, and reaches 78% assuming an OR of 1.27 and a sample size as the one we analysed. The relatively tight LD of the two SNPs (r  = 0.87) suggests that association between rs10774671 and MS susceptibility is driven by rs11352835. Therefore, minor differences in the level of LD in cohorts with distinct geographic origin might account for the contrasting results previously observed for the splice-site polymorphism. In this respect, it is also worth mentioning that insufficient sample size is regarded as a major explanation for several inconsistent associations reported in the scientific literature. This consideration may apply to previous studies on the role of OAS1 in MS, which relied on relatively limited patient/control samples (Fedetz et al. 2006; O’Brien et al. 2010). Here we genotyped larger MS/HC cohorts compared to previous studies; our data confirm that OAS1 represents an MS susceptibility gene, and refine the association signal to a variant different from, but in linkage with, the one previously described to associate with disease.

In summary, our population genetic analysis indicates that the OAS1 gene region harbouring the functional splice-acceptor polymorphism is neutrally evolving, and the variant itself shows no association with MS in Italians. Conversely, a new MS susceptibility variant was identified within or in proximity to the alternatively spliced exon 7, a region targeted by balancing selection. In line with previous reports (Grossman et al. 2010), our work demonstrates that population genetic approaches can be regarded as instrumental to the identification of selection targets and, consequently, of functional polymorphisms with an effect on phenotypic traits. This observation suggests that the ever-increasing availability of resequencing data, obviously including the 1000 Genomes Project (http://browser.1000genomes.org/), can be exploited to identify selection signatures at the genome-wide level and this information, in turn, may be used to prioritize variants to be typed in association analysis for complex diseases and traits.

Below is the link to the electronic supplementary material.

Supplementary Tables (DOC 696 kb)

Supplementary Figures (DOC 71 kb)