Abstract
Non-coding genetic variants play an important role in driving susceptibility to complex diseases but their characterization remains challenging. Here, we employed a novel approach to interrogate the genetic risk of such polymorphisms in a more systematic way by targeting specific regulatory regions relevant for the phenotype studied. We applied this method to meningococcal disease susceptibility, using the DNA binding pattern of RELA – a NF-kB subunit, master regulator of the response to infection – under bacterial stimuli in nasopharyngeal epithelial cells. We designed a custom panel to cover these RELA binding sites and used it for targeted sequencing in cases and controls. Variant calling and association analysis were performed followed by validation of candidate polymorphisms by genotyping in three independent cohorts. We identified two new polymorphisms, rs4823231 and rs11913168, showing signs of association with meningococcal disease susceptibility. In addition, using our genomic data as well as publicly available resources, we found evidences for these SNPs to have potential regulatory effects on ATXN10 and LIF genes respectively. The variants and related candidate genes are relevant for infectious diseases and may have important contribution for meningococcal disease pathology. Finally, we described a novel genetic association approach that could be applied to other phenotypes.
Similar content being viewed by others
Introduction
A vast majority of Single Nucleotide Polymorphism (SNP) associated with susceptibility to complex diseases identified through Genome Wide Association Studies (GWAS) are located in non-coding regions of the genome1. They have been hypothesized to affect gene regulation, notably via variation in transcription factor binding2,3. Despite the growing interest and resources available to study these polymorphisms, understanding their functional effect remains challenging for several reasons: (i) most associated SNPs are still identified through genome-wide genotyping arrays, which does not allow for all variants to be investigated but only tag SNPs linking a locus to a change in gene expression, (ii) studying the right cell type in the right environment is necessary to uncover the mechanism of action of a variant because gene expression and transcription factor binding varies across tissues and conditions4. In an attempt to address these challenges, we employed a reverse genetic approach to identify regulatory variants involved in the innate immune response.
We have recently identified genome-wide binding of RELA, one of the Nuclear Factor kappa B (NF-kB) members involved in the response to microbes, as well as gene expression data following microbial stimuli in nasopharyngeal epithelial cells5. In addition we also investigated epigenetic changes following a potent gram negative bacterial endotoxin, Lipopolysaccharide (LPS), in the same cells6.
We concluded that some of the potential regulatory regions identified in our previous study will be relevant in mounting an immune response against infectious pathogens with the following characteristics: (i) airborne transmitted, (ii) able to infect human epithelial cells, (iii) able to bind some of the receptors targeted in our previous study, and lastly (iv) shown to have a host-genetic susceptibility component. As such we identified Meningococcal Disease (MD) as a relevant example, complying with all the requirements mentioned above.
MD is a severe infection resulting in potentially lethal meningitis and sepsis. It is caused by a gram negative bacterium, Neisseria meningitides, which is transmitted through respiratory secretions. The bacterium colonizes the nasopharynx before crossing the epithelium to cause the serious invasive disease7. Interestingly, asymptomatic carriage in the nasopharyngeal mucosae is relatively common with the bacteria being detected in about 10% of the population8.
Neisseria meningitides has been shown to bind several pattern recognition receptors, TLR2, TLR4, TLR9 and NOD receptors9,10, leading to the activation of downstream signaling pathways. One of the main TF activated is the master regulator of the innate immune system, NF-kB11.
Previous studies have demonstrated that host genetic make-up is a risk factor for MD12 and a number of polymorphisms have been associated with susceptibility to the disease, notably in innate immunity genes13. Thus, host-pathogen interactions are decisive in the development of the disease, notably at the nasopharynx epithelium where epithelial cells are critical in detecting pathogens and organizing an efficient immune response14.
Finally, our group has been involved in previous genome wide association studies (GWAS) for MD susceptibility15,16, therefore we have access to well-characterized cohorts for this disease.
Briefly, our approach consisted of identifying regulatory regions in response to bacterial stimulation of pharyngeal epithelial cells which were then used to perform targeted sequencing in cohorts of healthy individuals and MD patients followed by further validation in three European cohorts. This strategy allowed us to identify two novel SNPs, rs4823231 (P-value = 9.58 × 10−5, OR = 0.73) and rs11913168 (P-value = 3.46 × 10−3, OR = 0.77) showing association with genetic susceptibility to MD.
Results
Selection of regulatory regions relevant for airborne bacterial infection
We have previously identified RELA genome-wide binding sites as well as gene expression in Detroit 562 cells in response to different microbial stimuli5 and in FaDu cells under LPS stimulation. Both of these lines are pharyngeal epithelial cells. In addition, we have identified H3K27ac changes following LPS stimulus in both cell lines6. The regions detected were particularly relevant for infectious respiratory diseases, especially bacterial infection, and were selected for targeted sequencing (further details on region selection can be found in Methods).
In total, 9,551 genomic regions were selected, covering 9,943,597 basepairs (bp) (see Supplementary Fig. S1A for an example of the regions covered). Expectedly, Gene Ontology analysis on the nearby genes revealed that these regions were highly enriched for the “immune response” as well as “response to other organism”, which demonstrates biological relevance to infection (Supplementary Fig. S1B).
Using the Nimblegen technology, these regions were used to design a custom probes set as baits to specifically capture targeted fragments for sequencing. The capture designed consisted of 16,784 probes (Supplementary Table S1) resulting in 8,923,013 bp of the initial targets submitted covered (89.7%) and an estimated 9,569,494 bp covered by targeted sequencing (96.2%). Only 37 of the regions submitted could not be included (mainly due to sequence repeat) and showed no coverage in the probes designed (Supplementary Table S2).
Targeted sequencing in meningococcal disease cases and controls
First, we performed targeted sequencing of the selected regions in 238 MD cases and 237 controls from Western Europe. Several hybridization reactions were tested for enrichment by qPCR and showed satisfactory results (Supplementary Fig. S2A). Moreover, the sequencing reads were also highly specific for the regions covered in the capture (Supplementary Fig. S2B). This indicates that the targeted sequencing experiments were successful. Depth of coverage was satisfying across all samples (Supplementary Fig. S2C) with an average of 61.8x and the majority of the capture regions (16,357 out of 16,784) were characterized with sufficient coverage and adequate mapping.
Variant calling was performed using GATK Unified Genotyper, following the best practices17 and identified 92,929 Single Nucleotide Polymorphisms (SNP). Stringent filtering steps were then applied (more details in Methods) and finally 30,542 SNPs remained for the association analysis (Supplementary Fig. S3A). Samples were also filtered (Supplementary Fig. S3B,C), removing the samples with first degree relatives (N = 10), PCA outliers (N = 5) and samples with insufficient call rates (<95%, N = 0). Ultimately, 460 samples remained to be considered for the association study (235 cases and 225 controls).
Genetic association analysis was performed using a logistic regression model and adjusting for population heterogeneity by correcting the P-value with the genomic inflation factor (GC = 1.049). The top 15 SNPs showing the most significant association are reported in Table 1 (see Supplementary Table S3 for the full list).
Validation of top SNPs and Meta-analysis
Selected SNPs from the discovery phase were genotyped either by mass array with the Sequenom technology or by GWAS array using the previously generated data15. They were typed in 3 validation cohorts from Western Europe, Spain and the UK.
Using the same stringent criteria as the discovery stage, quality control filtering was done on the SNPs as well as the samples. Independent association analysis was performed separately in the three validation cohorts using the same method as in the discovery cohort. A meta-analysis including all four cohorts was then computed (see Supplementary Tables S4 and S5 for the association results from the Sequenom and GWAS array genotyping respectively).
The top SNPs validated with a meta-analysis P-value < 0.05 are reported in Table 2. We then focused on the top two SNPs – rs4823231 and rs11913168 – in the rest of the study as they showed much lower P-values.
We identified the most significant association in the meta-analysis for rs4823231, the P-value for association was 9.6 × 10−5. This variant showed a strong protective effect of the minor allele (T) in the discovery cohort with an odds ratio (OR) of 0.56 and P-value of 1.12 × 10−3. Although the associations in each validation cohorts did not reach significance, they all showed a similar trend with comparable P-values, frequencies as well as ORs (Table 2). The associations of the risk allele (T) across stages were in the same protective direction captured with the meta-analysis OR of 0.73 (Table 2).
The SNP rs4823231 is located in an intron of the ATXN10 gene on chromosome 22. In order to interrogate potential regulatory effects of this SNP, we looked for it in the HaploReg data base. The region where the SNP is located was identified as an enhancer in 14 tissues based on histone marks data from the Roadmap Epigenetic18. Our ChIP-seq for H3K27ac in two pharyngeal epithelial cell lines also revealed high deposition of this mark at the SNP location (Fig. 1A). This suggest that this region is active in this relevant cell type. Moreover, we noticed that the SNP is located very near a RELA binding site induced by microbial stimulation in these same cell lines (Fig. 1A).
Genomic data for rs4823231 and motif disruption. (A) UCSC genome browser view of the RELA and H3K27ac ChIP-seq data generated in two nasopharyngeal epithelial cell lines. (B) Potential disruption of transcription factor motifs induced by rs4823231. The SNP position is indicated by the red boxes. Ref: reference, Alt: alternate.
Our panel was designed to include RELA binding sites located near a differentially expressed gene or a region of LPS-induced H3K27ac changes. In this case, the FBLN1 gene located about 100 Kb upstream of rs4823231 was found to be up-regulated in one RNA-seq experiment after M tri-DAP stimulation of Detroit 562 cells and LPS stimulation of FaDu cells (Supplementary Fig. S4A). However, the duplicate RNA-seq experiment performed later on and RT-qPCR on these conditions did not replicate this change in expression (Supplementary Fig. S4B). In fact, results for this SNP from HaploReg revealed that rs4823231 is an eQTL associated with the expression of the ATXN10 gene in many tissues (based on GTEx data). This association was particularly significant in the oesophagus mucosae tissue (P-value = 2.6 × 10−18), which is the closest tissue to our model (Supplementary Table S6). A closer look at the GTEx data19 from which this results were extracted, we observed that individuals carrying the genotype T/T have a higher level of ATXN10 expressed in their oesophagus mucosae compared to those with the genotype A/A (Supplementary Fig. S4C). When we looked at the expression of ATXN10 in our data, we found that this gene was highly expressed in both Detroit 562 and FaDu cells (average FPKM~325 in both cell lines in RNA-seq experiments and easily detectable by RT-qPCR – Supplementary Fig. S4D) but there was no change in expression after microbial stimulation. A further inspection of the Human Protein Atlas20 revealed protein expression in most tissues available in this database (Supplementary Fig. S4E).
As rs4823231 is not located inside a RELA binding sites identified in our study, it is less likely to have an effect on RELA binding itself. However, it is located very near and could therefore influence binding of a RELA co-factor. HaploReg returned this SNP to potentially alter 4 transcription factor motifs, DMRT4, DMRT7, E2F and FOXA. In addition, RegulomeDB tool also indicated rs4823231 to disrupt two factors motifs, DMRT4 and IRX3 (Supplementary Fig. S5). Concerning the DMRT4 motif, the SNP is located at a position that is less conserved where any of the A, T or C bases can be found thus it is less likely to have an effect on this factor binding. This is further captured by the score given to the reference and alternate alleles in HaploReg which is very similar. The same observation is true for DMRT7 motifs where the variant occurs at a variable base position. However, IRX3 motif alteration reported by regulomeDB is interesting as rs4823231 is located at a position where a base C or T is present, suggesting that the major allele A would disrupt the motif (Fig. 1B). Moreover, the E2F and FOXA motifs match identified in HaploReg are also of interest as they showed that the major allele A (base T highlighted in Fig. 1B as the motif matches the reverse complement strand) could create a motif for these factors to bind that location.
The second top SNP more significantly associated with MD susceptibility in our study was rs11913168. This variant showed association with a P-value of 8.69 × 10−3 in the Western Europe discovery cohort which was replicated in the British validation cohort with a P-value of 3.7 × 10−2. The Spanish cohort showed a weak association (P-value = 0.53) while the Western European cohort revealed sign of association but not very strong (P-value = 0.21). Nonetheless, all cohorts showed a similar trend leading to a meta-analysis P-value of 3.5 × 10−3 and an OR of 0.77 suggesting a protective effect of the minor allele A.
The SNP rs11913168 is located in an intergenic region on chromosome 22 (Fig. 2). According to HaploReg, this SNP is under sequence constraint (by GERP and SiPhy methods) and this is also observed in the UCSC genome browser with positive PhyloP score at this particular location (Supplementary Fig. S6). HaploReg reported that the SNP is located in a region with enhancer mark in 14 cell types, DNase sites in 29 tissues and 10 bound proteins which could be seen in the UCSC genome browser (Supplementary Fig. S6). In our genomic data, we also found that this SNP was located in a region with H3K27ac marks in both Detroit 562 and FaDu cells (Fig. 2A). Notably, in FaDu cells, this region exhibited increased H3K27ac signal after LPS stimulation as well as a bigger region without acetylation signal (Fig. 2B). This can be linked to nucleosome remodeling leading to a larger nucleosome-depleted region and transcription factor binding, characteristic of an active enhancer21. In addition, rs11913168 was located within a RELA binding site in both cell lines.
Genomic data for rs11913168. RELA and H3K27ac ChIP-seq data following LPS stimulation of Detroit 562 and FaDu cells are represented. The SNP of interest is indicated with a black line. The lower panel (B) consists of a magnification of the dashed box in the upper panel (A).
The SNP rs11913168 is situated between two genes, HORMAD2 and LIF (Fig. 2). Interestingly, LIF was found to be up-regulated following stimulation of Detroit 562 cells with microbial stimuli as well as stimulation of FaDu cells with LPS (Fig. 3A). This was further validated experimentally in FaDu cells after LPS treatment (Fig. 3B). Moreover, using publicly available Hi-C data from Jin et al.22 we found evidences of chromatin interaction between the region containing rs11913168 and regions around the LIF gene (Supplementary Fig. S7A). Taking into consideration these data, we hypothesized that the region where the SNP is located is an enhancer of LIF and that chromatin looping would bring this regulatory region and the gene close to each other. We searched available resources listing known NF-kB targets (see Methods) but we could not find LIF on any of these lists. Thus we investigated it in our system. We pre-treated the FaDu cells with a widely used NF-kB inhibitor, BAY 11–7082, before LPS stimulation to inhibit RELA. Indeed, we observed that this pre-treatment abolished RELA activation (Supplementary Fig. S7B) as well as up-regulation of known NF-kB target genes TNF and NFKBIA (Supplementary Fig. S7C) following treatment. In the same conditions, LIF up-regulation following LPS stimulation was also completely suppressed (Fig. 3B), suggesting that this gene is regulated by NF-kB. Therefore, binding of RELA at the regulatory region where rs11913168 is located could be potentially regulating LIF expression.
LIF expression and NF-kB regulation. (A) The bar graph consists of the average log2FC between two RNA-seq experiments following 4 stimuli in Detroit 562 cells (green – see legend) or LPS stimulation in FaDu cells (orange) compared to the control condition without any treatment. (B) LIF expression following LPS stimulation in FaDu cells without (DMSO) or with pre-treatment with BAY 11–7082 (BAY inhibitor) compared to the control condition without treatment. Each dot represent the results from one experiment and the black dash is the mean between duplicates. Significance was tested with a 2-sample T-test for which the P-value is reported.
Based on our ChIP-seq data, we noticed that the Detroit 562 cells that we used for our genomic analyses were heterozygous (G/A) for rs11913168 and this was further validated by Sanger sequencing (Supplementary Fig. S8A). We investigated any allelic bias for one allele or the other in our RELA and H3K27ac data. By counting the reads carrying the G or A allele, we noticed that there were more reads covering the G allele. After ensuring that there was no strand bias (Supplementary Fig. S8B), we tested this statistically and found that the G allele was indeed significantly enriched in the RELA ChIP-seq (P-value = 2.8 × 10−13) and H3K27ac ChIP-seq in control (P-value = 2.9 × 10−3) and LPS (P-value = 0.018) conditions (Table 3). This suggests an allelic bias where the transcription binds preferentially to the base G at that position than to the base A and this histone mark follows the same trend. The same analysis on input DNA also revealed more reads containing the base G compared to A, although this did not reach statistical significance probably due to the very small number of reads detected in this data-set.
Discussion
In this study, we used genomic data generated in nasopharyngeal epithelial cells stimulated with microbial components, modelling airborne bacterial infection, to design a capture probes panel for targeted sequencing. Association analysis and further validation by genotyping in three independent cohorts revealed two new SNPs associated with MD susceptibility in European populations: rs4823231 and rs11913168. Although these two SNPs showed some evidence of association, the P-value from the meta-analysis would not withstand adjustment for multiple testing if applied.
The SNP with the smallest P-value in this study was rs4823231 (P-value = 9.6 × 10−5), showing a protective effect of the minor allele (T) (OR = 0.73). The risk effects of this variant across the different cohorts was variable, especially between the discovery and validation cohorts. It likely to be due to the variation in allele frequencies in controls in the discovery cohort against the three validation cohorts. Yet, these percentages remained coherent and in line with the known frequencies in European populations reported in public data bases. Indeed, dbSNP data base reported the following frequencies in the European population: 82.5% A and 17.5% T in the HapMap data and 83.4% A and 16.6% T in the 1000 Genome project data. Moreover, all three validation cohorts showed the same trend with the minor allele T being less common in controls than in cases, with similar effect sizes and P-values. The limited sample size of these cohorts could explain why the P-values did not reach a more significant level. This is especially true for the Spanish and the British cohorts as they contained a small number of cases and controls respectively. The data base dbSNP also revealed differences in terms of allele frequencies across populations. Indeed, Africans were much less likely to carry the T allele, with HapMap and the 1000 Genome Project data reporting frequencies of 0.8% and 1.8% respectively. Our analysis showed that the T allele has a protective effect, thus the fact that this allele is uncommon in this population may contribute to the epidemic scale of meningococcal disease in Africa7.
One of the candidate gene on which rs4823231 could have an effect is ATXN10 as it was reported by the GTEx consortium19 that this genetic variant is an eQTL for this gene. The ATXN10 gene was first described in neurons where it is highly expressed and involved in growth and apoptosis23. However, this gene has been later shown to be expressed in many tissues, as supported by our search on the Human Protein Atlas. Recently, research on HeLa cells showed that decreased expression of ATXN10 leaded to defect in cytokinesis24, thus further supporting its role in the cell cycle in various cell types. There is preliminary evidence suggesting a link between ATXN10 and the NF-kB pathway. Warner et al. conducted an RNAi screen in HEK cells to identify genes involved in the NF-kB response to NOD2 stimulation, an intracellular PRR targeted by certain motifs of bacterial peptidoglycan. They reported that knock-down of ATXN10 leaded to a reduced NF-kB activity measured by a decrease in reporter expression following stimulation of NOD225. Moreover, they further showed that the same knock-down led to a reduced IL-8 production following stimulation26. The NF-kB pathway plays an important role in survival and apoptosis11, thus this could explain how ATXN10 is also involved in these processes. Overall, these evidences suggest that high level of expression of this gene is beneficial for cell division and survival, especially following microbial detection. This is in line with the eQTL effect seen for rs4823231 and our association results. Indeed, the T allele for this SNP was associated with a higher level of ATXN10 and this allele also showed a protective effect in our analysis with the healthy controls being more likely to carry this allele than the MD patients (see Supplementary Fig. S9 for a schematic). These changes in expression could have critical consequences on tissue integrity due to the role of ATXN10 and could in turn influence bacterial infection outcome. However, the mechanism by which rs4823231 leads to changes in ATXN10 expression remains elusive. We can hypothesize that it could be due to variation in TF binding at this location as we found that the SNP potentially impacted motifs for IRX3, E2F and FOXA. The possible effect on E2F binding is of particular interest as these factors are involved in the cell cycle and apoptosis. Moreover, some members of the E2F family are repressors27 which could explain that recruitment at this locus via the binding motif created by the major allele A would result in repression of ATXN10 expression. In addition one of the E2F member, E2F1, has been shown to interact with RELA following LPS stimulation28.
The second most significant SNP, rs11913168 (P-value = 3.46 × 10−3), revealed a protective effect of the A allele compared to the more common G allele with an odds ratio of 0.77. The British and Western European validations cohorts were in line with the discovery cohort, especially the British cohort for which a strong association was reported. The Western European validation cohort showed a similar effect size although it did not reach the same level of significance, probably due to the limited number of cases available in this population. However, the effect of the SNP was not validated in the Spanish cohort which may be due to variation in frequencies across populations, with the minor allele (A) being less common in Spanish individuals. This is further supported by in silico down-sampling of the UK cohort to similar numbers of cases and controls than the Spanish cohort revealing an association P-value for rs11913168 of 0.061 which is much more significant than the association in the Spanish cohort. Overall, the frequencies observed in our study were in line with those reported in dbSNP where the HapMap data showed 89.9% of the major allele G and 10.1% of the minor allele A while the 1000 Genome project reported 89.4% and 10.6% for allele G and A respectively in Europeans.
Concerning rs11913168, a potential gene regulated by this polymorphism could be the LIF gene located around 30 Kb downstream. In addition to our genomic data as well as chromatin interaction evidence, we found this gene to be up-regulated following microbial stimuli in epithelial cells. Similarly, an earlier report showed an increased LIF expression in rat trachea after treatment with LPS29. This suggests that LIF gene product plays a role in defense against pathogens. Indeed, LIF encodes a cytokine of the IL-6 superfamily, it has various functions, including signaling to immune cells such as T-cells and monocytes30 which is particularly important in the case of infection. Moreover, we showed that LIF looks to be regulated by the NF-kB pathway and we have data suggesting that RELA binding at rs11913168 would be influenced by the nucleotide present at this particular location. Taking this into account, we can hypothesize that the variation in RELA binding at the rs11913168 locus could influence the expression of the LIF gene downstream, leading to fluctuations in this cytokine’s secretion. This would in turn have consequences on the recruitment and communication with specialized immune cells needed to fight the pathogen (see Supplementary Fig. S9 for a schematic). However, differential RELA binding as well as its effect on LIF regulation warrant further experimental validation. Notably, the allelic bias observed in the ChIP-seq data in Detroit 562 cells should be further quantified. These cells consist of a cancer cell line and could contain more than two copies of each chromosome, therefore the bias observed may not reflect the actual binding preference of the transcription factor.
From previous GWAS performed on the same cohorts, several SNPs were identified as associated with meningococcal susceptibility, they were located close to the CFH and CFHR3 genes15,16. Although two regions at the CFH locus were included in the custom capture used in this study, they did not overlap with the variants identified in the GWAS (Supplementary Fig. S10). As a result, the previously associated SNPs were not investigated in our study. Moreover, we observe that the RELA and H3K27ac ChIP-seq signal around these regions is low in both epithelial cell lines studied. This is consistent with the expression of CFH which is produced primarily by the liver31. Other variants have been associated with MD susceptibility13, however, a closer observation revealed that none of them overlapped the capture designed.
Conclusion
In this study, we employed a novel method to discover new regulatory polymorphisms driving genetic susceptibility to MD. By using genomic data generated in a relevant cell type – nasopharyngeal epithelial cells – and under relevant stimulation – bacterial stimuli – we designed a capture for targeted sequencing of cohorts of MD patients and healthy controls. Following validation, we were able to identify two new putative regulatory variants that showed association with MD susceptibility.
Methods
Custom capture design
The capture used in this study was based on data generated previously in Detroit 562 cells5,6, available in the NCBI’s Gene Expression Omnibus repository under the following accession number: (i) GSE91018 for RELA ChIP-seq data under microbial stimuli, (ii) GSE91019 for the RNA-seq data following microbial stimuli and (iii) GSE104635 for the H3K27ac ChIP-seq data under LPS stimulation. In addition, RELA ChIP-seq, RNA-seq and H3K27ac ChIP-seq data generated in the FaDu cells treated with LPS is available on demand.
Briefly, the regions selected to design the capture were the following: (i) RELA binding sites located within 200 Kb of a differentially regulated genes following stimulation with bacterial stimuli (LPS, Pam2CSK4 or M tri-DAP) or TNFa in Detroit 562 cells, (ii) H3K27ac regions that showed a signal consistently higher under LPS stimulation compared to control, (iii) RELA binding sites within 200 Kb of differentially expressed genes under LPS treatment in FaDu cells and (iv) H3K27ac regions with a higher signal in LPS than Control ChIP-seq in FaDu cells. At the time of the capture design, some duplicates of these experiments were not available, thus some of the gene expression data was not replicated, as detailed in the results section for FBLN1.
The coordinates of the selected regions were submitted to Nimbelgen via their website for the design of the custom capture and the hybridization probes were then ordered.
Targeted sequencing
Genomic DNA libraries from 475 Western European samples (238 cases and 237 controls) were prepared previously with the NEBNext Library Prep kit (New England Biolabs) according to the manufacturer’s instructions. After library building, each sample was loaded on a DNA 1000 chip and run on the Bioanalyzer (Agilent) for quality check and quantification. Based on this concentration, equal amounts (58.82 ng) of DNA from each of the 17 samples (either 8 controls and 9 cases or 9 controls and 8 cases) that made up a multiplexed capture reaction were pooled together, resulting in 1 µg of total DNA/Capture pool. Hybridization was then performed according to the SeqCap EZ system (Nimbelgen) protocol (SeqCap EZ Library SR User’s Guide v5.1) as follow. The pooled DNA was mixed with 5 µg of COT DNA, 1 nmol of universal oligo (MP-HE-1) and 1 nmol of Terminator oligo pool matching our custom indexes. The mix was dried up in a desiccator at 60 °C for 20 to 30 minutes and resuspended in 7.5 µl of 2X SC Hybridization Buffer and 3 µl of SC Hybridization component A. DNA was denatured by incubating the tubes for 10 minutes at 95 °C before being transferred to a 4.5 µl aliquot of biotinylated custom capture and incubated at 47 °C for 50–52 hours. The samples were then washed according to the manufacturer’s instructions. Post-capture PCR was performed with the Phusion High Fidelity master mix (New England Biolabs). Briefly, 150 uL of the 2X master mix was added to the 50 µl of beads from the capture, together with 3 µl of each primers (MP_PE_POST 1 and 2) and 94 µl of water. Six aliquots of 50 µl were loaded into a PCR strip and the following program was applied: 1 min at 98 °C + 18 × (30 seconds at 98 °C + 30 seconds at 60 °C + 30 seconds at 72 °C) + 1 minutes at 72 °C. The samples were then purified with 0.9X Ampure beads and loaded on a DNA 1000 chip to be run on the Bioanalyzer for quality check and quantification. Hybridization reactions were diluted to 10 nM and 2 reactions were pooled together in equal amount for sequencing on one lane. Libraries were sequenced on an Illumina Hi-seq (paired end 2 × 101 reads) lane. The sequences of the primers and oligos used can be found in Supplementary Table S7.
Hybridization quality control by qPCR
Selected hybridization reactions were tested by qPCR for enrichment of negative and positive loci as follows. Pooled genomic DNA libraries (Pool) as well as post-capture PCR products (Post Capture) were diluted to 5 ng/µl and 1 µl was used per reaction. Primers for 2 positive loci (NFKBIA and TNF) and 2 negative loci (ACTB and chr12D) were used, sequences can be found below. Quantitative PCR with LightCycler 480 SYBR Green I Master mix (Roche) was performed according to the manufacturers’ instructions in triplicates. The following PCR program was used: 5 minutes at 95 °C; 45 × (10 seconds at 95 °C, 1 minutes at 65 °C, 30 seconds at 72 °C). The fold enrichment for each locus was calculated using the following formula: 2^(average Ct value of the Pool sample – average Ct value of the Post Capture sample).
Discovery phase: variant calling and association analysis
Sequencing reads were mapped to the targeted regions of the human genome build hg19 using bwa32, followed by local realignment and base quality score recalibration of the reads. The following analyses were then performed with GATK 3.517. Depth of coverage as well as interval quality control were run on each sample with the modules DepthOfCoverage and DiagnoseTargets respectively. SNPs were called in all samples from the discovery cohort together using the Unified Genotyper module. To obtain high quality variants, SNPs with: (i) a quality score < 50, or (ii) QD < 2, or (iii) FS > 60, or (iv) MQ < 40, or (v) haplotype score > 13, or (vi) MQRankSum < −12.5, or (vii) ReadPosRankSum < −8 were further removed.
The rest of the analyses was performed in PLINK v1.0733. Mono-allelic SNPs were removed and further filtering was performed: SNPs with (i) minor allele frequency (MAF) < 0.5%, or (ii) genotyping rate < 95%, or (iii) deviation from Hardy-Weinberg equilibrium in controls with a P-value < 0.0001, or (iv) significantly different genotyping rate between case and controls were removed (Supplementary Fig. S3).
Sample filtering was performed as follow: (i) samples with first degree relatives identified with identity by descent analysis were removed (keeping only the sample with the highest genotyping rate per pair), (ii) Principal Component Analysis (PCA) outliers were removed based on PCA1 and PCA2 (PCA analysis was performed using the GWAS data available to us, extracted for the 475 samples of the discovery cohort), and (iii) samples with a call rate < 95% were filtered out (Supplementary Fig. S3).
Association analysis was performed on the final set of variants (30,542 SNPs) between 235 case and 225 control samples remaining using a logistic regression model in PLINK (–logistic). In addition, the option–adjust was used to correct for population heterogeneity with the genomic inflation factor (GC). For this population, GC reported was 1.049 thus the P-values reported from the association were corrected with this factor (GC P-values).
Validation phase: SNP selection and mass array genotyping
Based on the results from the discovery phase, SNPs were selected for validation if: (i) their association P-value < 0.005, (ii) they were not in linkage disequilibrium with other SNPs (r2 < 0.9), keeping one SNP per LD block, (iii) they were not present on the GWAS array used to type the validation cohorts, and (iv) their MAFs were not too divergent from databases or too low to be replicated (<4%). Finally, 42 SNPs were selected for validation and submitted to design a custom Sequenom Multiplex MassArray assay. One SNP failed to be designed so 41 SNPs were left and genotyped in three validation cohorts according to the manufacturer’s instructions.
Stringent thresholds of 95% call rate for SNPs and 95% genotyping rate for samples were applied to filter the data. Fourteen SNPs were filtered out at this step; 27 SNPs remained for the association analysis. The samples consisted of 3 cohorts containing the following number of samples: 501 in Western European (Dutch together with Austrian) cohort (350 cases and 332 controls), 1182 in the Spanish cohort (262 cases and 920 controls) and 827 in the British cohort (639 cases and 189 controls). The three cohorts were first analyzed for association separately in PLINK v1.07. The same additional SNP filtering steps were applied: SNPs with (i) MAF < 0.5%, or (ii) genotyping rate < 95%, or (iii) deviation from Hardy-Weinberg equilibrium in controls with a P-value < 0.0001, or (iv) significantly different genotyping rate between case and controls were removed. Two SNPs were filtered out, so 25 SNPs were left for the association analysis. The same procedure was used to perform the association test with a logistic regression model and correcting for population heterogeneity. The GC reported were 1.478, 1.215 and 1.000 for the Western European, British and Spanish cohort respectively.
Validation by GWAS array typing
As some of the SNPs called during the discovery phase were included in the GWAS array used on the same cohorts previously15, they were not included in the Sequenom assay and the data from the GWAS array (typed and imputed genotypes) was used instead. The selected SNPs with an association P-value in the discovery phase < 0.05 were extracted from the GWAS data in the Western European (excluding the samples present in the discovery cohort), Spanish and British cohorts. The following number of samples were available: 2541 in Western European cohort (147 cases and 2394 controls), 1332 in the Spanish cohort (422 cases and 910 controls) and 5178 in the British cohort (475 cases and 4703 controls). The same SNP filtering in PLINK was done: SNPs with (i) MAF < 0.5%, or (ii) genotyping rate < 95%, or (iii) deviation from Hardy-Weinberg equilibrium in controls with a P-value < 0.0001, or (iv) significantly different genotyping rate between case and controls were removed. This resulted in 744 SNPs to analyze for association, the analysis was performed as mentioned above. The GC coefficients reported were 1.000, 1.100 and 1.113 for the Western Europe, UK and Spain cohorts respectively. Each validation cohort was analyzed separately and a meta-analysis together with the discovery cohort was performed. In order to further select the SNPs to follow up, we focused on those variants that were typed rather than imputed, present in all three validation cohorts, and with a discovery P-value < 0.01 which resulted in 9 SNPs remaining.
Meta-analysis
In order to combine and analyze the discovery together with the validation cohorts, a meta-analysis was performed using the PLINK function–meta-analysis and including the 4 cohorts (Discovery-Western Europe, Validation-Spain, Validation-UK and Validation-Western Europe). The value I given by the program was used as a measure of the heterogeneity between cohorts, thus when I was low (<5), the fixed-effects P-value and OR was used while when I was high (>5), the random effect P-value and OR was reported.
Evidences for significance
The following resources were used to further investigate the two SNPs reported in this study: db SNP: https://www.ncbi.nlm.nih.gov/projects/SNP/. Known alleles and frequencies for the SNPs in the European population were inquired in this data base using the SNP rs ID as search word. HaploReg (version 4.1): http://archive.broadinstitute.org/mammals/haploreg/haploreg.php. Search was conducted using the SNPs rs ID. Regulome DB (version 1.1): http://www.regulomedb.org/index. The list of SNPs showing a meta-analysis P-value < 0.05 was submitted to the web interface and a score for each SNP was returned as well as information described in the text. HiView: http://www.unc.edu/~yunmli/HiView/. This website allows browsing of the data on chromatin interaction in IMR90 fibroblast cell line from Jin et al.22. Search was conducted using a SNP rs ID as input, returning the fragment of DNA containing the SNP as bait, highlighted in blue, and the interactions observed with the surrounding regions (Supplementary Fig. S8). Data bases of known NF-kB targets: http://bioinfo.lifl.fr/NF-KB/ and https://www.bu.edu/nf-kb/gene-resources/target-genes/. The two lists were browsed for LIF, to examine whether it is a known NF-kB target. UCSC genome browser: http://genome.ucsc.edu/. Our own genomic data (RELA and H3K27ac ChIP-seq) were loaded to examine the signal at each SNP location. In addition, the following public tracks were added: UCSC genes, ENCODE transcription factor binding sites, layered H3K27ac mark and PhyloP score.
Allelic bias
In order to investigate allelic bias, the genomic data generated in Detroit 562 and FaDu cells were used. For each cell line, the reads from all RELA ChIP-seq experiments performed previously were pooled and the subset of them located 2 Kb around the SNPs identified in our study were extracted using samtools view and samtools merge. The same was done for H3K27ac ChIP-seq in Control and LPS condition as well as Input DNA. Subsequently, the newly extracted bam files were sorted (samtools sort), duplicates were removed (samtools rmdup) and they were indexed (samtools index) in order to be loaded onto the IGV software. Each SNP was then examined and considered for heterozygosity in both cell lines.
Rs11913168 was found to be heterozygous in Detroit 562 cells and the number of reads carrying the G and A alleles, mapping on each of the plus or minus DNA strand was extracted for each condition. Strand bias was investigated with a Fisher exact test using the 2 × 2 table calculator at http://vassarstats.net/index.html. The Fisher one-tailed P-value is reported for each data set in Supplementary Fig. S8B. Actual allelic bias was then evaluated with a Binomial test against the probability of each allele occurring in 50% of cases, using the binom.test command in R. P-values from this test are reported in Table 3.
Cell culture and treatment
FaDu cells were purchased from ATCC and cultured in RPMI medium (Gibco) supplemented with 10% fetal bovine serum performance (Gibco), 100 U/ml penicillin and 100 ug/ml streptomycin (Gibco), and 1 mM sodium pyruvate (Gibco). LPS from E. coli B4:111 (Sigma) at 1 μg/mL was used to stimulate the cells. Inhibition of NFKB was achieved using BAY 11–7082 (ChemCruz) at 50 uM (stock solution at 50 mM in DMSO) for one hour prior and at 45 uM (dilution 9/10 in the medium containing the ligand) during the treatment. RELA activation as well as gene expression were performed as described in Supplementary Table S7.
RELA activation assay
After treatment of the FaDu cells, the cytoplasmic and nuclear proteins were extracted using the NE-PER kit (Pierce) and quantified with the Coomasie Protein quantification kit (Pierce). Ten to 15 ug of nuclear proteins were used to test RELA activation using the NFκB p65 transcription factor assay kit (Pierce) according to the manufacturer’s instructions.
RT-qPCR
FaDu cells were treated with LPS or fresh medium for 80 minutes and total RNAs were extracted. One microgram of was used for reverse transcription using QuantiTect Reverse Transcription kit (Qiagen) followed by qPCR with LightCycler 480 SYBR Green I Master (Roche) according to the manufacturers’ instructions. The following PCR program was used: 5 minutes at 95 °C; 45 × (10 seconds at 95 °C, 1 minutes at 65 °C, 30 seconds at 72 °C). The ΔΔCt method was applied to determine the fold change of the genes of interest under LPS compared to the control, normalized with 2 housekeeping genes (ACTB, GAPDH). Sequences of the primers used for qPCR are shown below.
Use of human samples
All samples have been collected under the Ethikkomission of the Medizinische Universitat Graz (Reference Number 24–116 ex 11/12), methods were carried out in accordance with relevant guidelines and regulations and all experimental protocols were approved by the same Review board. Informed consent was obtained from all subjects or, if subjects are under 18, from a parent and/or legal guardian.
Ethical approval and informed consent
All samples have been collected under country-specific Institutional Review Boards. The following review boards approved the protocols: St Mary’s Research Ethics Committee (EC3263) and Medical Ethics Committee Erasmus MC (MEC 154.304/1996/139, MEC 190.703/2000/79 and MEC 217.867/2002/226).
A common clinical protocol agreed by the EUCLIDS Clinical Network and approved by all ethics committees was implemented at all hospitals. All clinical staff were trained in the project’s procedures and specified criteria were used for clinical definitions and assignment of patients to diagnostic categories. Written informed consent was obtained from a parent or legal guardian for each patient before inclusion in the study.
Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.
References
Edwards, S. L., Beesley, J., French, J. D. & Dunning, A. M. Beyond GWASs: illuminating the dark road from association to function. American journal of human genetics 93, 779–797, https://doi.org/10.1016/j.ajhg.2013.10.012 (2013).
Reddy, T. E. et al. Effects of sequence variation on differential allelic transcription factor occupancy and gene expression. Genome research 22, 860–869 (2012).
Kilpinen, H. et al. Coordinated effects of sequence variation on DNA binding, chromatin structure, and transcription. Science 342, 744–747 (2013).
Albert, F. W. & Kruglyak, L. The role of regulatory variation in complex traits and disease. Nature reviews. Genetics 16, 197 (2015).
Borghini, L., Lu, J., Hibberd, M. & Davila, S. Variation in Genome-Wide NF-kappaB RELA Binding Sites upon Microbial Stimuli and Identification of a Virus Response Profile. Journal of immunology (Baltimore, Md.: 1950), https://doi.org/10.4049/jimmunol.1800246 (2018).
Borghini, L., Hibberd, M. & Davila, S. Changes in H3K27ac following lipopolysaccharide stimulation of nasopharyngeal epithelial cells. BMC Genomics 19(1) (2018).
Stephens, D. S., Greenwood, B. & Brandtzaeg, P. Epidemic meningitis, meningococcaemia, and Neisseria meningitidis. Lancet 369, 2196–2210, https://doi.org/10.1016/S0140-6736(07)61016-2 (2007).
Caugant, D. A., Tzanakaki, G. & Kriz, P. Lessons from meningococcal carriage studies. FEMS microbiology reviews 31, 52–63 (2007).
Mogensen, T. H., Paludan, S. R., Kilian, M. & Ostergaard, L. Live Streptococcus pneumoniae, Haemophilus influenzae, and Neisseria meningitidis activate the inflammatory response through Toll-like receptors 2, 4, and 9 in species-specific patterns. J Leukoc Biol 80, 267–277 (2006).
Girardin, S. E. et al. Nod1 detects a unique muropeptide from gram-negative bacterial peptidoglycan. science 300, 1584–1587 (2003).
Hayden, M. S. & Ghosh, S. Signaling to NF-κB. Genes & development 18, 2195–2224 (2004).
Haralambous, E. et al. Sibling familial risk ratio of meningococcal disease in UK Caucasians. Epidemiology and infection 130, 413–418 (2003).
Wright, V., Hibberd, M. & Levin, M. Genetic polymorphisms in host response to meningococcal infection: the role of susceptibility and severity genes. Vaccine 27, B90–B102 (2009).
Parker, D. & Prince, A. Innate immunity in the respiratory epithelium. American journal of respiratory cell and molecular biology 45, 189–201, https://doi.org/10.1165/rcmb.2011-0011RT (2011).
Davila, S. et al. Genome-wide association study identifies variants in the CFH region associated with host susceptibility to meningococcal disease. Nature genetics 42, 772–776, https://doi.org/10.1038/ng.640 (2010).
Martinón-Torres, F. et al. Natural resistance to Meningococcal Disease related to CFH loci: Meta-analysis of genome-wide association studies. Scientific reports 6, 35842 (2016).
DePristo, M. A. et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nature genetics 43, 491–498 (2011).
Kundaje, A. et al. Integrative analysis of 111 reference human epigenomes. Nature 518, 317–330 (2015).
Aguet, F. et al. Genetic effects on gene expression across human tissues. Nature 550, 204–213 (2017).
Uhlen, M. et al. Towards a knowledge-based human protein atlas. Nature biotechnology 28, 1248–1250 (2010).
Calo, E. & Wysocka, J. Modification of enhancer chromatin: what, how, and why? Molecular cell 49, 825–837 (2013).
Jin, F. et al. A high-resolution map of the three-dimensional chromatin interactome in human cells. Nature 503, 290–294, https://doi.org/10.1038/nature12644 (2013).
Waragai, M. et al. Ataxin 10 induces neuritogenesis via interaction with G‐protein β2 subunit. Journal of neuroscience research 83, 1170–1178 (2006).
Tian, J. et al. Aurora B-dependent phosphorylation of Ataxin-10 promotes the interaction between Ataxin-10 and Plk1 in cytokinesis. Scientific reports 5 (2015).
Warner, N. et al. A genome-wide siRNA screen reveals positive and negative regulators of the NOD2 and NF-κB signaling pathways. Science signaling 6 (2013).
Warner, N., Burberry, A., Pliakas, M., McDonald, C. & Núñez, G. A genome-wide small interfering RNA (siRNA) screen reveals nuclear factor-κB (NF-κB)-independent regulators of NOD2-induced interleukin-8 (IL-8) secretion. Journal of Biological Chemistry 289, 28213–28224 (2014).
Attwooll, C., Denchi, E. L. & Helin, K. The E2F family: specific functions and overlapping interests. The EMBO journal 23, 4709–4716 (2004).
Lim, C. A. et al. Genome-wide mapping of RELA(p65) binding identifies E2F1 as a transcriptional activator recruited by NF-kappaB upon TLR4 activation. Molecular cell 27, 622–635, https://doi.org/10.1016/j.molcel.2007.06.038 (2007).
Ulich, T. R. et al. Intratracheal injection of LPS and cytokines. V. LPS induces expression of LIF and LIF inhibits acute inflammation. American Journal of Physiology-Lung Cellular and Molecular Physiology 267, L442–L446 (1994).
Nicola, N. A. & Babon, J. J. Leukemia inhibitory factor (LIF). Cytokine & growth factor reviews 26, 533–544 (2015).
Ferreira, V. P., Pangburn, M. K. & Cortés, C. Complement control protein factor H: the good, the bad, and the inadequate. Molecular immunology 47, 2187–2197 (2010).
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25, 1754–1760 (2009).
Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. The American Journal of Human Genetics 81, 559–575 (2007).
Acknowledgements
We thank Kumar Vikrant for his advices on the association analysis as well as Phua Zai Yang for his help on the library preparation and hybridization. The Next Generation Sequencing Platform and the Research Pipeline Development team at the Genome Institute of Singapore made the sequencing possible as well as data processing and mapping of the reads. In addition, we acknowledge A*STAR for the SINGA scholarship awarded to LB. Finally, we would like to thank all children and parents who participated in this study. S.D., E.P., M.H. and L.B. were supported by the Agency for Science Technology and Research of Singapore (A*STAR). This work has been partially funded by the European Seventh Framework Programme for Research and Technological Development (FP7) under EUCLIDS Grant Agreement No. 279185. The UK meningococcal cohort was established with support from Meningitis Research Foundation through grants to Imperial College London. The Research from Newcastle partners was supported by the National Institute for Health Research Newcastle Biomedical Research Centre based at Newcastle Hospitals NHS Foundation Trust and Newcastle University. The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health. The ESIGEM Research group activities were supported by grants from Instituto de Salud Carlos III (Proyecto de Investigación en Salud, Acción Estratégica en Salud: proyecto GePEM PI16/01478) (A.S.); Consellería de Sanidade, Xunta de Galicia (RHI07/2-intensificación actividad investigadora, PS09749 and 10PXIB918184PR), Instituto de Salud Carlos III (Intensificación de la actividad investigadora 2007–2016), Convenio de colaboración de investigación (Wyeth España-Fundación IDICHUS 2007–2011), Convenio de colaboración de investigación (Novartis España-Fundación IDICHUS 2010–2011), Fondo de Investigación Sanitaria (FIS; PI070069/PI1000540) del plan nacional de I + D + I and ‘fondos FEDER’ (F.M.T.). The Swiss Pediatric Sepsis study was funded by grants from the Swiss National Science Foundation (342730_153158/1), the Swiss Society of Intensive Care, the Bangerter Foundation, the Vinetum and Borer Foundation, and the Foundation for the Health of Children and Adolescents. The Western Europe Meningococcal Study was supported by grants no 8842, 10112 and 12710 of the Oesterreichische Nationalbank (Austria), grants A3-16.K-8/2008-11 and A3-16.K-8/2006–9 of the Department for Science and Research of the Styrian federal government (Austria) and the non for profit association ‘In Vita’, Graz (Austria).
Author information
Authors and Affiliations
Consortia
Contributions
S.D. designed the study, L.B. performed the experiments and the analyses with occasional help from E.P., M.H. provided support and comments of the results. A.B., V.W., R.dG., J.H., M.E., S.A., C.F., E.C., A.P., T.K., F.M.T., W.Z., M.L., M.H., S.D., E.P., M.V.dF., L.S., F.S., A.S. and L.C. obtained funding and designed the clinical study. All the EUCLIDS authors have contributed to patient recruitment and inclusion and clinical assessment. L.B. wrote the manuscript with the help of S.D.
Corresponding authors
Ethics declarations
Competing Interests
The authors declare no competing interests.
Additional information
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Borghini, L., Png, E., Binder, A. et al. Identification of regulatory variants associated with genetic susceptibility to meningococcal disease. Sci Rep 9, 6966 (2019). https://doi.org/10.1038/s41598-019-43292-6
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598-019-43292-6
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.
