Authors: Hyejeong Hong, Katharina V. Schulze, Ian E. Copeland, Manasa Atyam, Kendra Kamp, Neil A. Hanchard, John Belmont, Tamar Ringel-Kulka, Margaret Heitkemper, Robert J. Shulman
Categories: Article, Carbohydrate Digestive Enzymes, Common Variants, Inheritance, Irritable Bowel Syndrome, Rare Variants
Source: medRxiv
Authors: Hyejeong Hong, Katharina V. Schulze, Ian E. Copeland, Manasa Atyam, Kendra Kamp, Neil A. Hanchard, John Belmont, Tamar Ringel-Kulka, Margaret Heitkemper, Robert J. Shulman
Irritable Bowel Syndrome (IBS) is characterized by abdominal pain and alterations in bowel pattern, such as constipation (IBS-C), diarrhea (IBS-D), or mixed (IBS-M). Since malabsorption of ingested carbohydrates (CHO) can cause abdominal symptoms that closely mimic those of IBS, identifying genetic mutations in CHO digestive enzymes associated with IBS symptoms is critical to ascertain IBS pathophysiology. Through candidate gene association studies, we identify several common variants in TREH, SI, SLC5A1 and SLC2A5 that are associated with IBS symptoms. By investigating rare recessive Mendelian or oligogenic inheritance patterns, we identify case-exclusive rare deleterious variation in known disease genes (SI, LCT, ALDOB, and SLC5A1) as well as candidate disease genes (MGAM and SLC5A2), providing potential evidence of monogenic or oligogenic inheritance in a subset of IBS cases. Finally, our data highlight that moderate to severe IBS-associated gastrointestinal symptoms are often observed in IBS cases carrying one or more of deleterious rare variants.
Irritable Bowel Syndrome (IBS) is a common disorder affecting 10–15% of the global population; both adults and children.^1–5^ IBS is characterized clinically by abdominal pain associated with a change in bowel habits and categorized by bowel IBS-Constipation (IBS-C), IBS-Diarrhea (IBS-D), and IBS-Mixed bowel patterns (IBS-M).^3,6,7^ IBS is associated with diminished quality of life among both adults and children.^8,9^
Multiple studies have sought to determine the etiology of this complex disorder. In some cases there is evidence of altered immune function systemically as well as in the gut,^10–13^ altered gut barrier function,^13^ gut dysbiosis,^14,15^ and visceral hypersensitivity^16^ among other factors. Prior enteric infection is a risk factor for post-infectious IBS.^17,18^ Additionally, results from twin studies suggest a potential genetic and environmental influence to the expression of IBS.^19,20^
Several lines of evidence suggest that variations in digestive enzyme genes influence the risk of IBS. Genome-wide association studies (GWAS) have identified approximately 26 genetic variants associated with the risk of IBS in 25 genes inclusive of ASB4 and SPATA5.^21–24^ Although previous GWAS provide a valuable first insight into candidate genetic loci of IBS, the functional effect of most genetic variants from GWAS remains largely unexplained given the polygenic nature of IBS. In addition, studies have implicated rare recessive variation in genes that cause congenital deficiencies of carbohydrate (CHO) digestive enzymes,^25,26^ but the potential role of CHO digestive enzyme deficiencies in IBS is poorly defined. There is substantial clinical overlap between symptoms of IBS and the classical presentation of individuals with congenital enzymatic deficiency of these key digestive enzymes.^25,27^ Both monoallelic (homozygous recessive) and biallelic variation (compound heterozygous in trans) have been reported to cause congenital deficiency in sucrase-isomaltase (SI),^28^ maltase-glucoamylase (MGAM),^29^ lactase (LCT),^30^ trehalase (TREH),^31^ Solute Carrier Family 5 Member 1 (SLC5A1),^32^ Solute Carrier Family 2 Member 5 (SLC2A5),^33^ and hereditary fructose intolerance (Aldolase, Fructose-Bisphosphate B [ALDOB]). Given the heritability of IBS (0–57%)^28^ and the phenotypic parallels observed in individuals with congenital enzymatic deficiencies, we hypothesized that common and rare variation within the aforementioned digestive enzymes could contribute to IBS symptoms and diagnosis in a subset of patients.^34^ We undertook targeted deep-sequencing using next-generation technologies of these seven genes in a well-phenotyped IBS cohort inclusive of 687 adult and pediatric cases and 439 non-IBS controls from three regional sites in the United States. The generated data were then interrogated to identify IBS-associated common variants through candidate gene association studies as well as rare recessive Mendelian or oligogenic inheritance that would likely result in congenital deficiency.
We used genomic DNA that had been collected from saliva from 229 children with IBS and 130 healthy control (HC) children at Texas Children’s Hospital (TCH) and from blood from 458 adults with IBS and 309 HCs from the University of Washington (UW) and the University of North Carolina at Chapel Hill (UNC) (Figure 1). Our IBS sample cohort included a total of 1,126 participants (687 IBS cases and 439 HCs) with a median age of 27.0 (range, 7–83). The sample was predominantly female (74.5%, n= 839) and white (87.4%, n= 691) (Table 1). To assess whether genetic ancestry could confound our results, we conducted multidimensional scaling (MDS) using single nucleotide variant (SNV) data from the 1000 Genomes Project compared with common variation from individuals within the IBS cohort (see Methods). Aside from 60 individuals who did not report their race, our cohort included persons who self-identified as Native American or Alaskan Native (n= 3), Asian (n= 28), Black (n= 55), Hispanic or Latino (n= 103), White (n= 832), and Multiple or Other (n= 45). We found that self-reported ethnicities were highly consistent with genetic ancestry, showing 93.85% concordance (Figure 2). To account for the population stratification evident in our cohort, we included MDS components as covariates in our association tests. We did not, however, see a significant difference in genetic ancestry between cases and controls (p > 0.05, Fisher’s exact test), nor among IBS subtypes.
Association testing was conducted on IBS cases versus HCs for each candidate gene with an additive genetic model adjusted for sex and the first two MDS components. Five SNVs within TREH were associated with an IBS clinical diagnosis (Table 2). Of those, an increased number of minor alleles at rs45472704 and rs45529131 was positively associated with the IBS clinical phenotype with adjusted odds ratios (aORs) of 2.41 (95% confidence interval [CI]= 1.26–4.60) and 2.36 (95% CI= 1.23–4.52), respectively. The number of minor alleles at rs2277296 and rs2277297 were negatively associated with the IBS phenotype (aOR= 0.76, 95% CI= 0.62–0.93); these two SNVs were in strong linkage disequilibrium (LD, r^2^= 1), suggesting a single effect with dependent variants. Via GTEx portal, we further found that both rs2277296 and rs2277297 impact TREH expression levels in tissues of gastrointestinal (GI) tract. Each is an expression quantitative trait locus (eQTL) where the minor allele is associated with lower TREH expression in tissues, including the terminal ileum and transverse colon (Supplemental Figure 1). These two loci also showed strong evidence for regulation of transcription in the GI tract and are associated with regulatory motif alteration with Calcyclin-binding protein (CACYBP), according to RegulomeDB and HaploReg v4.1.^35,36^ In IBS subtype analyses, we found that the number of minor alleles at eight loci (rs17748, rs7928371, rs2276065, rs10790256, rs12225548, rs10892251, rs2277296, and rs2277297) in TREH are, cumulatively, significantly associated with decreased risk of IBS-C (aORs= 0.57–0.59, 95% CIs from [0.38–0.84] to [0.40–0.87]) (Table 2). Finally, the associations between genotypes and IBS-associated phenotypes were externally validated using the United Kingdom BioBank (UKBB) dataset (Supplemental Table 1a). From an independent panel from UKBB data, we also found that rs10892251, rs2277296, and rs2277297 were significantly associated with decreased sugar consumption behavior (i.e., less likely to have sugar added to tea) in the European population and metadata (β= −0.01, p< 0.01; Supplemental Table 1b).
We observed two variants in the SI gene (rs75172324 and rs62280366) in 124 cases and 439 controls. The number of minor alleles of each locus had negative associations with the odds of IBS-C subtype (aOR= 0.63, 95% CI= 0.46–0.86 and aOR= 0.64, 95% CI= 0.46–0.88, respectively), suggesting that those with these variants may have a lower likelihood of constipation but potentially higher odds of diarrhea. Likewise, these two variants are also significantly associated with the increased risk of diarrhea (β= 0.15, p <0.01; Supplemental Table 1c) in the UKBB European population.
SLC5A1 rs130406 was observed in 124 cases and 439 controls, and the number of minor alleles was associated with the reduced odds of IBS-C (aOR= 0.68, 95% CI= 0.51–0.91), adjusted for sex and first six MDS components. In addition, rs130406 was associated with IBS in the UKBB European-GWAS (β= 0.05, p= 0.005; Supplemental Table 1d).
SLC2A5 rs113665082 and rs111429581 were also found to be in strong LD (r^2^=1) and were identified in 108 cases and 439 controls from the subset analysis. Given that IBS-U is defined as an unspecified subtype due to a lack of available specification data, we observed strong evidence of an association with IBS-U risk, with an aOR of 6.8 (95% CI=2.09–21.94) after controlling for sex and the first six MDS components.
Rare deleterious variants were identified by assessing a deleterious scoring matrix combining damage and conservation prediction scores as well as potential monogenic or oligogenic inheritance. We further investigated the relationship between rare deleterious variants and the severity of IBS symptoms, including abdominal pain, diarrhea, constipation, and bloating based on Rome questionnaires rated on a 0–4 scale (Supplemental Table 2), subsequently rescaled to a 0–10 Likert-like scale. A single IBS case (RS16038_S63_L002) had biallelic ALDOB rare variants both predicted to affect splicing (NM_000035.4:c.799+6G>A and NM_000035.4:c.324+8C>G) with Transcript-inferred Pathogenicity (TRaP) scores of 0.134 and 0.249, respectively. Both have consensus benign variant classifications in ClinVar with respect to classical Hereditary Fructose Intolerance. This affected individual reported severe constipation with a score of 8 out of 10. No controls had ALDOB candidate rare variants. One IBS case (B-184_S147_L004) had biallelic rare missense variants in LCT (p.Lys88Asn and p.Arg1463Lys). Neither are present in the public databases used in this study and both are predicted to be damaging using the Combined Annotation-Dependent Depletion (CADD) scores (21.2 and 10.8, respectively), but are less clear with the rare exome variant ensemble learner (REVEL). This individual tended to have alternating diarrhea and constipation with severity scores of 8 and 6, respectively. No controls had similar candidate variants in LCT. We also observed one IBS case (RS2165_S355_L008) with biallelic variants in SLC5A1 (NM_000343.4:c.1666–5T>C, TRaP score 0.127, p.Pro269His, CADD score 24.3, and REVEL score 0.59). This individual reported severe diarrhea (10 out of 10). No rare biallelic deleterious SLC5A1 variants were observed in HCs. Rare and damaging biallelic missense variants were identified in four cases in SI, while no similar variants were observed among HCs (Table 3 and Supplemental Table 3). These four IBS cases with SI variants tended to present with moderate to severe abdominal pain, diarrhea and bloating (as rated on a scale of 6–10, Table 3).
Using damage prediction score and conservation prediction score matrices, we identified rare and damaging biallelic missense variants in three IBS cases, and two of these three cases (RS16176_S113_L003 and 2–143_S273_L006) were identical to the variants independently identified by CADD and REVEL scores. We further validated that three variants were well conserved and predicted to be damaging to the resulting protein in ClinVar: p.Val577Gly was classified as pathogenic, p.Val371Met as likely benign, and p.Tyr975His as having conflicting interpretations (Table 4 and Supplemental Table 4). In these three IBS cases with SI variants, moderate to severe levels of abdominal pain and bloating were commonly reported (scored from 6 to 10). In addition, we identified rare damaging biallelic variants in MGAM in two IBS cases, and one of these missense variants (p.Thr1628Met) was classified as benign in ClinVar (Table 4 and Supplemental Table 4) and both cases reported moderate constipation (6 out of 10). None of the HCs had MGAM candidate rare variants. In none of these cases were the variants detected by genetic testing sufficient to give a diagnosis of an inherited enzyme deficiency.
We also wished to address the possibility of a cumulative intolerance to dietary carbohydrate via digenic or oligogenic inheritance. Using stringent filters for deleterious variants, we identified five cases with a rare coding or splicing variant in at least two known enzyme deficiency genes (Table 5 and Supplemental Table 5). A complex pattern of IBS symptoms was reported in these five IBS cases, with an alternating diarrhea and constipation bowel pattern being the main symptoms in three cases. None of HCs had these combinations of variants. Considering SLC5A2 and MGAM as candidate disease genes, another nine cases had at least one variant from a known disease gene and these two candidate genes; while four HCs also had such combinations (Supplemental Table 5). Importantly, we identified that five IBS individuals had digenic or oligogenic recessive inheritance patterns in SI and MGAM; whereas only one HC did (Supplemental Table 5).
To assess the involvement of rare variant burden in IBS pathophysiology, we employed the kernel-based multiple regression model to identify rare variation associated with IBS case-control phenotypes; however, the burden of rare variants (minor allele frequency [MAF] ≤ 0.05) in ALDOB, TREH, LCT, MGAM, SI, SLC2A5, and SLC5A1 was not associated with IBS case-control phenotype (Table 6).
Although malabsorption, particularly of CHO, can cause symptoms (e.g., abdominal pain, bloating) that mimic those of IBS, there are few studies exploring the potential role of CHO digestive enzyme deficiencies in IBS. Here, we identify several common and rare variants in CHO digestive enzymes associated with IBS clinical phenotypes and symptoms in adults and children. Our results highlight (i) several common variants in TREH, SI, SLC5A1 and SLC2A5 are associated with IBS, and the IBS-C and IBS-U phenotypes; (ii) case-exclusive rare deleterious variation in known disease genes, such as SI, LCT, ALDOB, SLC5A1 as well as candidate disease genes such as MGAM and SLC5A2 provided potential evidence of monogenic or oligogenic inheritance in a subset of IBS cases; and (iii) moderate to severe IBS-associated GI symptoms were often observed in IBS cases carrying one or more of deleterious rare variants.
Traditionally, these enzyme deficiencies are considered rare, genetically inherited in a recessive fashion, and mainly the result of “loss-of-function” mutations.^37–40^ Although it is believed that congenital SI deficiencies cause symptoms recognized in infants or toddlers, more recent data suggest that the presentation may be more subtle in the case of the newly recognized allelic variants.^41,42^ Indeed, previous data from pediatric studies suggest that SI deficiency can be found in some children with recurrent abdominal pain, because enzyme deficiencies may be more common than previously recognized.^43,44^ In a previous IBS study in adults, a low CHO diet reduced abdominal pain and stooling frequency in some IBS-D patients, underscoring the potential role of mutations and variations in function producing IBS-like symptoms.^45^
SI exhibits a wide α-glucosidase activity profile and cooperates with MGAM in digesting α-1,4 linkages, the major glycosidic linkages in starchy foods.^46,47^ SI and MGAM are paralogues that both possess GH31 family catalytic domains and a Trefoil (P-type) domain.^48,49^ We identified the presence of rare SI deleterious variants with or without MGAM are more likely associated with an increased risk of IBS-associated GI symptoms, including abdominal pain, diarrhea, and bloating. In addition, two variants were previously reported as either pathogenic or benign in Clinvar. Although Clinvar classifications were not made with respect to an IBS phenotype, it is likely these variants contribute to a partial loss of function and impact digestion. Screening for patients with rare variants such as SI with or without MGAM may be beneficial for identifying subgroups of patients and providing personalized treatment; dietary reduction of sucrose and starch and/or enzyme supplementation may be particularly beneficial for such individuals.^38,50^
TREH converts trehalose to glucose.^51^ In one study, trehalase deficiency, an autosomal recessive trait, was associated with abdominal pain, distension, flatulence, vomiting and diarrhea after trehalose ingestion.^52^ Our study found that several common variants in TREH were associated with an IBS diagnosis. Importantly, the odds of having IBS-C is 42% less likely according to the number of minor alleles at either rs2277296 or rs2277297. Such associations are perhaps mediated by, in part, downregulation of TREH expression because rs2277296 and rs2277297 are eQTLs in diverse tissue types in the GI tract, such as esophageal mucosa, terminal ileum, colon, pancreas, and liver, as shown by GTEx. For example, individuals carrying two copies of minor alleles (C/C) at rs2277296 may have lower TREH expression in the small intestine and reduced trehalase activity, leading to an increased risk of having an IBS-D phenotype. Moreover, a subset of common variants in TREH (rs10892251, rs2277296, and rs2277297) may be evidence of dietary selection pressures within certain populations that have traditionally increased consumption of sugar-sweetened beverages. This supposition is explained partly by the UKBB data because minor alleles of these three variants were also associated with less sugar consumption.
We also examined genetic variations in SLC5A1 and SLC2A5 which are involved in glucose, galactose, and fructose intestinal uptake.^53–55^ Cellular uptake of glucose is driven by two glucose transport sodium driven glucose cotransporters (SGLT) and the facilitative glucose transporters (GLUT) encoded by the solute carrier genes SLC5A1 and SLC5A2.^56,57^ Loss of SGLT results in glucose-galactose malabsorption, a rare congenital autosomal recessive disorder characterized by severe diarrhea starting at birth.^53,55^ Our common variant analysis detected that having at least one copy of a minor allele at rs130406 in SLC5A1 was associated with 32% decreased risk of having an IBS-C phenotype. In addition, one IBS individual with digenic inheritance in SI and SLC5A1 reported moderate diarrhea, while other IBS symptoms were reported as mild. Together, these data suggest that variants in SLC5A1 increase the risk of an IBS-D phenotype.
Although there is conflicting clinical and experimental evidence on a pathophysiological role for the fructose transporter GLUT5, most studies highlight a role of SLC2A5 gene as a major intestinal transporter for fructose and signaling molecule for the induction of down-stream genes encoding fructolytic and gluconeogenic enzymes.^58–62^ Studies from animal and cell line models support that fructose may positively regulate SLC2A5/GLUT5 expression via transcriptional and posttranscriptional modification.^63–66^ However, ex vivo studies of human intestinal biopsies did not show significant genetic effects of SLC2A5 on fructose intolerance.^67^ To understand the role of fructose intolerance in the context of IBS, we tested our hypothesis that genetic variation in ALDOB in conjunction with SLC2A5 synergizes the association with the risk of fructose intolerance with GI symptoms, such as abdominal pain and diarrhea. Pathogenic mutations that strongly affect the function of ALDOB, a key enzyme for metabolizing fructose, are clinically evident in very few people.^68^ First, we identified that individuals with two common variants, rs113665082 and rs111429581 in SLC2A5, have an approximately seven times higher odds of having IBS-U. Additionally, we identified that one HC had digenic recessive inheritance patterns in ALDOB and SLC2A5; whereas no IBS case individuals did. Considering fructose intolerance has monogenic causes, our data shows that only one IBS case has rare recessive variant in ALDOB with moderate to severe constipation; whereas no HCs did. Together, due to its rarity, our data is limited to determine the clinical associations between the risk of IBS and monogenic or digenic inheritance in ALDOB with or without SLC2A5. Further studies are required to establish a role for these genetic variants in ALDOB and SLC2A5 in heightening the risk of IBS by modulating fructose metabolism; if validated, they may provide potential targets for therapeutic dietary intervention.
Lactase persistence - the ability to digest milk lactose during adulthood – has a genetic origin; the frequency of lactase persistence is high in Northern European populations and decreases across East Asia and Southern Africa.^69^ In the present study, only two rare missense variants were identified in one IBS individual; while no HC contain these variants. Due to the dominance of European populations in our study, our ability to identify lactase-associated variants may be limited.
Our study has several limitations. First, we did not include all potential predispositions to an IBS phenotype - for example, diet, previous enteric infection, and/or stress.^70^ Rather, we focused on genetic variants associated with GI symptoms in IBS to increase understanding how CHO digestive capacity contributes to IBS manifestations. Second, we acknowledge that rare variants are commonly defined as a genetic variant with a MAF of < 0.01, which was the definition applied to our primary analysis. Since monogenic variants of common complex disorders like IBS can be very frequent, we applied a MAF of ≤ 0.05 in our subsequent monogenic and digenic variant association analyses. Finally, we conducted targeted genomic sequencing that focused on a panel of CHO digestive enzyme genes instead of using whole genome profiling, such as whole-genome or whole exome sequencing. Although these latter methods may create a greater opportunity to find pathogenic variants, targeted sequencing is an invaluable tool for mutation detection associated with greater sequencing depth allowing us to test our hypothesis that CHO malabsorption-related gene variants may play a role in IBS symptom generation. Furthermore, since targeted sequencing can screen a large number of samples with fast turnaround and reduced costs, it is clinically relevant in practice.
Our study has many strengths. First, we identified novel common and rare variants in CHO digestive enzyme genes and their association with IBS phenotypes. Such findings were previously unknown and have clinical significance, underlining the potential clinical value of genotyping. We anticipate that sequencing or genotyping efforts in individuals that have a family history of deficiencies in CHO digestive enzymes and/or early onset of IBS may facilitate precision medicine initiatives and impact treatment regimen. For example, our findings support evidence of close evolutionary and functional relationships in SI and MGAM, suggesting a role for disaccharide gene aggregation in IBS pathophysiology. For people with oligogenic inheritance with rare variants in SI and MGAM, avoiding specific dietary CHO may be beneficial. Second, due to the phenotypic complexity of IBS, we formed a consortium of clinical centers to diversify our study population, which improves statistical power. Data harmonization in the consortium enabled collection of greater racial and ethnic representation from three regional sites in the United States with a large sample size and identification of IBS subtype-defining SNVs. To avoid potential sources of bias, we used the same sequencing platform. Further, the IBS population was well-phenotyped using Rome II and III screening questionnaires. This helped maintain a consistent and accurate definition of the IBS group and minimized misclassification. Finally, the associations of common variants with IBS clinical phenotypes were externally validated using public databases, including UKBB, the largest population-based cohort in the United Kingdom, and GTEx portal, a large-scale consortium collecting tissue-specific gene expression data. Such cross-validated findings with a spectrum of phenotypes strengthen the clinical impact of candidate genes and variants on the increased risk of IBS. There is a need for future research to delve into the impact of these variants on CHO enzymatic function, the interplay between monogenic and polygenic factors in IBS risk, and the potential advantages of clinical trials for personalized dietary interventions to optimize patient care and therapeutic strategies.
This study used banked genomic DNA isolated from saliva and blood samples from well characterized cohorts of adults and children with IBS and HCs from three collaboration TCH, UW, and UNC (Figure 1). Institutional Review Board approval from each institute that included informed consent from participants (and assent from children) was obtained. For the IBS cohort, we included all participants who met the following (i) children and adults with well-phenotyped IBS (IBS-C, IBS-D, and IBS-M) based on Rome II (adults) and Rome III (children) screening criteria^7^. We excluded patients who bowel surgery; documented GI disorders (e.g., Crohn’s disease); a serious chronic medical condition (e.g., diabetes); chronic conditions with GI symptoms (e.g., cystic fibrosis) or celiac disease. The HC cohorts had no GI complaints. Details, including IBS study case definitions, are published for the TCH cohort ^15,71^, UW cohort,^72,73^ and UNC cohort^74^. A total of 1,126 samples were sequenced at the Laboratory for Translational Genomics at the USDA/ARS Children’s Nutrition Research Center at Baylor College of Medicine.
Genomic DNA was sheared to 300–400 base-pair (bp) with a sonicator (Covaris, Fischer, San Diego) and the fragments purified. The sheared DNA from each sample was used to generate a paired-end library (Illumina TruSeq DNA Sample Prep Kit). After DNA quantitation and quality control (Biorad Experion), 500 ng of library products were hybridized to the custom capture bait material. A custom Nimblegen SeqCap EZ Choice library (Roche) was constructed to capture the full candidate genes, including 5’ untranslated region (UTR), all exons, introns, and 3’UTR. Seven candidate IBS genes were SI, MGAM, LCT, TREH, SLC5A1, SLC2A5 and ALDOB. Samples were multiplexed (48 samples per lane) to facilitate high-depth coverage (mean coverage of >450X) of a total target length of 667 kilobase-pair (kb) on an Illumina Hiseq 2000 using 100 bp paired end reads.
Paired-read FASTQ files were aligned to the hg19 human reference genome using BWA-MEM (version 0.7.12)^29^ with shorter split hits marked as secondary (-M parameter). The functions SortSam, MarkDuplicates, and BuildBamIndex from the Picard toolkit (version 1.84)^30^ were used to sort aligned reads by coordinates, mark duplicate reads, and create index files for each alignment file. Realignment of reads was performed using the RealignerTargetCreator and IndelRealigner functions from the Genome Analysis Toolkit (GATK, version 4.0.5.1)^31^. To recalibrate base calls, GATK’s functions of BaseRecalibrator, AnalyzeCovariates, and PrintReads were used.
After removal of unmapped reads with samtools (v1.2)^32^ and subsequent indexing of cleaned alignment files with samtools’ index function, variants were called with GATK’s HaplotypeCaller. Following the creation of chromosome specific database files with GATK’s function GenomicsDBImport, joint genotyping was performed across all sequenced samples using the GenotypeGVCFs function from the same toolkit. Resulting chromosome specific joint variant call format (VCF) files were merged with GATK’s GatherVcfs function. GATK version 4.0.4.0 was used for all joint calling related functions (GenomicsDBImport, GenotypeGVCFs, and GatherVcfs).
Alongside variants within targeted gene regions the resulting joint-called VCF file also contained off-target variant calls, albeit at low coverage. Subsequent association tests were only run on targeted gene regions, the boundaries of which were determined by the first and last instance of 10X median coverage surrounding genes of interest across all sequenced samples. Per base coverage for each sample was calculated using bedtools’ coverage function (v2.26.0)^33^ on aligned and fully processed sequencing reads. All functions listed were run with default parameters unless otherwise stated.
In PLINK (v1.9)^75^, a total of 9,161 SNVs were called on seven targeted genes and filtered by Hardy-Weinberg Equilibrium (p < 1e^−6^), minor allele frequency (MAF ≥ 0.01), call rate (> 0.95), and linkage disequilibrium (LD, r^2^ < 0.2 in 50 bp windows with a 5 bp slide). Multiallelic loci were removed before merging IBS study samples and those from the 1000 Genomes Project Phase 3 based on the remaining variants’ dbSNP rsIDs (build 151), which were added with SnpSift’s annotate function (version 4.3t).^35^ Identity-by-state was calculated with PLINK cluster constraint (CC) parameters on the filtered 1,449 SNVs. Using MDS, we found no evidence of batch effects or any global influence of age or the other available demographic data from each cohort when plotting the first three dimensions. We further used the first three MDS components to calculate the Euclidean distance between every IBS sample and each group of individuals from the five 1000 Genomes superpopulations (‘AFR’: African, ‘AMR’: Admixed American, ‘EAS’: East Asian, ‘EUR’: European, and ‘SAS’: South Asian); an individual’s ancestral background was assigned based on the nearest superpopulation. To assess concordance, self-reported ethnicity labels were mapped to superpopulations as ‘Native American / Alaska Native’ = EAS, ‘Asian’ = EAS or SAS, ‘Black’ = AFR, ‘Hispanic / Latino’ = AMR, ‘Native Hawaiian / Pacific Islander’ = EAS, and ‘White’ = EUR. Individuals who did not report their ethnicity and those who identified as “Multiple / Other” were omitted from concordance calculations.
Variants in each of the seven genes were identified among 1,126 samples comprising of 687 IBS case and 439 HC samples. For each of the identified and filtered 1,449 SNVs, associations between genotype and clinical phenotype (IBS cases and HCs) were explored using an additive genetic regression model accounting for sex and the first two MDS components using PLINK (v1.9). Significant variants (MAF ≥ 0.01, p <0.01) were then reassessed across association tests using up to six MDS components. A subtype analysis was also performed for each of four IBS categories by comparing the allele frequencies between individuals in the subtype (IBS-D n= 206, IBS-C n=124, IBS-M n=47, IBS-U n=186) and the control samples, i.e., excluding cases in the other subtypes. QQ and Manhattan plots were generated using qqman in R (v1.0.136) for each analysis to determine deviation from expected. Results were filtered for SNVs that met a significance threshold of p < 0.01. We further queried the Broad Institute GTEx portal (https://gtexportal.org) to access tissue-specific gene expression data and identify SNVs that impact expression levels of these seven genes in GI tissues.
The Pan-Ancestry GWAS of UKBB (https://pan.ukbb.broadinstitute.org/docs/technical-overview) is a publicly available resource of imputed common variant association scores that provides a multi-ancestry analysis of 7,221 phenotypes using a generalized mixed model association testing framework, spanning 16,119 GWASs performed using the UKBB cohort. Because there is no specific IBS phenotype that is defined in the same way as our primary research cohort, we reviewed the available UKBB phenotypes for phenotype overlap and for possible relevant endophenotypes (https://pan.ukbb.broadinstitute.org/downloads). We selected clinical phenotypes as listed in Supplemental Table 1a.
The frequency of each candidate variant was compared with MAF from the Genome Aggregation Database (gnomAD).^76^ A SNV was deemed to be rare if it had a MAF < 0.01 in gnomAD. An allelic-ratio threshold of 30 was used to as a quality control step to reduce false positive variation. We focused on non-synonymous variants.
VarCards was used to functionally annotate the in situ likelihood of variants being protein-damaging, and included annotations from Combined Annotation-Dependent Depletion (CADD),^77^ Sorting Intolerant From Tolerant (SIFT),^78^ likelihood ratio test (LRT),^79^ Polymorphism Phenotyping v2 (Polyphen2),^80^ Functional Analysis through Hidden Markov Models (FATHMM)^81^ and Protein Variation Effect Analyzer (Provean).^82^ ‘Damaging’ variants were inferred from the following annotation “damaging”, “possibly damaging”, “probably damaging”, “disease causing”, and CADD scores > 12. Each variant was then given a final ‘damage prediction score’, representing the proportion of damage prediction algorithms that supported a ‘damaging’ designation for the variant in question; specifically a variant was considered damaging if it had a damage prediction score higher than three of six assessment criteria. Evolutionary conservation of non-synonymous sites, annotated by VarCards, was used for orthogonal assessment of ‘damaging’, using Genomic Evolutionary Rate Profiling (GERP),^83^ phylogenetic P-values (phyloP),^84^ SIte-specific PHYlogenetic analysis (SiPhy)^85^ and PHylogenetic Analysis with Space/Time models (phastcons).^86^ ‘Conserved’ variants were inferred from annotations as >12 for SiPhy, >4.4 for GERP, >1.6 for phyloP, and > 0.5 for phastcons. Each variant was then given a final ‘conservation prediction score’ representing the proportion of conservation algorithms supporting strong conservation. A variant was classified as ‘conserved’ if it had a conserved prediction in at least two out of four assessment criteria. In addition, we also compared variant annotations using seqR^87^ as confirmation.
In order for biallelic variants to be considered as disease-causing candidates, variants were required to be seen together only in cases (i.e., not biallelic in control samples) and have been characterized as pathogenic and evolutionarily conserved. This was viewed as a conservative measure to reduce false positives given the lack of parental data to determine whether variants were in the cis or trans configuration. Candidate variants were also queried in Clinvar^88^–a large database of curated putative disease-causing variation.
VCF files for each gene (n= 7) and individual (n= 1,126, IBS n= 687, HC n= 439) were merged. The merged file was uploaded into an anVIL workspace (https://anvil.terra.bio/) and variants were annotated for filtering in seqr (https://seqr.broadinstitute.org/). Focusing on variants included in transcripts and predicted splice sites and filtering for variants with MAF ≤ 0.05, resulted in a list of 2,314 variants in 803 individuals. Cases were further prioritized by the occurrence of ≥ 2 non-synonymous variants. Finally, to increase potential interpretability of cases, we focused on variants which were more likely to be deleterious by filtering for CADD scores ≥ 10 and rare exome variant ensemble learner (REVEL)^89^ scores ≥ 0.3.
In exploratory work we considered the effect of known protein motifs in our seven candidate genes in the calibration of CADD and REVEL deleteriousness score filtering. Each protein sequence was annotated with its corresponding pfam domains. The resulting list of domains was then used to select sequences in all additional human proteins annotated with those pfam domains resulting in a set of selected pfam intervals corresponding to all those domains. CADD and REVEL scores for all possible variants and all variants reported in gnomAD exome database were collected. Protein domain annotations and deleteriousness scores were obtained in the UCSC Table Browser (https://genome.ucsc.edu/cgi-bin/hgTables) in the GRCh37/hg19 reference. Distribution of scores was analyzed in R.^90^
VCF files corresponding to each gene (n= 7) including all individuals (n= 1,126) were parsed in PLINK1.9 (https://zzz.bwh.harvard.edu/plink) to produce binary input and site information files for association testing. Rare variant kernel-regression-based association tests were performed the R package SKAT.^91^ Three marker sets were (i) the full marker set consisting of both common and rare variants (n= 9,161 variant positions); (ii) rare variants only (n= 7,627); and (iii) rare variants with high deleteriousness scores (n= 260).
We investigated the relationship between rare genetic variants and severity/frequency of IBS symptoms, including abdominal pain, diarrhea, constipation, and bloating. IBS phenotype data were collected based on a Rome II screening questionnaire in adults and Rome III screening questionnaire in children, with severity ratings ranging from 0 to 4 (Supplemental Table 2). Subsequently, IBS symptom severity scores were converted on a 0–10 likert-like scale with 0 being no symptom and 10 the worst symptom that could be imagined. Because the statistical power of classical single-variant-based association tests for low-frequency and rare variants is extremely low, we report absolute scores that reflect the severity/frequency of IBS symptoms identified by individuals found to have rare variants, instead of testing the effects of candidate rare variants on IBS symptom severity.