Gene-Gene Interaction and Functional Impact of Polymorphisms on Innate Immune Genes in Controlling Plasmodium falciparum Blood Infection Level

Genetic variations in toll-like receptors and cytokine genes of the innate immune pathways have been implicated in controlling parasite growth and the pathogenesis of Plasmodium falciparum mediated malaria. We previously published genetic association of TLR4 non-synonymous and TNF-α promoter polymorphisms with P.falciparum blood infection level and here we extend the study considerably by (i) investigating genetic dependence of parasite-load on interleukin-12B polymorphisms, (ii) reconstructing gene-gene interactions among candidate TLRs and cytokine loci, (iii) exploring genetic and functional impact of epistatic models and (iv) providing mechanistic insights into functionality of disease-associated regulatory polymorphisms. Our data revealed that carriage of AA (P = 0.0001) and AC (P = 0.01) genotypes of IL12B 3′UTR polymorphism was associated with a significant increase of mean log-parasitemia relative to rare homozygous genotype CC. Presence of IL12B+1188 polymorphism in five of six multifactor models reinforced its strong genetic impact on malaria phenotype. Elevation of genetic risk in two-component models compared to the corresponding single locus and reduction of IL12B (2.2 fold) and lymphotoxin-α (1.7 fold) expressions in patients'peripheral-blood-mononuclear-cells under TLR4Thr399Ile risk genotype background substantiated the role of Multifactor Dimensionality Reduction derived models. Marked reduction of promoter activity of TNF-α risk haplotype (C-C-G-G) compared to wild-type haplotype (T-C-G-G) with (84%) and without (78%) LPS stimulation and the loss of binding of transcription factors detected in-silico supported a causal role of TNF-1031. Significantly lower expression of IL12B+1188 AA (5 fold) and AC (9 fold) genotypes compared to CC and under-representation (P = 0.0048) of allele A in transcripts of patients' PBMCs suggested an Allele-Expression-Imbalance. Allele (A+1188C) dependent differential stability (2 fold) of IL12B-transcripts upon actinomycin-D treatment and observed structural modulation (P = 0.013) of RNA-ensemble were the plausible explanations for AEI. In conclusion, our data provides functional support to the hypothesis that de-regulated receptor-cytokine axis of innate immune pathway influences blood infection level in P. falciparum malaria.


Introduction
Infection with Plasmodium falciparum is still a major health problem worldwide, causing about 225 million new malaria cases each year [1]. Majority of the patients with parasite infection present febrile symptoms while a subgroup develops life-threatening complications such as severe malarial anemia (SMA) or cerebral malaria (CM). Malaria parasitemia which affects disease severity and transmission is controlled by a balance between depletion of red blood cells and immune clearance by T-helper cells, B lymphocytes, and cytokines as shown in murine model of Plasmodium chabaudi infection [2][3][4]. The role of innate immunity to restrict parasite growth in rodent model is further demonstrated by a recent study which shows when parasite dose saturates the capacity of innate response; experimentally enhanced innate immunity can control parasite density indirectly by depletion of RBCs [5]. In humans, hemoglobin degradation and heme detoxification by the obligate intracellular malaria parasite results in the formation of hemozoin (HZ) which along with malarial glycosylphosphatidylinositol (GPI) prime innate immune response by the production of pro-inflammatory cytokines through toll-like receptor mediated signaling [6,7]. Two important functions of innate immunity in the defense against the parasite are (i) it triggers a battery of pro-inflammatory cytokines inhibiting rapid parasite growth and thereby limiting the onset of malaria pathology and (ii) it determines the type and efficiency of subsequent parasite specific adaptive immune responses through the cytokine mediators at later stages of infection [8][9][10]. An account of the etiological components of innate immunity and the mechanistic framework of their interactions will offer the much needed therapeutic alternatives for control of the global burden of malaria.
Genetic studies correlating malaria susceptibility with host immunity have a long history. Majority of the studies on the role of host genetics in determining malaria susceptibility have utilized population based case-control design and candidate or genome wide markers [11][12][13][14][15][16][17][18][19][20]. A general concern about the association study is that the alliance between disease phenotype and genetic loci remains tentative unless a functional causality is established. The major objective of this study is to understand the interactions and functional contribution of host genetic factors associated with a measurable phenotype, namely parasitemia, to extract an improved insight of host regulatory mechanisms controlling the blood infection level in P. falciparum mediated malaria. We have published a genetic association study of this phenotype analyzing fourteen single nucleotide polymorphisms (SNPs) located on genes encoding toll-like-receptor (TLR)-2, 4 and 9 and tumor necrosis factor-a (TNF-a) and lymphotoxin-a (LTA) and reported significant correlation of blood parasitemia with TLR4 non-synonymous and TNF-a promoter polymorphisms in Indian patients with mild malaria [21]. Here we extend this analysis by investigating genetic association between polymorphic variability of IL12B, which encodes the IL-12p40 subunit of IL12, with blood parasitemia. IL12 is a heterodimeric pro-inflammatory cytokine with pleiotropic effects, acting as a potent immune-regulatory molecule and hematopoietic growth factor in infections caused by Plasmodium parasites [22]. It is produced by phagocytes and dendritic cells in response to pathogens through toll-like and other extracellular receptors. Physiologically the most important target cells of IL12 are hematopoietic progenitors, NK cells and T cells for which it induces proliferation and production of type-1 cytokines (e.g. IFNc). IL12 and IFN-c together enhance activation and production of T h 1 associated classes of immunoglobulin from B cells [23]. IL12 is composed of a 35 kD subunit encoded by IL12A and 40 kD subunit encoded by IL12B. Polymorphisms in genes encoding both the subunits have been reported in a wide range of immune and inflammatory diseases including malaria [24][25][26][27][28][29][30][31][32][33][34].
Given the complexity of parasite biology and host immune system, it is unlikely that genetic variation of a single locus would provide an adequate explanation of inter-individual differences of host immune response which results in diverse clinical manifestations. To this end, identification of gene-gene interactions could enhance the power and accuracy of predicting disease outcome of a complex disorder. We have collated the genetic data on IL12B from the present report with that on toll-like receptors and cytokine loci from our previous study [21] to model the possible genetic interactions that may account for the differences in blood parasitemia in malaria patients. The functional relevance of receptor-cytokine epistatic models has been captured by analysis of receptor-genotype dependence of cytokine expression in vivo.
For a better description of genetic architecture of disease susceptibility and unambiguous identification of factors responsible for both causality and predisposition to a disease, functional appraisal of disease-associated polymorphisms is essential. There is widespread recognition that differences in gene expression may be an important source of phenotypic diversity in complex diseases [35][36][37] and that non-coding polymorphisms contribute to the variance and etiology of a trait by regulating the expression of nearby genes [38,39]. To explore the plausible regulatory mechanisms exerted by cytokine SNPs we have characterized the allele specific events by studying their transcriptional differences in terms of reporter gene activities and allelicexpression-imbalance (AEI). Our study provides detailed insights into molecular effects of cis-regulatory variants in controlling cytokine gene expression in P. falciparum mediated malaria. However it underscores the possibility that this complex trait involves even more complex regulatory intricacies than previously anticipated.

Patient recruitment and Laboratory measures
A total of 293 mild malaria patients (age = 16-37 years) with Plasmodium falciparum infection were recruited from The Calcutta School of Tropical Medicine between September 2008 and January 2009 following WHO guidelines after obtaining the written informed consent from each study participant [40]. Patients with bacteremia, measles, acute lower respiratory tract infection, severe diarrhea with dehydration and other chronic or severe diseases such as cardiac, renal or hepatic diseases, HIV/ AIDS were excluded from the study [21,40]. Appropriate approvals have been obtained from Institutional ethical committees of University of Calcutta and The Calcutta School of Tropical Medicine, India. Parasitemia status of each patient was determined during their first visit to the clinic using Giemsa-stained blood smears and oil immersion microscopy. Detailed description of study samples and procedure for parasite enumeration can be found in our earlier report [21]. For genetic epidemiology, the patient pool from our previous study was used while for functional analyses case samples were enrolled (N = 64, age = 4-16 years) and registered in the year 2010 during July to September, from Calcutta National Medical College & Hospital, Kolkata, India as per WHO 2006 guidelines [40]. The blood samples for genetic and functional analyses were collected before any medical interventions.

Genotyping IL12B polymorphisms
Genomic DNA was extracted from peripheral blood leukocytes using a QIAamp DNA Blood Kit (Qiagen, Hilden, Germany). IL12B promoter (rs17860508) and 39UTR (rs3212227) polymorphisms were assayed through PCR followed by an allele specific restriction enzyme digestion. The primer pairs for genotyping were listed in Table S1. PCR amplification was carried out in 10 ml reaction mixtures containing 1 U AmpliTaq Gold TM DNA polymerase, 1.5 mM MgCl 2 , 250 mM of each dNTP, 5 pmoles of each primer (Sigma Aldrich, St. Louis, MO). The cycling parameters consisted of an initial denaturation at 96uC for 5 mins, followed by 40 cycles of denaturation at 96uC for 30 sec, annealing at 60uC for 30 sec, extension at 72uC for 30 sec, and then completed with a final extension at 72uC for 5 mins using a thermal cycler (Applied BiosystemsH GeneAmpH PCR System 9700). PCR amplicons were digested with AluI (rs17860508) and TaqI (rs3212227) (Fermentas International Inc.) at 37uC and 65uC respectively according to manufacturer's protocol. Digestion patterns were analyzed by electrophoresis in agarose gel. The accuracy of the PCR-RFLP assay was confirmed for each locus by randomly selecting 10% of total samples and sequencing the PCR products on ABI Prism 3100 Genetic Analyzer using Big-Dye Terminator v3.1 (Applied Biosystems, Foster City, CA).
Statistical analysis. All the statistical tests were performed using SPSS version 10.0 and R-program (version 2.0 package). Associations between alleles/genotypes and parasite estimate were considered ''strong'' and ''statistically significant'' in cases where the P-value,0.05.
(a) Single Gene Analysis: Distribution of allele and genotype frequencies and genetic risk assessment of IL12B polymorphisms. The allele and genotype frequencies for IL12B SNPs were estimated by gene counting and Hardy-Weinberg equilibrium was evaluated by x 2 test using Haploview (http://www.broad.mit.edu/mpg/haploview/). Pairwise comparisons between log-transformed parasitemia and genotypic classes for IL12Bpro and +1188 loci were analyzed by model based independent t test and ANOVA F-ratio. The results were further validated using non-parametric Wilcoxon Rank sum and Kruskal Wallis H tests. Bonferroni correction was done for multiple testing. For allele based risk assessment, the unprocessed parasitemia data was previously dichotomized into high (N = 88) and low (N = 205) parasitemic groups employing the two component mixture model and expectation-maximization algorithm and a threshold value of 7520 parasites/ml or 3.88 log-parasitemia was derived as the cutoff. Detailed procedure of this strategy was described earlier [21]. The proportion of the risk allele for each locus was compared between high and low parasitemic groups using a two-way contingency table.
(b) Multi-loci Analyses: Multifactor Dimensionality Reduction analysis. We adopted the multifactor dimensionality reduction (MDR) method upon merging the IL12B genotype data with those from our previous genetic epidemiologic study for identification of non-additive epistatic interactions [21,[41][42][43][44][45][46]. MDR is a non-parametric approach which converts multiple variables into a single attribute, thereby changing the representation space of the data. To run MDR, eight out of sixteen loci from the assembled dataset were selected as they showed an increasing/decreasing trend of log-parasitemia distribution across three genotypes which allowed us to code the genotypes into distinct risk groups. The ''low'', ''moderate'' and ''high'' risk genotypic classes for each locus were designated as ''0'', ''1'' and ''2'' while the ''low'' and ''high'' parasitemic groups were assigned as ''0'' and ''1'' respectively to depict the disease status. Since an equal number of sample size in high and low risk groups was a prerequisite for execution of MDR, the sample size of our analysis was kept to eighty eight, the number that represented in the high parasitemic group. Of the 205 individuals belonging to low parasitemic group, 100 distinct files each comprising of 88 individuals were generated using random numbers as seeds using R version 2.0 package to minimize the chance of detecting spurious results due to any bias. This was run on MDR software to obtain 100 independent outputs (version 2.0 beta 8.2). From eight genetic attributes (SNPs) namely TLR4Thr399Ile (rs4986791), IL12B+1188 (rs3212227), TNF-1031 (rs1799964), TNF-857 (rs1799724), TNF489 (rs1800610), LTA80 (rs2239704), LTA252 (rs909253), and TLR9P545P (rs352140), the possible multilocus genotypes were represented in a multidimensional contingency tables and classified as ''high risk'' if the ratio between two groups exceeded a threshold tested over the entire training data set. After reducing the dimensionality, all factor combinations were evaluated for their ability to classify the disease status in the training dataset and the best combination of factors with the minimum prediction error (1-testing balance accuracy) was calculated from the test dataset. Then the entire data was partitioned into 100 different subsets for cross-validation, from which ninety nine out of 100 subsets were assigned as a training balanced set (99/100) while the remaining one hundredth dataset (1/100) was termed as testing balanced set. To assign statistically meaningful gene-gene combinations, 100 MDR files were permutated 10,000 times using MDR permutation testing software (version 1.0 beta 2.0). The models with P value,0.05 and .95% average cross-validation consistency (CVC) were regarded as the best predictive models. To infer the genetic risk of the multifactorial models, we used odds ratio (OR) as a measure of association between genetic polymorphisms and blood infection intensity. The value of ORs was logarithmically transformed and standard errors were derived from the confidence intervals (CI). Genotype counts of the constituent loci between high and low parasitemic groups were used to calculate ORs and their corresponding 95% CIs for each gene-gene interaction model. A forest plot (blobbogram) was used to present the ORs and their 95% CIs estimated from the two-way interaction models where the odds ratio was denoted by a dot and the width of the horizontal line represented the 95% CI for the estimated OR.

RNA extraction from PBMC fraction and cDNA synthesis
Peripheral blood mononuclear cells (PBMCs) were isolated from fresh whole blood collected from malaria patients using Histopaque 1077 (Sigma Aldrich, St. Louis, MO) double-gradient density centrifugation. Total RNA was extracted with TRIH Reagent (Sigma Aldrich, St. Louis, MO) and dissolved in DEPC treated water (Bioline). RNA (1 mg) samples were treated with 2 U of RNAse free DNAse I (Fermentas Life Sciences, UK) and incubated at 37uC for 30 mins to remove DNA contamination, reverse transcribed using random hexamers and High Capacity cDNA Reverse Transcription kit (Applied Biosystems Inc.). The optimized condition for cDNA preparation was 10 mins at 25uC, 120 mins at 37uC followed by heating at 85uC for 5 mins in a thermal cycler (Applied BiosystemsH GeneAmpH PCR System 9700) and stored at 220uC.

Quantitative Real Time PCR analyses
To scrutinize the genetic impact of TLR4 non-synonymous polymorphism (rs4986791) on pro-inflammatory cytokine (IL12B and LTA) production, we compared the expression of IL12B and LTA genes under different TLR4Thr399Ile genotypic background. Since TLR4Ile-Ile and TLR4 Thr-Ile genotypes together and independently showed significantly low parasitemia compared to that of TLR4Thr-Thr, we pooled individuals having genotypes Ile-Ile with Thr-Ile genotypes and examined the pattern of cytokine gene expression with respect to patients harboring TLR4Thr-Thr genotype. The relative expression of IL12B was estimated for three +1188 genotypes by qRT-PCR using 18S rRNA as an endogenous control. The sequences of oligonucleotides used in the expression analyses were listed in Table S1. A 1:10 fold dilution of cDNA samples were used as the template and all qPCR reactions were carried out in a 10 ml reaction volume with 5 ml of SYBR Green Master Mix (Applied Biosystems) with optimized concentrations of specific primers using Applied Biosystems 7900HT Fast Real-Time PCR System. The thermal cycler was programmed for an initial denaturation step of 5 min at 95uC and followed by 40 thermal cycles of 30 sec at 95uC, 30 sec at 60uC and 30 sec at 72uC. The experiments were carried out in triplicate including the non template controls each time. Specificity of PCR amplification for each primer pair was confirmed by melting curve analysis [47] and the relative quantification (RQ) data was analyzed between groups using the DDCt method. Statistical analysis on relative expression levels was estimated using Mann-Whitney U-test (with exact probabilities) for independent samples using SPSS (ver. 10.0., SPSS Inc. Chicago, IL).

Allelic-Expression-Imbalance analysis
To examine the differential expression of IL12B 39UTR A and C alleles, the peak heights of this transcribed SNP were compared in patients heterozygous (AC) for IL12B+1188 polymorphism by resequencing both genomic DNA (gDNA) and complementary DNA (cDNA) counterpart. IL12Bexp primers (Table S1) were used to amplify and resequence the 167 bp IL12B region encompassing the polymorphic site. Sequencing was carried out on ABI Prism 3100 Genetic Analyzer using BigDye Terminator (BDT) v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA) with 0.8 ml of the PCR product, and 5 pmoles of the primer in a 10 ml reaction mixture. The raw sequence files (.ABI) were analyzed by PeakPicker Software, originally developed and kindly provided by Dr. T. Pastinen's group [48]. Genomic DNA sequence from each sample was selected as reference and the corresponding cDNA was aligned with a default cutoff of 70%. The SNP was selected manually and the peak heights of the alleles were estimated using PeakPicker. The significant difference of IL12B+1188 A/C allele ratio between the gDNA and cDNA was tested through a bivariate Sign test in SPSS version 10.0.

Construction of reporter fusion plasmids
(a) For TNF-a promoter assay. Five TNF-a promoter haplotypes pertaining to four SNPs (TNF-1031, TNF-857, TNF-308 and TNF-238) were selected for functional evaluation. The 1.2 kb fragment was amplified using TNFa_MluI_F.P and TNFa_XhoI_R.P (Table S1) where MluI and XhoI recognition sequences were appended into the primers for efficient cloning into pGL3 vector. TNF-a haplotypes-T-C-G-G (Hap1), T-T-G-G (Hap2) and C-C-G-G (Hap3) were directly cloned into TA vector pTZ57R/T (InsTAclone TM PCR Cloning Kit, Fermentas) using the above primers by amplifying gDNA from patients harboring the homozygous genotypes for all four loci. Haplotypes: T-C-A-G (Hap4) and T-C-G-A (Hap5) were cloned by amplifying gDNA from patients harboring genotypes T/T-C/C-G/A-G/G and T/ T-C/C-G/G-G/A respectively, followed by screening the clones with restriction endonucleases NcoI for TNF-308 locus and BamHI (NEB Inc., Bethesda) for TNF-238 respectively. DNA amplification protocol was conditioned with an initial denaturation at 95uC for 5 min, followed by 35 cycles of denaturation at 95uC for 1 min, annealing at 60uC for 1 min, extension at 72uC for 1 min and a final extension of 15 min at 72uC in a thermal cycler (Applied BiosystemsH GeneAmpH PCR System 9700). The unique haplotypes from the TA constructs were subsequently double-digested, purified using QIAquick Gel Extraction Kit (Qiagen, Hilden, Germany) and subcloned in the upstream of pGL3-Basic vector (kindly gifted by Dr. Susanta Roychowdhury). The integrity of each polymorphic locus and directionality of all the resulting constructs were confirmed by sequencing using both vector and insert specific internal primers. Next, the reporter gene expressions were measured in three different cell lines under the presence and absence of endotoxin lipopolysaccharide (LPS) stimulation using Luciferase Reporter Assay System (Promega, Madison, WI).
(b) For IL12B39UTR assay. IL12B 39UTR sequence (1047 bp) encompassing the A+1188C polymorphic site was amplified using gene specific primers (Table S1). XhoI and NotI (Fermentas International Inc.) restriction endonuclease recognition sites were incorporated in the oligos for directional cloning. Allele specific IL12B 39UTR containing the target sites for the microRNAs was cloned into pTZ57R/T vector (InsTAclone TM PCR Cloning Kit, Fermentas) followed by subcloning into pSiCHECK2 dual luciferase reporter plasmid (kindly gifted by Dr. Soma Banerjee) in frame to the 39 end of the cDNA encoding Renilla reniformis luciferase to produce pSiCHECK2-IL12B+1188A and pSiCHECK2-IL12+1188C constructs. Approximately 100 bp upstream and downstream sequences flanking the 70 nucleotide pre-miR sequences (www.genecards.org/) were amplified for hsa-miR-545, hsa-miR-1284, hsa-miR-23a and hsa-miR-23b using appropriate primers (Table S1) and cloned into pRNAU6.1 vector (kindly gifted by Dr. Soma Banerjee). The amplicons were inserted into pRNAU6.1 vector by double digestion with BamHI and HindIII (New England Biolabs, UK Ltd.). The orientation of the recombinant vectors was checked by sequencing. Allele specific transcript stability and miRNA-mediated IL12B 39UTR expression were examined with SYBR Green based qRT-PCR analysis and by Dual-Luciferase Reporter (DLR) Assay System (Promega, Madison, WI) respectively.
(a) TNF-a promoter haplotype assay. HepG2 cells (1.0610 5 ) were plated in 12-well plates with DMEM/10% FCS medium and were transiently transfected with 0.2 mg/ml of each of the TNF-a promoter constructs (Hap1-Hap5) or empty pGL3-Basic vector cotransfected with 0.5 mg/ml pSV-b-galactosidase vector (Promega, Madison,WI) using Lipofectamine 2000 Reagent (Invitrogen, Life Technologies, Carlsbad, USA) according to the manufacturer's protocol. pGL3-Basic vector (Promega, Madison, USA) was used as control for luciferase assay and b-galactosidase vector (Promega, Madison, USA) was used as transfection control.
To determine the optimum concentration of LPS (Sigma Aldrich, St. Louis, MO) for stimulation assay, a dose response experiment was performed with cells transiently transfected with TNF-Hap1 (T-C-G-G) followed by a stimulation of LPS solutions with five different concentrations such as 250 ng/ml, 500 ng/ml, 1 mg/ml, 1.5 mg/ml to 2 mg/ml and luciferase activities were measured. 0.2 mg of each of the five reporter plasmids (Hap1-Hap5) were transiently transfected. Thirty six hours post transfection, cells were stimulated with optimized dose of LPS for 4-6 hr and cellular extracts were prepared according to the manufacturer's instructions. Luminescence was measured as relative light units (RLU) in GloMax 20/20 Luminometer (Promega, Madison, USA) using 15 ml of cell supernatant. To normalize the luciferase activity, total protein concentration in each lysate was measured by the standard Bradford (BioRad) method. Transfection normalization was performed by b-galactosidase assay colorimetrically at 420 nm using O-nitrophenyl-b-D-galactopyranoside (ONPG) as substrate (Promega, Madison, WI) in JASCO V-630 spectrophotometer. Measurements of mean 6 S.D were taken in triplicate. Fold change represents the difference in promoter activity for variant TNF-a haplotypes in comparison to the wild-type promoter (Hap1). All the assays were done in triplicates and repeated at least for three times. Similar experiments were also performed in HEK293 and HCT116 cell lines both under LPS stimulated and unstimulated conditions.
(b) IL12B mRNA stability assay by actinomycin- [1,4,7,10,13] oxatetraazacyclohexadecin-10-yl]-4,6-dimethyl-3-oxo-3H-phenoxazine-1,9-dicarboxamide), a polypeptide containing antibiotic, can inhibit transcription by binding tightly and specifically to double-helical DNA [49,50]. It has been extensively used as a highly specific inhibitor for the formation of new de novo mRNA synthesis in both prokaryotic and eukaryotic cells [51,52]. In this experiment, HepG2 cells were transiently transfected with plasmid DNA containing the IL12B 39UTR fused to pSiCHECK2 vector, as described above. To determine the optimal concentration of actinomycin-D (SRL, Mumbai, India) in our experimental system, a range of concentrations of actinomycin-D (0.1-10 mg/ml) were studied. Twenty-four hours post-transfection with wild-type +1188AUTR construct, actinomycin-D was added into the media and the cells were harvested at different time points, i.e., 0, 1, 2, 4, 6, 18 and 24 hrs. Total RNA was extracted from these cells using TRIzolH Reagent (Invitrogen Life Technologies) and IL12B mRNA level was determined by RT-PCR and normalized with GAPDH and 18S rRNA genes. RNA isolated when actinomycin-D was just added was denoted as the 0-hr time point and served as negative control. It was compared with samples obtained at different time points following actinomycin-D addition to select the adequate incubation time where maximum change in expression could be found for both alleles. To compare relative mRNA stability of IL12B 39UTR A and C containing constructs, HepG2 cells were transfected and de novo transcription was inhibited with the optimal concentration of actinomycin-D. Cells were collected after 24 hrs of actinomycin-D treatment and the mRNA levels of the transfected gene (IL12B-39UTR/Luc) were determined by quantitative real-time PCR.
(c) IL12B 39UTR-microRNA interaction assay. HepG2 cells (1610 5 ) were transiently transfected with 0.1 mg/ml of the IL12B+1188 wild type (A) and variant allele (C) containing pSiCHECK2 dual reporter constructs alone and cotransfected with either 0.2 mg/ml of empty miR-vector (pRNAU6.1) and/or with the pre-miRNA (hsa-miR-545, hsa-miR-1284, hsa-miR-23a and hsa-miR-23b) cloned in pRNAU6.1. Forty eight hours later, cells were lysed; Firefly and Renilla luciferase activities were evaluated using dual luciferase reporter assay system (Promega, Madison, WI) in a GloMax 20/20 Luminometer. Renilla luciferase activity was normalized with respect to Firefly luciferase activity and total protein produced was estimated by Bradford method. All the transfection assays were done in triplicate and repeated thrice. Measurements of mean 6 S.D were taken in triplicate. The pairwise comparisons between each mRNA-miRNA interactions were tested using Student-t test in GraphPad Prism (http://www.graphpad.com).
Bioinformatic analysis: SNP based DNA and RNA structural rearrangements (a) TNF-1031 polymorphism associated DNA two dimensional configurations. DNA folding was performed with MFOLD version 3.0 (www.bioinfo.math.rpi.edu/,mfold/ dna/) using the DNA energy rules in structure-melting simulations, with default parameters [53,54], i.e., having 1 M Na+ concentration (and 0.0 M Mg 2+ ) at 37uC. Outputs were generated in the form of structure plots based on minimum free energy (MFE) calculations, single strand (ss) frequency count and energy dot plots [55]. In addition, to identify the transcription factor binding sites (TFBS) in the genomic sequence; TFBind (http:// tfbind.hgc.jp/) [56] was used to scan the 101 bp long sequence flanking the promoter SNP (rs1799964) located at 51 st position. This positional pattern detection tool is able to attain high sensitivity and specificity of detection by capturing the dependencies between nucleotide positions within the TFBS.
(b) IL12B+1188 derived RNA structural ensemble prediction. In an attempt to improve our understanding of the structure and function relationship of the IL12B 39UTR polymorphism on RNA turnover and stability, we performed an in silico structural analysis including the wild type and variant form of the locus on IL12B mRNA. An mRNA is unlikely to adopt a single and stable conformation, but it exists in a population of structures [57]; therefore instead of evaluating the minimum free energy (MFE) values of biological sequences, those of randomized sequences are now taken into account. As RNA secondary structure might affect the functionality of component regulatory elements, RNAfold program (http://rna.tbi.univie.ac.at/cgi-bin/ RNAfold.cgi) was employed to predict secondary structures under default parameters and the folding temperature was fixed at 37uC [55,58,59]. In addition to the minimum free energy approach using RNA MFOLD and RNAfold, we adopted a partition function calculation method which evaluates the entire ensemble of possible RNA conformations for a given sequence and computes the effect of single nucleotide polymorphisms on RNA structural ensemble using SNPfold (http://ribosnitch.bio.unc.edu/ snpfold/SNPfold.html). The structures were also generated with Sfold (http://sfold.wadsworth.org/) algorithm that computes equilibrium partition functions for all substrings of an RNA sequence based on the Turner thermodynamic parameters [60,61]. Base-pair distances were used [62,63] as the basis for evaluating dissimilarity between structures where the ''centroid'' in the entire structure ensemble space has the shortest total base-pair distance [62].

Results
The patients participating in this genetic association study were recruited according to WHO 2006 guidelines [40]. All of them displayed mild febrile symptoms, the summary of demographic, clinical and laboratory characteristics of the study participants were presented earlier [21]. In short, the crude parasitemia data derived from thick and thin blood smears deviated from normal behavior (skewness = 3.671), it was processed with log-transformation (skewness = 0.433) and regressed with respect to the covariate, age (P value of correlation between age and parasite load was 0.074). For allele based genetic analyses, cases were stratified into high and low parasitemic groups on the basis of a statistically derived threshold value of 7520 or 3.88 in terms of crude and logparasite counts per ml respectively. The above threshold value was obtained employing a two component mixture model and Expectation Maximization algorithm [21]. No significant differences were observed in the distribution of age (P = 0.2632), gender (P = 0.3728), Hemoglobin count (P = 0.2971) and in average body weight (P = 0.3901) between the high and low parasitemic groups. Only peripheral parasitemia was found to be significantly different (high = 15227611014; low = 271561705, P,0.0001) between the groups.

Genetic analysis with IL12B polymorphisms
We studied the genetic variants located in the promoter and 39UTR region of the gene encoding IL12B by PCR based RFLP method. Two biallelic markers; rs17860508 (designated as IL12Bpro): a complex promoter insertion/deletion polymorphism -/G/GC/TTAGA/TTAGAG situated at 22703 bp upstream of the transcription initiation site and rs3212227 (designated as IL12B+1188): an A to C substitution in the 39UTR region at +1188 position of IL12B gene were selected to examine their probable association with the blood infection intensity in two hundred and ninety three patients exhibiting Plasmodium falciparum mild malaria from Eastern India.
Genetic association between IL12B SNPs and peripheral parasite load: Genotype & allelic IL12Bpro (MAF: 0.35760.002) and IL12B+1188 (MAF: 0.40460.002) were found to be highly polymorphic in the study population and there was no deviation from Hardy-Weinberg equilibrium (P values = 0.3934 and 0.6799 for IL12Bpro and IL12B+1188 respectively). The distribution of genotypes for IL12Bpro and IL12B+1188 loci were summarized in Table 1. For model dependent association analysis, the crude parasitemia data collected from malaria patients was log transformed and subjected to linear regression with respect to age in order to obtain the standardized residuals which were then used to examine the influence of genotypes at each locus on the parasite load under different genetic models. To find out the genetic risk conferred by a SNP, IL12Bpro and IL12B+1188 were independently evaluated as covariate using three genetic models viz., dominant ((11+12) vs 22), recessive (11 vs (12+22)), and additive (11 vs 12 vs 22), where 1 represents as the risk allele using independent t-test and ANOVA (SPSS v10.0). The results of parametric analyses were substantiated by Wilcoxon Rank Sum and Kruskal Wallis H tests for two and three way comparisons respectively using the crude parasitemia data. Box-plot diagrams ( Figure 1A & B) displayed a gradual decline in parasitemia across three genotypes of IL12B+1188 locus in which AA (11), AC (12) and CC (22) had the highest, intermediate and lowest median values respectively. Significant differences were obtained for the comparison 11 vs 12 vs 22 using ANOVA test (P value = 0.001) and for the following pairwise genotype comparisons namely (11+12) vs 22 (P value = 0.001); 12 vs 22 (P value = 0.014) and 11 vs 22 (P value = 0.0001) using independent t-test presented in Table 1. The P values in the parentheses were obtained after Bonferroni correction. Nonparametric statistical tests yielded similar trends in the results ( Table 1). The mean log parasitemia of 11 (AA = 3.6860.45) genotype and the heterozygote genotype (AC = 3.5860.41), separately and jointly, differed significantly (P value = 0.001) from that of 22 (CC = 3.4060.35), while the difference between mean parasite loads representing genotypes 11 and 12 was not statistically significant (P value = 0.153) suggesting a dominant effect of 1 (A) allele over allele 2 (C) for the IL12B+1188 locus to be the most likely explanation. Similar statistical tests were conducted on IL12Bpro (rs17860508), however none of the comparisons were significant for this polymorphism. The linkage disequilibrium estimate (r 2 ) between these two loci that were 17,249 bp apart was 0.06 (P = 0.3052) suggesting an independent segregation of the promoter and UTR polymorphisms in our population.
The allele based risk assessment for IL12B+1188 was performed by stratifying the patients into ''low'' (N = 205) and ''high'' (N = 88) parasitemic groups. A detail of this procedure was described elsewhere [21]. The proportions of two alleles A and C of IL12B+1188 compared using a two way contingency table were significantly different between the groups (P,0.001). The odds ratio (OR) of comparison (1.887, 95% CI = 1.3-2.78) was due to an excess of allele A in the ''high'' parasitemic group (Table 2). Taken together, the results of the genetic association study supported a strong influence of IL12B+1188 polymorphism on P. falciparum infection load in our samples at both genotype and allelic levels.

Evaluation of gene-gene interactions: Multilocus analyses and MDR
The genetic etiology of complex diseases is considered to involve interactions among multiple genetic variants and environmental conditions. Multifactorial Dimensionality Reduction [41] was performed to explore potential gene-gene interactions by combining the IL12B genotype data of this study with our previous one [21]. We identified five two-factors (Model I-V) and one threefactor (Model VI) gene-gene interaction models with an average cross-validation consistency (CVC) .95% and a P value,0.05 using MDR permutation testing program. The gene-gene partners of MDR models were summarized in Table 3. The interacting models that satisfied the above criteria, included TLR4Thr399Ile (rs4986391) and/or IL12B+1188 (rs3212227) in combination with either TNF-1031 (rs1799964) or LTA80 (rs2239704) or TLR9P545P (rs352140). The average prediction errors (PE) for these models ranged from 0.305 to 0.366. To estimate the genetic risk under a two-factor model, the proportions of individuals harboring risk genotypes for partner loci, were compared between Table 1. Genetic association of IL12B promoter and 39UTR polymorphisms with parasitemia at genotypic level. the low and high parasitemic groups. Since IL12B+1188 was present in five out of six genetic models the odds ratio of the second locus was estimated under IL12B+1188 risk (AA+AC) background. The comparison of odds ratio (OR) and 95% confidence intervals (CI) of two-way and the constituent interacting partners were shown in the Forest plot (Figure 2A).  (Figure 2A; ID: G, I and H). For Model II, odds ratio for LTA80 (OR = 2.113, 95% CI 1.0-4.44, P value = 0.047) under TLR4 risk background (Thr-Thr) was also increased (Figure 2A; ID: J). Notably, the difference in mean parasitemia for LTA80 between the genotype comparison (11+12) vs 22 was not significant (OR = 1.232, 95% CI 0.861-1.766, P value = 0.074) when the locus was analyzed alone (Figure 2A; ID: A). In short, except for TLR9P545P, all two-way epistatic models yielded higher risk of developing high parasitemia compared to the single marker analyses. Model VI: TLR4Thr399Ile/IL12B+1188/LTA80 did not modify the risk significantly due to reduction of sample size in the high parasite group when three loci were pooled. TLR4 encodes an extra-cellular receptor present in the antigen presenting cells (APCs) and perceives the malaria antigens to trigger a complex cascade of downstream signaling events which ultimately culminates into the production of several pro-inflammatory cytokines including IL12 and LTA [64,65]. To translate functional impact of the above MDR derived gene-gene interaction models (Model I and II) in terms of TLR4 mediated signaling efficacy, we compared the expression levels of IL12B and LTA genes under TLR4Thr399Ile risk (11; N = 16) and non-risk (12+22; N = 33) genotype backgrounds using quantitative RT-PCR analyses. Patients harboring heterozygous and the rare homozygous genotypes (Thr-Ile and Ile-Ile) showed significantly increased (2.2 fold, Mann Whitney P value = 0.001) expression of IL12B gene ( median DC t 6S.D = 18.8060.87) compared to those with homozygous (Thr-Thr) genotype ( median DC t 6S.D = 19.6360.74) as shown in Figure 2B. The difference in LTA expression between TLR4Thr399Ile risk and non-risk genotypic backgrounds did not cross the threshold of statistical significance by non-parametric Mann Whitney U test (1.7 fold, P value = 0.097), though the expression of LTA gene ( Figure 2C

Functional analysis: TNF-a promoter and IL12B 39UTR polymorphisms
In addition to provide functional support to the statistically inferred epistatic models, we attempted to delineate the molecular role of TNF-a promoter and IL12B 39UTR polymorphisms using a combination of methods including quantitative RT-PCR, allelicexpression-imbalance, mRNA stability assay, luciferase based promoter (TNF-1031) and 39UTR (IL12B+1188) activity analyses, bioinformatic prediction of putative mRNA-miRNA binding interactions and SNP based DNA, RNA secondary structure assessment.
(a) Efficiency of TNF-a promoter haplotypes under LPS stimulated and unstimulated conditions: Reporter gene assay and Bioinformatic support. To explore whether TNF-1031 polymorphism influenced the basal rate of transcription of the TNF-a gene, haplotypes pertaining to a 1.2 kb fragment spanning the proximal promoter SNPs at positions 21031, 2857, 2308 and 2238 were subjected to reporter gene based assay. Notably, haplotypes were derived using Haploview and out of the 16 possible haplotypes, only six were observed in our population with haplotype T-C-G-G (1-1-1-1) showing the highest frequency in both high and low parasitemic groups. This reference haplotype as well as four other haplotypes that differed from the former at a single nucleotide position were cloned. The frequency of the variant haplotypes were .3% in the patient samples. Figure 3A showed the map of the promoter haplotypes cloned and assayed for reporter gene activity in HepG2 cells both in presence and absence of LPS, which was known to be one of the parasiteassociated-molecular-patterns (PAMPs) recognized by TLR molecules [66]. Fold changes in variant promoter (Hap2-5) activities with respect to reference haplotype Hap1 were indicated in Figure 3B. All four variant haplotypes (Hap2-5) showed decrease in luciferase activity both in presence (250 ng/ml) of LPS (27-84%) and in absence (12-78%) of the antigen ( Figure 3B). The maximum reduction in luciferase activity [78% (2LPS) and 84% (+LPS), P,0.001] was observed for Hap3 construct (C-C-G-G) which harbored the risk allele C for TNF-1031, the locus that showed significant genotypic and allelic association with blood parasite infection [21]. Our comparative promoter assay identified Hap1 and Hap3 as the strongest and weakest promoter among the five different TNF-a haplotypes in HepG2 cell line in three independent experiments. Similar experiments were performed in HEK293 and HCT116 in triplicate to compare the promoter activities of Hap1 and Hap3 ( Figure 3C). A consistent reduction of the reporter gene expression under stimulated and unstimulated conditions (80%) for Hap3 was noted in both the cell lines. The control pGL3-Basic vector showed very low levels of relative luciferase activity in both the unstimulated and LPS-stimulated cells.   Scanning of a 101 bp region flanking 50 bp in the either direction of TNF-1031 locus by TFBind program revealed that T to C transition caused loss of binding of several transcription factors. A list of the transcription factors that putatively showed higher binding affinity for the T allele given by binding scores, corresponding consensus and signal sequences with relative nucleotide positions upstream of TNF-a coding sequence were presented in Figure 3D. The possibility of altered DNA-protein interactions in terms of changes in local secondary structure due to the sequence variation at 21031 was examined using MFOLD web server ( Figure 3E and F). Out of the seventeen possible computed foldings, thirteen (76.5%) with T at 21031 position were found to be base paired while twenty (83.3%) out of twenty four possible conformations had variant C allele locked in stems.
The minimum free energy (MFE) conformations constructed using 50 bp local DNA sequence encompassing TNF-1031 were shown for wild type and variant alleles ( Figure 3E & F). Results of these in-silico comparisons indicated a relative reduction of accessibility of a transcription factor when the risk allele was present in the sequence. Taken together, the polymorphism at 21031 of TNF-a caused a loss of promoter activity as well as transcription factor binding as per our reporter gene assay and bioinformatics analyses respectively.
(b) Impact of IL12B+1188 (rs3212227) genotypes and alleles on gene expression: a quantitative evaluation. IL12B+1188 polymorphism showed strong genetic association with the blood parasite level in the malaria patients both at genotype and allele levels. To test the hypothesis that this 39UTR polymorphism may modulate IL12B transcript turnover, relative expression of IL12B-mRNA was quantified using the cDNA synthesized from mRNA extracted from peripheral blood mononuclear cell (PBMC) fraction of whole blood collected from individuals harboring three different genotypes namely AA (N = 24), AC (N = 28) and CC (N = 12) using real-time quantitative PCR. Expression of IL12B was normalized for each individual separately using 18S rRNA as an endogenous control. The median and distribution of IL12B relative expression in terms of DCt for IL12B+1188 AA, AC and CC genotypes was presented in Figure 1C. No significant difference (P value = 0.165) in transcript levels of IL12B gene was observed between genotype groups representing AA ( median DC t 6S.D = 17.9861.40) and AC ( median DC t 6S.D = 18.7461.67). On the other hand, the median level of IL12B expression of AA (5 fold; P value = 0.005) and AC (9 fold; P value = 0.0001) genotypic categories were significantly lower both separately as well as collectively (7 fold, P value = 0.0001) compared to that ( median DC t 6S.D = 15.4961.25) of non-risk CC genotype group ( Figure 1C & D).
To reinforce the hypothesis that A+1188C change affects the transcription levels of IL12B gene in vivo, we assayed if there was any quantitative difference of transcripts represented by A and C alleles in patients harboring heterozygote genotype. The test for allelic-expression-imbalance (AEI) was carried out by amplification and resequencing of 167 bp fragment that covered IL12+1188 from cDNAs and gDNAs of twenty eight individuals with AC genotype. The peak heights representing A (green peak) and C (blue peak) alleles from individual electropherograms of paired samples were determined using PeakPicker software ( Figure 4A). The average A to C peak height ratio ranged from 0.42-1.05 and 0.82-1.07 ( Figure 4B) in cDNA and gDNA respectively and the comparison of these ratios was statistically significant (P value = 0.0048) estimated by nonparametric Sign test ( Figure 4C). The allele specific expression of IL12B in the peripheral blood mononuclear cells (PBMCs) of malaria patients by 39UTR polymorphism may be attributed to (i) differential transcript stability, (ii) influence of microRNAs and (iii) conformational change in RNA secondary structure.
IL12B+1188 allele-based mRNA stability assay: Inhibition of transcription by actinomycin-D In order to directly assay differences in mRNA stability between IL12B 39UTR alleles, we determined mRNA levels in cells transfected with plasmids in which 1047 bp region from IL12B 39UTR representing A or C (+1188) allele was cloned under Renilla luciferase gene with SV40 constitutive promoter in the vector pSiCHECK2 ( Figure 5A). Twenty four hours posttransfection, HepG2 cells were treated with actinomycin-D (5 mg/ml) to suppress de novo transcription. Figure 4D represents the dose response and time kinetics of actinomycin-D treatment that enabled us to select for the optimum treatment conditions for the drug. mRNA levels at 0 and 24 hours after the drug treatment were quantified by real-time PCR using IL12B and GAPDH (as endogenous control) specific primers. As shown in Figure 4E, the IL12B mRNA levels had declined 3.2 and 1.6 folds respectively for A (DC t at 0 hr = 3.5760.39 vs DC t at 24 hr = 5.2360.17) and C (DC t at 0 hr = 3.4460.68 vs DC t at 24 hr = 4.1260.47) alleles compared to their respective starting levels suggesting that A+1188C polymorphism differentially influenced IL12B mRNA stability.

Influence of microRNA regulators targeting IL12B 39UTR
To identify the putative microRNA binding sites in the 39UTR of IL12B gene we employed a consensus approach by using different miRNA target prediction softwares, viz. TargetScan, MiRanda, MicroCosm, miRBase and RegRNA. A comprehensive microRNA map of IL12B 39UTR predicted consensusly by at least three different databases was presented in the Figure 5B. Four candidate microRNAs from this list were selected for cell based reporter gene assay as they were located in the closest proximity of the polymorphic site (hsa-miR545 and hsa-miR-1284) and conserved among species (hsa-miR23a and hsa-miR-23b). Their respective seed sequences were shown in Figure 5C.
The precursor microRNA (pre-miR) sequences for each of these four regulatory RNAs were extracted from GeneCard database followed by amplification of a region encompassing 100 bp upstream and downstream of the mature miRNA and cloning in appropriate orientation into pRNAU6.1 vector. To verify whether A+1188C SNP alters the binding and regulation incurred by the above candidate microRNAs, we transfected each of these pre-miR constructs in HepG2 cells along with the 1047 bp segment of IL12B 39UTR containing either +1188A or +1188C variant cloned into the 39UTR of luciferase gene in pSiCHECK2 plasmid. Reporter plasmid co-transfected with empty miR-pRNAU6.1 vector served as the control and luciferase expression was measured 48 hours post-transfection. All the experiments were done in triplicate and repeated thrice and the change of luciferase gene expression was denoted as normalized RLU (in percentage) relative to the control. Co-transfection of hsa-miR-23a (33.33%60.08, P value = 0.0043) and hsa-miR-23b (22.22%60.04, P value = 0.0179) with pSiCHECK2-IL12B-+1188A resulted in significant reduction in Renilla luciferase expression ( Figure 5D). Similar reduction of reporter gene expression was observed when hsa-miR-23a (36.6%61.15, P value = 0.0001) and hsa-miR-23b (20%62.06, P value = 0.003) were transfected with pSiCHECK2-IL12B-+1188C ( Figure 5E). On the other hand, luciferase activity from neither pSiCHECK2-IL12B-+1188A nor pSiCHECK2-IL12B-+1188C construct was altered by hsa-miR-545 (A = 4.35%60.4 and C = 4.08%60.5) and hsa-miR-1284 (A = 12.61%60.5 and C = 4.1%60.6). Similar trend of hsa-miR-23a and hsa-miR-23b mediated reduction of luciferase activity from pSiCHECK2-IL12B-+1188 A & C constructs were observed in HCT116 and U937 cells (data not shown). Taken together, our results suggest that hsa-miR-23a and hsa-miR-23b bind to IL12B 39UTR, however, none of them exerts any allele specific influence over the IL12B expression.

Allele dependent conformational changes in RNA secondary structure
There is increasing evidence that cis-regulatory RNA structural elements within mRNAs mediate post-transcriptional gene regulation by determining several aspects of the mRNA life cycle such as stability, localization, and translational efficiency [67,68]. However, RNA structures are dynamic and evolved to adopt multiple conformations forming an ensemble that may be best described by a partition function defined as the probabilities of all possible base pairings [69]. To predict if A+1188C polymorphism provides any localized effect on the structural ensemble of IL12B mRNA, we reconstructed foldings of 101 nucleotide RNA sequences using four different bases (A, G, U and C) in the central 51 st position, keeping 50 nucleotides upstream and downstream sequences with respect to +1188 site identical. The partition function matrices computed for all four sequences using SNPfold were represented as dot-plots where probability of each base-pairing was denoted by a dot (Figure 6A to D). The pairwise Pearson correlation coefficients of A to G, A to U and A to C were 0.998, 0.925 and 0.751 respectively indicating a significant (P value = 0.013) modulation of overall structural assembly for A to C change ( Figure 6A to D).
Mountain plots were generated for 101 nucleotide sequences harboring A and C alleles ( Figure 6E & F). A mountain plot represents the secondary structure in a plot of height versus position where the height m(k) is given by the number of base pairs enclosing the base at position k. The red, green and blue curves represent minimum free energy structure, thermodynamic ensemble of RNA structures and the centroid structure respectively. The overlap between the red and green lines in Figure 6E suggested that structure produced by A allele was closer to that predicted by MFE. The blue plateau in Figure 6F indicated the central tendency of the C allele to remain in open conformation. Deviation from the base line in the position versus entropy graph of C allele shown in lower panel signified a higher randomness and positional entropy in the local structures involving C allele ( Figure 6F). This also gave rise to a higher Ensemble Diversity (ED) in the structural constellation generated by C at position 51 (ED A = 20.73 vs ED C = 34.06). Taken together, the in silico modeling of RNA secondary structure indicated an overall conformational dissimilarity in the base-pairing probabilities of the RNA thermodynamic ensemble marked by A and C alleles which could result in an alternative cis-regulatory functions mediated by IL12B 39UTR polymorphism.

Discussion
Malaria epidemiology studies have shown that, within human populations, a high degree of variation exists between individuals with respect to malaria susceptibility phenotypes, including parasite load, disease incidence, severity, and the magnitude and type of immune responses to malaria antigens [70]. Immune response to malaria has been an object of extensive investigation aiming at understanding the genetic regulation since an inappropriate immune response leads to uncontrolled parasite replication and is detrimental to host fitness [71]. A large body of the literature has demonstrated that protection against severe P. falciparum malaria is provided by distinct genetic traits, while relatively less is known about the mild events which clearly constitute a substantial fraction of global socioeconomic burden of malaria [72][73][74]. To understand the role of host genetics in the susceptibility and clinical course of mild malaria, we in this study, primarily rely on reconstruction of genetic interactions among tolllike receptors (TLRs) and pro-inflammatory cytokine genes together with an indepth functional appraisal of determinants associated with variation in blood parasite level.
Over the past several years, genetic epidemiology studies technologically have progressed from investigating few candidate markers to interrogating thousands of variants in genome-wide association studies (GWAS) while the data analysis mostly focuses on detecting single SNP effects [75,76]. Looking beyond the boundary of additive inheritance of SNPs using marker by marker approach may offer a better exercise to identify genetic determinants and biological mechanisms involved in disease etiology. Keeping this in mind, we have used MDR to reduce the dimensionality of the multilocus genetic data pertaining to immune response components to identify combinations of polymorphisms associated with the risk of high parasite load. MDR introduced by Ritchie, et al. is a computational strategy to detect and characterize gene-gene interactions that may be associated with disease susceptibility [41]. So far, MDR and its extensions such as GMDR, FAM-MDR, MDR-PDT, EMDR [77][78][79][80] have identified many interacting genetic variants underlying a wide variety of complex human diseases such as Alzheimer's disease, asthma, atrial fibrillation, autism, bladder cancer, hypertension, nicotine dependency, prostate cancer, schizophrenia, sporadic breast cancer and type II diabetes [81][82][83][84][85][86][87][88][89]. In our study single SNP and multifactor dimensionality reduction analyses together identified three genetic loci which served as important factors in controlling blood parasitemia in mild malaria. These loci include TLR4Thr399Ile, TNF-1031 and IL12B+1188. Two additional loci viz. LTA80 and TLR9P545P appear in combination with TLR4Thr399Ile and IL12B+1188 suggesting that the formers fail to show statistically significant association in independent analyses because of their minor effects  (2) the effect of miRNAs. Co-transfection of the empty pRNAU6.1 (+) vector with pSiCHECK2-IL12B+1188 construct was set as 100% and reductions in luciferase expression in presence of four miRNAs were measured in relation to this. hsa-miR-23a and hsa-miR-23b resulted in significant reduction in luciferase activities for both pSiCHECK2-IL12B+1188A and pSiCHECK2-IL12B+1188C constructs. Statistical significance was measured with t-test. (*) indicates the P values and percentage reduction to be statistically significant. doi:10.1371/journal.pone.0046441.g005 on trait variation (Table 3) [90]. Notably, the minor allele frequencies of both LTA80 and TLR9P545P were approximately 40% in our samples [21]. The presence of IL12B+1188 in five of six models and elevation of odds ratio of TLR4Thr399Ile, TNF-1031 and LTA80 in combination with IL12B+1188 are suggestive of a pivotal genetic influence of IL12B on malaria phenotype. This is supported by the fact that different malaria phenotypes such as hyperparasitemia, severe malarial anemia or cerebral malaria are often found to be associated with IL12B gene [34,[91][92][93][94] while the genetic alliance between parasitemia with LTA is seldom unequivocal [95][96][97]. A study conducted among 198 individuals belonging to 34 families living in Burkina Faso using the pedigreebased generalized multifactor dimensionality reduction approach, identified statistical interactions among immune genes including IL12B 39 untranslated region, IL12Bpro, LTA+80 for mild malaria, maximum parasitemia or asymptomatic parasitemia [98] underscoring the usefulness of modeling gene-gene interaction in genetic dissection of complex diseases ( Table 3).
The association of high parasitemia with IL12B+1188 (Table 1) and significantly low expression of the cytokine in individuals harboring risk genotypes (Figure 1) indicate that the susceptible individuals are deficient in controlling parasite replication due to inadequate level of IL12B transcript in mild malaria patients. Critical role of IL12B in malaria has been demonstrated by other studies that show IL12 production is inversely associated with disease severity in human malaria [94,99,100] and the molecule is extremely effective in correcting malarial anemia in murine model [101,102]. Reduced expression of IL12B and LTA under TLR4 risk genotype background may be attributed to altered TLR4 signaling ( Figure 2). To the best of our knowledge, this is the first report that provides functional validation of statistical epistatic models demonstrating in vivo suppression of cytokine gene expression due to genetic deficiency of TLR4 signaling. The attenuated TLR4 signaling by the variant receptors may be ascribed to altered distribution of electrostatic surface potential as shown by homology modeling in our previous study [21]. However, it should be borne in mind, a comprehensive profiling of transcriptome and network analyses are necessary to identify additional genetic partners involved in the pathway.
Pinpointing phenotypically causal variants from a large fraction of SNPs that show statistical association to diseases remains a major challenge. Only a handful of disease-associated variants occur at non-synonymous (nsSNPs) sites while a majority are located on the non-coding genomic regions (ncSNPs) adjacent to protein coding genes suggesting that the latter may be involved in transcriptional regulation. Global expression quantitative trait loci (eQTL) analyses in yeast, mice, and humans have detected significant levels of cis and trans-eQTL that simultaneously regulate a large fraction of the transcriptome [103][104][105][106][107][108]. We previously reported that the polymorphism located at TNF-1031 displayed significant association with peripheral parasite load [21]. Here we showed that any variation from the TNF-a wild-type haplotype resulted in reduced promoter efficiency and the maximal reduction was observed for the haplotype that was altered at 21031 locus. The dysfunctioning of TNF-a promoter due to 21031 polymorphism may be attributed to allelic differences in binding of transcription factors (TFs). A number of transcription factor binding sites (TFBSs) were abolished in presence of the variant allele (C) which assumed a sterically inaccessible locked-in conformation as depicted in the bioinformatic outputs ( Figure 3D-F). Chromatin immunoprecipitation or electrophoretic mobility shift assays (EMSA) are necessary to validate in-silico prediction of altered binding affinity between putative TFs and promoter SNP. Reduction of transcriptional activity due to single nucleotide variation in TNF-a promoter has been reported for the polymorphisms located at 2863 [109], 2376 [110] and 2238 [111].
Allelic-expression-imbalance (AEI) phenomenon which serves as an integrative quantitative measure of any and all cis-acting regulatory variants [112] is highly context-specific [113] and widespread not only in humans and mice but in most organisms [114][115][116][117][118][119][120]. It can affect the transcriptome by altering stability [121], processing efficiency [122], isoform expression of mRNAs [123,124] and inducing epigenetic changes [125,126]. The unequal expression of IL12B between marker genotypes and allelic imbalance detected in malaria patients' PBMC samples strongly suggest that IL12B+1188 polymorphism executes a cisregulatory function. To interpret the probable mechanisms of IL12B+1188 mediated AEI, we have investigated allele specific transcript stability, differential interaction of IL12B 39UTR with candidate microRNAs and allele-dependent modulation of RNA sub-optimal structures. The remarkable modulation of local RNA secondary structures due to A to C transversion observed in this study may affect the transcript stability directly or through the association of RNA-binding proteins (RBPs) indirectly [69]. Many RNA binding proteins have both sequence-specificity and RNAsecondary structure binding preferences and are known to coregulate functionally related transcripts [127]. From SNP-targeted studies, it has been estimated that 20% of the measured transcripts show 1.5 fold differences between alleles while 30% show 1.2 fold differences [118,128]. However, it is still obscure how small allelic imbalance may cause a phenotype and how accurate the bioinformatic predictions are with actual experimental evidences.
Recent studies have also shown that 39UTRs contain recognition motifs for microRNAs (miRNAs) which play important generegulatory roles by pairing to their target mRNAs to direct posttranscriptional repression. The most accurate predictors of miRNA target sites relies on conserved matches to the ,7 bp seed region near the 59end of the microRNA and also on the accessibility of target sites. The four candidate microRNAs (hsa-miR-545, hsa-miR-1284, hsa-miR-23a and hsa-miR-23b) predicted in common from multiple databases target IL12B 39UTR sequence located proximally to A+1188C. Our data showed hsa-miR23a and hsa-miR-23b interacts with IL12B 39UTR to repress gene expression without any allelic bias. This may be attributed to our reductionist approach to candidate microRNA prediction and limited ability of the current prediction algorithms that calculate target efficacy based on interactions between the mRNA with itself and the mRNA with a miRNA. This is of particular importance because mRNA-miRNA interaction occurs in a complex cellular environment in which mRNAs and miRNAs are likely to be bound by cellular RNA-binding proteins, which are currently impossible to account for in silico. Our findings exemplify that diversity of posttranscriptional gene regulation may extend beyond microRNAs and underscores the importance of characterization of structural cis-regulatory elements and their interaction partners in the context of mRNA stability. Recently, there has been renewed interest in identifying the effects of disease-associated noncoding SNPs on changes in RNA structure [129,130] which may regulate gene expression at virtually every possible stage ranging from local chromatin remodeling to mRNA translation [131][132][133]. Taken together, our findings convincingly demonstrated inadequate gene expression of pro-inflammatory cytokines, namely IL12B and TNF-a attributable to cis-regulatory polymorphisms together with loss of efficacious TLR4 mediated signaling due to nsSNP contribute to uncontrolled parasite growth in P.falciparum malaria in an epistatic manner. In the study of genetics of complex diseases, to determine which of the multitude of variants carried by an individual are responsible for a given phenotype remains a massive task. Our study emphasizes the need to evaluate genegene interaction and biological credibility of hundreds of common non-coding variants with low effect size as powerful complementary research strategies to illuminate the genetic architecture of common diseases.