Associations between Single Nucleotide Polymorphisms in Cellular Viral Receptors and Attachment Factor-Related Genes and Humoral Immunity to Rubella Vaccination

Background Viral attachment and cell entry host factors are important for viral replication, pathogenesis, and the generation and sustenance of immune responses after infection and/or vaccination, and are plausible genetic regulators of vaccine-induced immunity. Methods Using a tag-SNP approach in candidate gene study, we assessed the role of selected cell surface receptor genes, attachment factor-related genes, along with other immune genes in the genetic control of immune response variations after live rubella vaccination in two independent study cohorts. Results Our analysis revealed evidence for multiple associations between genetic variants in the PVR, PVRL2, CD209/DC-SIGN, RARB, MOG, IL6 and other immune function-related genes and rubella-specific neutralizing antibodies after vaccination (meta p-value <0.05). Conclusion Our results indicate that multiple SNPs from genes involved in cell adhesion, viral attachment, and viral entry, as well as others in genes involved in signaling and/or immune response regulation, play a role in modulating humoral immune responses following live rubella vaccination.

In the current candidate gene association study, we follow up on previously found genetic associations and also make use of recent major discoveries in the virology field, and explore the plausible role of selected cell surface receptor-, and attachment factor-related genes, such as MOG and poliovirus receptor-related gene family members, in the genetic control of immune response variations after live rubella vaccination. Our results from two independent study cohorts (discovery and replication) strongly suggest that genetic variants from these genes play a role in modulating humoral immune responses following rubella vaccination.

Study Participants
The study cohort was a large population-based sample of 2,250 healthy children, older adolescents, and healthy adults (age 11 to 40 years), residing in Rochester, MN, and San Diego, CA, with clinical and demographic characteristics previously reported [7,9,21,28,29,30].
The Rochester cohort comprised a sample of 1,145 individuals from three independent age-stratified random cohorts of healthy schoolchildren and young adults from all socioeconomic strata from Olmsted County, MN, enrolled between 2001 and 2009, as published elsewhere [7,9,21,28,29,30]. Eleven hundred and one parents agreed to let their children join the current rubella vaccine study.
Between July 2005 and September 2006, we enrolled an additional 1,076 healthy older adolescents and healthy adults (age 18 to 40 years, San Diego cohort) from armed forces personnel at the Naval Health Research Center (NHRC) in San Diego, CA, as previously described [30]. All subjects included in the current rubella vaccine study had a documented receipt of measlesmumps-rubella (MMR) vaccine. The Institutional Review Boards of the Mayo Clinic and the NHRC approved the study, and written informed consent was obtained from each subject, from the parents of all children who participated in the study, as well as written assent from age-appropriate participants.

Rubella Virus-specific Neutralization Assay (sICNA)
We used a modified version of the immunocolorimetric-based neutralization method described by Chen et al. [31], optimized to a high-throughput micro-format, to measure rubella virus-specific neutralizing antibodies [25,26]. Subjects' sera were heat-inactivated for 1 hour at 56uC. Sera were serially diluted in two-fold, in triplicate for each dilution, beginning from 1:12.5 through 1:100 (to yield a final volume of 30 mL per dilution), using phosphatebuffered saline (PBS, pH 7.4) supplemented with 1% fetal bovine serum (FBS). Rubella virus stock (vaccine virus HPV77) was diluted to a working concentration of 1.2610 3 plaque-forming units (PFU)/mL, and was added (30 mL) to an equal volume of diluted serum (or diluent as in the case of virus-only control), yielding a final serum dilution series of 1:25 through 1:200. The plate was incubated for 1.5 hour at 37uC, 5% CO 2 . Fifty microliters of each mixture were used to inoculate confluent Vero cell monolayers (cultured in flat-bottom 96-well plates) and the cells were incubated for 1 hour at 37uC, 5% CO 2 . After the incubation period, DMEM supplemented with 5% FBS and 50 mg/mL Gentamicin (Gibco; Invitrogen, CA, USA) was added to each well and the plate was further incubated for 72 hours at 37uC, 5% CO 2 . The plates were developed using an indirect immunocolorimetric method and rubella E1-specific monoclonal antibody (Centers for Disease Control and Prevention/CDC, GA, USA), as previously described [31]. Optical density (OD) values were measured using a measurement/reference wavelength pair of 450 nm/630 nm. The neutralization titer was calculated as the highest dilution at which the input virus signal was reduced by at least 50% within the dilution series (neutralization titer 50/NT 50 ). The Loess method of statistical interpolation was used to estimate interpolated NT 50 values from observed values [25,32]. The intraclass correlation coefficient (ICC), based on log-transformed estimates from samples with repeated NT 50 measurements was, 0.89, demonstrating a high degree of reproducibility in the assay [25,26].

SNP Selection and Candidate SNP Genotyping
The SNP selection and genotyping methods described herein are similar or identical to those published for our previous genetic association studies [7,8,9,28,29]. The genotyping effort included follow-up on 197 previously identified genetic associations (p, 0.05) between SNPs in immune response genes (cytokine and cytokine receptor genes, Toll-like receptor and signaling genes, vitamin D and vitamin A receptor family genes, antiviral effector and other innate genes) and immune measures after rubella vaccination. A total of 571 additional SNPs on the GoldenGate 768-plex panel were filled with candidate tagSNPs selected from eight newly discovered viral receptor genes (MOG, PVR, PVRL1, PVRL2, PVRL3, PVRL4, BTN2A1, BTN3A1) using the tagSNP selection approach based on linkage disequilibrium (LD), as previously described [7,8,9,28,29,33].
The 768 SNPs were genotyped using a custom designed 768plex Illumina GoldenGate TM assay (Illumina Inc., San Diego, CA) following the manufacturer's instructions. The BeadStudio 2 software was used to call genotypes. The genotyping quality was high, with a mean SNP call rate of 99.7% and 99.8% (for the Rochester cohort and the San Diego cohort, respectively), and a mean subject call rate of 99.7% and 99.8%. The estimated concordance of duplicates was 96%-100%.
For the Rochester cohort, after removing samples that failed genotyping and/or QC checks, there were 1,039 samples and 606 SNPs available for analysis, including 887 samples and 555 SNPs for the Caucasian-only analysis. For the San Diego cohort, we had data (available for analysis) on 989 subjects and 611 SNPs, including 542 samples and 565 SNPs for the Caucasian-only analysis.

Statistical Methods/Analysis
Summaries of the subject characteristics were obtained within the two study cohorts, counts and percentages for categorical features and medians and interquartile ranges (25th and 75th percentiles) for continuous features. Principal components analysis of a large collection of independent SNPs was applied to estimate and define racial groups, as previously described [34,35]. As the Rochester cohort was predominantly Caucasian, the data for both cohorts was subset to this race. Hardy-Weinberg equilibrium (HWE) was tested within the Caucasian subsets of each cohort. Principal component analysis was repeated via EIGENSTRAT [36] on a collection of SNPs except those that were in low linkage disequilibrium (r 2 ,0.1) and that had HWE p-values.1610 23 .
Non-autosomal SNPs, or SNPs that had a minor allele frequency less than 0.01, were not included in these population stratification analyses, which produced additional race and cohort-specific population stratification eigenvectors.
Linear regression models, which modeled the quantitative neutralizing antibody outcome as a linear combination of covariates and SNP genotype variables, were used to assess associations between neutralizing antibody levels and SNP genotypes. An additive genetic model, which coded SNP genotypes into a single variable reflecting the number of minor alleles carried by each individual, was used for the primary tests of significance. As sensitivity analyses, dominant and recessive genetic effects were also examined for their associations with the outcome, but for most of the SNPs these alternative models did not prove to be superior compared to the additive genetic model used in our primary assessment of significance. Neutralizing antibody (NT 50 ) levels were analyzed as the phenotype of interest after logarithmic transformation in order to conform to linear models assumptions. SNP associations were assessed while adjusting for gender, assay batch effect, quartiles of age at enrollment, immunization age, time since last immunization to enrollment, and race-specific population stratification eigenvectors.
Associations were tested independently within each cohort, and the results were compared between the two cohorts. In order to determine which SNPs were potentially associated with neutralizing antibody titers across the two study cohorts, a fixed effects meta-analysis was performed for each SNP [37]. These metaanalyses used estimates of effect from the two study cohorts, and combined them using a fixed effects paradigm, to obtain a pooled estimate of effect, along with its pooled standard error and test of significance. Additionally, the meta-analyses provided a test of heterogeneity between cohorts. SNPs that did not display significant heterogeneity between the two cohorts (p.0.1) and that had meta-analysis p-values less than 0.05 were considered to be potentially associated with neutralizing antibody levels, and meriting further follow-up in future studies.
Analyses were carried out using the SAS software system (SAS Institute, Inc., Cary, NC) and the PLINK whole genome data analysis software package.

Demographic Characteristics and Immune Variables of the Study Subjects
The demographic characteristics of the study cohorts have been reported previously [7,9,25,26,27,28,29,30,38], and are summarized in Table 1. Briefly, out of the Rochester (discovery) cohort, 1,029 subjects met all inclusion, exclusion and QC criteria, and had genotyping and rubella neutralizing antibody data available for this study. Of these 1,029 subjects, 887 (85.4%) were Caucasians, and were included in the final analysis. Out of the San Diego (replication) cohort, 985 subjects met all inclusion, exclusion and QC criteria, and had genotyping and rubella neutralizing antibody data available. Of these 985 subjects, 542 (55.03%) were Caucasians, and were included in the final analysis. The median age at enrollment was 15 years (IQR, 13, 17) for the Rochester cohort, and 23 years (IQR, 22,27) for the San Diego cohort. The median neutralizing antibody level was 55.5 NT 50 for the Rochester cohort and 62.8 NT 50 for the San Diego cohort ( Table 1) Detailed characterization of the immune response variables in the study cohorts and their interrelationships has been published previously [25,26,27].

Associations between SNPs and Antibody Titers in Two Distinct Cohorts
We found three replicated genetic associations between SNPs tagging the poliovirus receptor gene (PVR, nectin-like protein 5, CD155) in high linkage disequilibrium/LD (rs1550551, rs2165538 and rs73936843; r 2 = 1, Fig. 1), and rubella-specific neutralizing antibodies (meta-analysis p-value = 0.008, Table 2). Although relatively rare, the heterozygous genotype of these PVR SNPs (no homozygous minor allele genotype was observed) was consistently associated with a significant decrease (. 2-fold) in neutralizing antibody titers after rubella vaccination in both cohorts. Similarly, the heterozygous genotype of two SNPs (rs78245864 and rs41290122, in LD; r 2 = 0.66, Fig. 1) in the PVRL2 gene (poliovirus receptor-related 2, herpesvirus entry mediator B, nectin 2, CD112) was associated with a significant increase (, 30%) in rubella-specific neutralizing antibody titers (meta-analysis p-values of 0.038 and 0.04 for rs78245864 and rs41290122, respectively, Table 2). Furthermore, five additional PVRL2 genetic variants (some in LD, Fig. 1) demonstrated some evidence for association with rubella-specific NT 50 titers post vaccination (meta-analysis p-value,0.044, Table 2), although some of the findings (association with high/low response) were not consistent between the two cohorts. Our meta-analysis also revealed potential associations between SNPs in the promoter regions of CD209/DC-SIGN (rs2287886) and IL6 (rs1880241) genes and rubella-specific neutralizing antibodies (meta-analysis pvalues of 0.008 and 0.027, for rs2287886 and rs1880241, respectively, Table 2). We found a replicated genetic association for a SNP (rs1153600) in intron 3 of the vitamin A receptor gene RARB that demonstrated a significant decrease in neutralizing antibody titers with the representation of the minor allele (metaanalysis p-value = 0.037, Table 2). The homozygous minor allele genotype of a SNP (rs16895223) in the putative rubella virus receptor gene MOG (myelin oligodendrocyte glycoprotein) was associated with a decrease in rubella-specific antibodies (metaanalysis p-value = 0.046, Table 2). In addition, BTN2A1 SNP rs1977198, IRF9 SNP rs17256713 and EIF2AK2 SNPrs4648212 also demonstrated evidence for association with rubella-specific NT 50 titers (p,0.049, Table 2).
The analysis results for all SNPs assessed for association with neutralizing antibody levels in our current study are presented in Table S1.

Discussion
Much knowledge has been accumulated on the role of genetic factors in the inter-individual immune response variation following immunization with measles-mumps-rubella (MMR) vaccine. A genome-wide transcriptional study (using mRNA-Seq) in high and low antibody responders to rubella vaccine identified HLA genes, and several innate immunity and inflammation-related genes (MEFV, Mediterranean fever gene, EMR3, EGF-like module containing, mucin-like, hormone receptor-like 3 gene, etc.) that discriminated between high and low humoral immune response following rubella immunization [4]. We and others have highlighted the importance of HLA polymorphisms, SNPs in cytokine and cytokine receptor genes, antiviral genes, toll-like receptors and pathway signaling genes, vitamin A and D receptor and cellular viral receptor genes (CD46 and SLAM for measles virus), and other important innate immune genes (TRIM genes, DC-SIGN, etc.) for inter-individual immune response variation after measles and rubella vaccination [2,5,6,7,8,9,10,21,22,23,24,28,29,39,40,41,42]. Because the immune response encompasses a system of interrelated components and biological processes, it is likely that multiple genetic factors and gene networks influence host response and vaccineinduced immunity (involving antibody response), including measles and rubella vaccines [5,40].
The purpose of the current study was to refine our previous findings through replication, and to examine, in two independent cohorts, the role of newly discovered host viral receptors, attachment factors and immune regulators in the adaptive immune response induction and/or maintenance following vaccination with a live rubella vaccine.
Our meta-analysis results provide evidence for the involvement of multiple genetic variants (some in LD) within the PVR-PVRL2 gene region (Fig. 1) on chromosome 19 (19q13.31-32) in the genetic control of neutralizing antibody response to rubella vaccine. The poliovirus receptor (PVR, CD155) and poliovirus receptor-related 2 (PVRL2, CD112) genes encode transmembrane glycoproteins/nectins of the immunoglobulin (Ig) superfamily, which are components of adherens junctions and serve as cellular entry receptors for poliovirus (PVR), pseudorabies virus, and certain mutant strains of herpes simplex virus (PRRL2) [17,18,19]. As virus-specific cell surface receptors, these proteins have important functional downstream effects on processes related to viral entry, such as viral replication, cell-to-cell spread, viral tropism, pathogenesis and antiviral immunity. For example, the Ala67Thr mutation in the poliovirus receptor gene (PVR) was previously associated with a higher risk for developing vaccineinduced and wild-type virus-induced poliomyelitis [43].
In summary, based on the literature and our findings, we speculate that one or more of the multiple PVR and PVRL2 genetic variants, associated with variation in rubella vaccinespecific neutralizing antibody levels, or other tagged causal SNP/ SNPs in the 19q13 region, influence the gene expression and/or protein function of CD155 (PVR) and CD112 (PVRL2). In turn, this may have downstream effects on rubella virus-specific humoral immunity by either modulating the unspecific (or specific) virus attachment/entry, replication and immune response priming or, more likely, by modulating Th1/Th2 cell differentiation and promoting more (or less) efficient antibody response [16].
The only cellular receptor for rubella virus identified to date is the myelin oligodendrocyte glycoprotein (MOG), an adhesion molecule that is important in nerve myelination and is implicated as a target antigen in the pathogenesis of several autoimmune demyelinating diseases (including multiple sclerosis/MS), and in the central nervous system (CNS) damage in congenital rubella syndrome [11]. However, rubella virus cell tropism and infectivity suggest that other, yet undiscovered, cellular receptors and adhesion factors capable of mediating infection exist for this pathogen. MOG is expressed predominantly in the central nervous system and is hardly detectable in immune system organs/tissues and cells such as spleen, thymus and PBMCs [11]. Our analysis identified two intronic polymorphisms (rs16895223, p = 0.046 and rs1977198, p = 0.027, Table 2) in the MOG gene and the highly homologous butyrophilin BTN2A1 gene, both on chromosome 6p22, that were associated with variations in rubellaspecific neutralizing antibody response. Host genetic variation in these or other tagged functional MOG/BTN2A1 SNPs may potentially impair or enhance viral entry into susceptible cells, alter virus replication/propagation and antigen abundance and modulate subsequent immune response. Alternatively, these genes may interfere in immune regulation, as demonstrated for the butyrophilins. The latter are co-stimulatory molecules that are implicated in T cell inhibition (negative regulation of T cell proliferation, activation, cytokine production), immune signaling and interaction with DC-SIGN on dendritic cells (BTN2A1 is identified as a ligand for DC-SIGN), and genetic variation in these genes have been associated with inflammatory diseases [48,49]. Of note, our results also demonstrate an association between a promoter SNP in DC-SIGN (rs2287886, p = 0.008, Table 2) and humoral immune response variations following rubella vaccination. DC-SIGN (Cell-Specific Intercellular adhesion molecule-3-Grabbing Non-integrin, CD209) is a C-type lectin that mediates dendritic cell function and activation of CD4 + T cells and binds multiple microorganisms, including viruses, by recognizing mannose type glycoproteins [20]. This molecule is exploited by many viruses for attachment and/or cellular entry: measles virus (MV); rift valley fever virus (RVFV); influenza A viruses; human immunodeficiency virus (HIV-1); hepatitis C virus (HCV); herpes simplex virus (HSV1); human cytomegalovirus (HCMV); Ebola virus; SARS coronavirus and dengue virus [20]. Polymorphisms in DC-SIGN (and, in particular, the promoter SNP rs2287886) have been associated with the clinical course and outcome of dengue virus infection, cytomegalovirus infection, tick-borne encephalitis virus (TBEV) infection and pulmonary aspergillosis [50,51,52]. Furthermore, the promoter SNP rs2287886 is functional and was demonstrated to influence DC-SIGN gene expression and HCMC infection in dendritic cells [52]. Interestingly, we found the same DC-SIGN promoter SNP to be associated with variations in measles vaccine-induced TNFa secretion in African-Americans and several other DC-SIGN SNPs that were associated with variations in measles-specific neutralizing antibody levels in Caucasians and African-Americans [21]. Therefore, it is likely that DC-SIGN, PVR, PVRL2 or other attachment factors and/or Figure 1. Locus zoom plot and schematic LD block structure of the PVR and PVRL2 genetic variants of interest. A. Haplotype block structure of the PVR and PVRL2 genetic variants, analyzed using Haploview software, version 4.2 (all SNPs presented were directly genotyped in the study). The r 2 color scheme is: white (r 2 = 0), shades of grey (0,r 2 ,1), black (r 2 = 1). The numbers report the r 2 value multiplied by 100. B. Locus zoom plot of PVR/PVRL2 genetic region on 19q13, showing multiple SNPs associated with rubella-specific neutralizing antibody levels. P-value is represented on the left Y-axis. The recombination rate is represented on the right Y-axis and as a line at the bottom of the graph. SNPs are represented as circles, color-coded for LD with respect to rs73936843 (the top three PVR SNPs rs73936843, rs1550551 and rs2165538 are in perfect LD and overlap, as represented on the graph   Immune outcomes for the Rochester cohort, presented as median neutralizing antibody titer in NT 50 6 IQR (inter-quartile range), as measured by the sICNA assay. c Ordinal p-value for analysis in the Rochester cohort only, after adjusting for gender, assay batch effect, quartiles of age at enrollment, immunization age, time since last immunization to enrollment, and race-specific population stratification eigenvectors. d Genotypes for the San Diego cohort, presented as homozygous major allele/heterozygous/homozygous minor allele. N indicates the number of subjects with the specific genotype. e Immune outcomes for the San Diego cohort, presented as median neutralizing antibody titer in NT 50 6 IQR (inter-quartile range), as measured by the sICNA assay. f Ordinal p-value for analysis in the San Diego cohort only, after adjusting for gender, assay batch effect, quartiles of age at enrollment, immunization age, time since last immunization to enrollment, and race-specific population stratification eigenvectors.
g Meta-analysis p-value after adjusting for gender, assay batch effect, quartiles of age at enrollment, immunization age, time since last immunization to enrollment, and race-specific population stratification eigenvectors; pooled estimate from the meta analysis showing the magnitude and direction of the estimated effect on the immune measure; homogeneity test p-value estimating the homogeneity between the two cohorts. Only SNPs with p-value below 0.05 in the meta-analysis and homogeneity test p-value above 0.1 are presented. doi:10.1371/journal.pone.0099997.t002 cellular viral receptors and pathway-related genes are involved in the genetic control of rubella vaccine-induced immune response heterogeneity after vaccination. Another important finding is the observed genetic association between a SNP (rs1153600) in intron 3 of the retinoic acid receptor, beta (RARB) gene and rubella-specific neutralizing antibody response (the homozygous minor allele and the heterozygous genotypes were associated with lower antibody titers in both cohorts). RARB is a receptor for retinoic acid, the biologically active metabolite of vitamin A, with involvement in cell growth, differentiation, cell signaling and gene regulation. Interestingly, vitamin A has direct effects on the immune system and regulates antigen presentation, lymphocyte homing, proliferation of B and T cells, T-helper differentiation, cytokine production, T cell activation and cytotoxicity, as well as enhancing the antibody response to vaccination, including the measles vaccine, oral polio vaccine, and tetanus and diphtheria toxoids [53,54]. Furthermore, our previous studies demonstrated associations between: the same RARB SNP (rs1153600, p = 0.034) and two other RARB SNPs and rubella (whole) virus-specific IgG levels (as measured by a chemiluminescent immunoassay) in 738 children/young adults following two doses of MMR vaccine [8]; and RARB SNPs/haplotypes and variations in measles-specific neutralizing antibody levels and cellular immune outcomes after MMR vaccination [41]. The findings from the current study and previous reports add evidence to the importance of vitamin A receptors and pathway in immune response regulation after immunization.
Other potential findings include SNP associations between polymorphisms in the IL6 gene (promoter region SNP rs1880241), IRF9 gene (39intergenic SNP rs17256713) and EIF2AK2 gene (intronic SNP rs4648212), which we previously found to be associated with whole rubella virus-specific IgG levels using single-SNP and/or multigenic assessment analysis [5,7]. Other SNP associations identified in our previous studies were not confirmed in the current study, possibly due to differences in study design, sample size and cohorts' characteristics, analytical approach, differences in immune response measurement (rubella whole virusspecific IgG levels vs. rubella virus-specific neutralizing antibody levels), and/or false-positive findings [4,5,6,7,8,9,10].
While our efforts have identified a collection of genetic variants that are of high interest as potentially playing a role in variability of neutralizing antibody response after rubella vaccination, it must be noted that some of the genetic associations may not be true positives. The limitations of our study reflect the fact that the characteristics of the study participants were different between the two study cohorts. This made it difficult to compare and/or combine the results from the two cohorts using a method that would typically be applied to control for false-positive findings. For this reason, we applied meta-analytic approaches to evaluate the degree of evidence of association present when evaluating combined results across the two cohorts, with the understanding that false-positive results would be less likely to show consistent associations between cohorts. Given the signals that we observed, and the biological plausibility, it is likely that at least several of the identified associations/genes do indeed play a role in influencing neutralizing antibody levels observed in vaccinated individuals. In addition, our analysis focused on additive genetic effects, and while the support for other genetic models over the additive model was not compellingly strong in our sensitivity analyses, it is possible that for some of the genetic variants an alternative genetic model may prove superior. Other limitations include the influence of unpredictable factors confounding analysis results (e.g., potential wild-type rubella virus exposure in deployed military personnel from the San Diego cohort).
The strengths of our study include the use of a clinically relevant immune outcome to assess vaccine response (i.e., rubella-specific functional/neutralizing antibody levels measured using a state-ofthe-art, standardized, high-throughput immune assay). A significant strength is the genetic association analysis approach, adjusting for known confounding variables, and the use of two independent cohorts (subsetted to Caucasian subjects only) to minimize the false-positive findings. The identified genes and genetic variants will be further fine mapped in order to identify candidate functional variants in the genomic regions tagged by the original SNP (if not functional). The p-values, and magnitudes of effect, obtained from the fine-mapping analyses will provide important insights into which of the variants is most likely to be causing the genetic association that was originally observed. Finally, the identification of likely causal genetic variants will be complemented with functional studies to reveal the SNP immediate functional effect and the downstream consequences/biological mechanisms for humoral immune response variation.
In summary, our results identified multiple genetic variants, primarily in genes related to viral attachment and entry, and/or immune regulation that are associated with inter-individual differences in rubella-specific neutralizing antibody response after vaccination. Such findings may indicate evidence for novel receptors used by the rubella virus for cell entry, and will assist in development of improved vaccines and vaccine formulations for achieving optimal immune response after vaccination.

Supporting Information
Table S1 SNPs assessed for association with neutralizing antibody levels after rubella vaccination in two different cohorts. (DOCX)