Multiple apical plasma membrane constituents are associated with susceptibility to meconium ileus in individuals with cystic fibrosis.

Variants associated with meconium ileus in cystic fibrosis were identified in 3,763 affected individuals by genome-wide association study (GWAS). Five SNPs at two loci near SLC6A14 at Xq23-24 (minimum P = 1.28 × 10(-12) at rs3788766) and SLC26A9 at 1q32.1 (minimum P = 9.88 × 10(-9) at rs4077468) accounted for ~5% of phenotypic variability and were replicated in an independent sample of affected individuals (n = 2,372; P = 0.001 and 0.0001, respectively). By incorporating the knowledge that disease-causing mutations in CFTR alter electrolyte and fluid flux across surface epithelium into a hypothesis-driven GWAS (GWAS-HD), we identified associations with the same SNPs in SLC6A14 and SLC26A9 and established evidence for the involvement of SNPs in a third solute carrier gene, SLC9A3. In addition, GWAS-HD provided evidence of association between meconium ileus and multiple genes encoding constituents of the apical plasma membrane where CFTR resides (P = 0.0002; testing of 155 apical membrane genes jointly and in replication, P = 0.022). These findings suggest that modulating activities of apical membrane constituents could complement current therapeutic paradigms for cystic fibrosis.

Cystic fibrosis is caused by mutations in the cystic fibrosis transmembrane conductance regulator (CFTR) gene[1]. CFTR is a chloride channel located on the apical membrane of epithelia, where ion conduction and solute trafficking contribute to the regulation of transepithelial fluid flow. Individuals with the same loss-of-function CFTR mutations have variable disease severity, and differentially affected CF-associated organs including lung, pancreas, liver, intestines, and vas deferens; thus additional features, including other genes (referred to as modifier genes) may affect disease pathophysiology. Approximately 15% of CF patients have severe intestinal obstruction at birth, a complication known as meconium ileus[2]. Meconium ileus develops in utero, and presents following birth with complete intestinal obstruction that requires either medical or surgical intervention. This neonatal complication is highly indicative of CF, occurs in either sex, displays notable heritability exceeding 88%[3], and is likely minimally affected by environmental influences.

The North American CF Gene Modifier Consortium has accumulated 3,763 participants with ‘severe’ (pancreatic exocrine insufficient) CFTR alleles and genome-wide genotype data at 543,927 SNP loci[4,5] (Table 1 and Online Methods). The definition of meconium ileus was consistent within the consortium and was recorded following rigorous chart review. A conventional GWAS for meconium ileus used a generalized estimating equations (GEE) model[6] to include collected sibling-pairs, and led to five genome-wide significant SNPs (P<5×10−8)[7] from two regions that include SLC26A9 on chromosome 1 and SLC6A14 on chromosome X (Fig. 1, Supplementary Fig. 1, Table 2; sex-specific results in Supplementary Table 1). CFTR was not a significant confounder or effect modifier when incorporated in the GWAS (Supplementary Fig. 2 and Supplementary Table 2), indicating SLC6A14 and SLC26A9 are independent contributors to meconium ileus. We then replicated the associations in SLC6A14 (min P=0.001) and SLC26A9 (min P=0.0001) with meconium ileus in an independent combined collection from North America and France (Table 2).

Table 1

CF consortium participants of the meconium ileus study

Consortium Center/Site	Sample Sizea	CFTR Genotype		Status of CF Patientsb
	n	Phe508del/Phe508deln (%)	Phe508del/Otheror Other/Othern (%)	MI(M:Fc)n	Non-MI(M:Fc)n
CGS(Canadian Gene Modifier StudyPopulation Samples)	1661	992 (59.7)	669 (40.3)	252(128:124)	1409(772:637)
GMS-Lung Study(UNC/CWRU, Extreme PhenotypeStudy)	1120	1116 (99.6)	4 (0.4)	177(99:78)	943(498:445)
GMS-Liver Study(Samples from North America Only)	80	58 (72.5)	22 (27.5)	23(14:9)	57(35:22)
TSS(Johns Hopkins University,Family-based Twin and SiblingStudy)	902	522 (57.9)	380 (42.1)	159(78:81)	743(397:346)
Total	3763	2688 (71.4)	1075 (28.6)	611(319:292)	3152(1702:1450)

Size of CF patient sample used for the meconium ileus GWAS and GWAS-HD discovery study is given by consortium site. The TSS site contained sibling pairs. There were 3,199 unrelated individuals among the total of 3,763 studied.

Meconium ileus (MI) status is given separately for each center.

Gender break-down is indicated for each patient category.

Figure 1

Meconium Ileus GWAS identifies genome-wide significant SNPs. Association analysis was performed on all SNPs with minor allele frequencies > 2% that passed QC criteria (Online Methods).

(a) Genome-wide Manhattan plot for meconium ileus. The black solid line corresponds to the genome-wide significance threshold[7] (P<5×10−8), and the black dashed line to the suggestive association threshold, expected once per genome scan (P <1/543,927=1.84×10−6). A total of five SNPs in two regions (SLC6A14 on chromosome X, and SLC26A9 on chromosome 1) have association evidence exceeding the genome-wide threshold. The SNPs, rs4077468 and rs4077469 are in perfect LD and appear as one SNP as they are separated by only 128bp.

(b) Regional plot for SLC26A9. LocusZoom viewer was used to display the association evidence around SLC26A9 based on NCBI Build 36/hg18. Symbol coloring reflects HapMap CEU LD r2 values with the most significant SNP. The significant SNPs, rs4077468, rs7419153 and rs12047830 (Table 2), occur 2.17 kb, 4.72 kb and 4.12 kb upstream, respectively, of the transcription start site. The significant SNP, rs7512462, occurs in intron 5. A gap occurs in the genomic sequence between the SLC26A9 and FAM72A genes in both NCBI36.3 and GRCh37 primary reference assemblies.

(c) Regional plot for SLC6A14. The association evidence around SLC6A14 was displayed as above. Rs3788766 (Table 2) is positioned 0.95 kb upstream of the transcription start site and is within the binding site of the CEBPB transcription factor as annotated by ENCODE[33] (not shown). The mRNA transcript corresponding to CXorf61 is downstream (3’) of SLC6A14.

Table 2

Discovery and replication of SNPs in SLC26A9 and SLC6A14.

					Discovery Sample (611:3152)d				Replication Samplea
									N. America(251:889)d		France(154:1078)d		Combinedb(405:1967)d
SNP	CHR	POS	GENE	RiskAllele	Risk AlleleFrequency		OR	P	OR	P	OR	P	OR	P
					MI	Non-MI
rs4077468c	1	204181380	SLC26A9	T	0.66	0.57	1.45	9.88×10−9	1.45	0.0005	1.27	0.0575	1.37	0.0001
rs7512462	1	204166218	SLC26A9	T	0.66	0.57	1.45	2.14×10−8	1.31	0.0134	1.20	0.2120	1.27	0.0063
rs7419153	1	204183932	SLC26A9	T	0.44	0.36	1.42	1.01×10−7	1.42	0.0007	1.20	0.1290	1.33	0.0004
rs12047830	1	204183322	SLC26A9	C	0.56	0.49	1.34	3.72×10−6	1.33	0.0054	1.27	0.0510	1.31	0.0007
rs3788766	X	115480867	SLC6A14	T	0.72	0.59	1.50	1.28×10−12	1.14	0.1328	1.47	0.0006	1.25	0.0011
rs5905283	X	115479909	SLC6A14	C	0.61	0.50	1.34	1.69×10−8	1.04	0.6030	1.47	0.0002	1.19	0.0074
rs12839137	X	115479578	SLC6A14	C	0.82	0.75	1.39	1.20×10−6	1.20	0.0666	1.15	0.3275	1.18	0.0386

The North American replication collection consisted of 1,140 individuals that correspond to the continuing enrolment at the North American sites. The French cohort consisted of 1,232 patients collected at 38 French CF centers (Online Methods). The French replication samples were genotyped genome-wide with the Illumina 660W-Quad BeadChip platform while the North American replication samples were genotyped on the six SNPs using Taqman® Assay-on-Demand Assays.

The combined replication P value was calculated using the inverse-variance method, most powerful when different studies have the “same direction” of effect. The corrected type 1 error for the replication study is 0.008, using the conservative Bonferroni correction for the 6 SNPs tested.

There is another significant SNP on chromosome 1 (rs4077469, P=9.88 × 10−9) but not included as it is separated by only 128bp and in perfect LD with rs4077468.

Meconium ileus status (MI:Non-MI) break-down for each site.

The signal intensity plots of the associated SNPs reflected autosomal- and X-associated SNPs at SLC26A9 and SLC6A14, respectively. Imputation analysis using MACH and minimac[8,9] identified the same regions of association as the genotyped SNPs (Online Methods, Supplementary Fig. 3). The five associated SNPs in SLC6A14 and two in SLC26A9 (Fig. 1b and 1c) are positioned just upstream of their respective transcription start sites such that binding of activating or repressing transcription factors may be affected as highlighted by ENCODE data[10] (data not shown). Neither SLC6A14 nor SLC26A9 coding regions exhibit evidence for CNVs; however, there is a gap in the reference sequence >10 kb upstream of the SLC26A9 locus.

The seven SNPs genotyped (Table 2) in the two genes account for <5% of the meconium ileus variation, estimated by pseudo R-squared[11], likely reflecting the common problems in association studies of locus heterogeneity and low power given the available sample size. Whereas conventional GWAS is often considered for complex disease mapping, modifier gene studies could incorporate disease etiology and pathobiology information to increase power and account for heterogeneity. To do so, we proposed application of GWAS with consideration of a hypothesis (and so is hypothesis-driven; GWAS-HD) to systematically prioritize SNPs for genome-wide analysis. The highest priority markers are also evaluated as a set to test the statistical significance of the hypothesis used for prioritization (Supplementary Fig. 4).

The GWAS-HD prioritization in this CF application is based on the knowledge that a major source of CF pathophysiology is impaired fluid and electrolyte flux in epithelia of CF-affected organs. The polarized epithelial layer forms a highly selective barrier between organ and ductal interfaces. Transepithelial ‘function’ is achieved by cell polarization whereby many determinants and regulators of fluid, solute and ion transport reside at the apical membrane alongside CFTR, with contributing features from basolateral surfaces. We have shown in a mouse model that CFTR function in the gastrointestinal epithelium is critical for preventing intestinal obstructions[12]. Thus, we hypothesized that with loss of CFTR, (genetic) variation in other apical membrane constituents could modify CF phenotypes, such as meconium ileus.

A list of 157 gene products (Fig. 2 and Supplementary Table 3) was annotated as localized to the apical plasma membrane using AmiGO[13] with Gene Ontology data[14,15]. CFTR and many solute transporters were included. However, the brush border membrane protein SLC6A14 was not listed, reflecting the high specialization of its corresponding intestinal cavity and a limitation of the GO annotation that we accepted without additional curation to avoid bias. In total, 3,814 GWAS SNPs were within ±10 kb of the boundaries of 155 genes (NCBI36/hg18); 2 genes were not tagged by any of the GWAS SNPs.

Figure 2

The apical membrane hypothesis identifies genes associated with meconium ileus. A list of 157 genes was annotated using the AmiGO tool[13] (version 1.7; March 28, 2010) based on the Gene Ontology data[14] (GO:00163245) using the cell location search phrase “apical plasma membrane” with restriction to Homo sapiens. In total, 3,814 GWAS SNPs are within ±10 kb of the boundaries of 155 genes (NCBI36/hg18). Two genes were not tagged by any of the GWAS SNPs; SLC6A14 was not annotated to the apical plasma membrane.

(a) QQ-plot of the apical SNPs in the discovery sample. The observed association statistics (red), and the statistics calculated from the 10,000 phenotype-permutated replicates are shown (light gray).

(b) Statistical significance of the apical membrane hypothesis in the discovery sample. Statistical significance (permutation P = 0.0002) was established via a sum statistic, summing the association evidence (Wald χ21 statistic) over all the 3,814 SNPs with the observed sum statistic displayed as a vertical line (red), and the 10,000 permutation-based sum statistics displayed as a histogram (light gray).

(c) QQ-plot of the apical SNPs in the replication sample. The SNPs in the French replication cohort were required to have MAF > 6% because of the reduced sample size (1232 × 6% ≈ 3763 × 2%). In total 3,420 GWAS SNPs are within ±10 kb of the boundaries of 154 apical genes.

(d) Statistical significance of the apical membrane hypothesis in the French replication sample (permutation P=22/1,000=0.022).

To implement the GWAS-HD for meconium ileus using the apical hypothesis, we first prioritized the genome-wide markers by assigning the 3,814 SNPs of the apical genes to a high priority group and all remaining genome-wide SNPs to a low priority group. We then performed two statistical analyses (Supplementary Fig. 4 and Online Methods). The first, analogous to a conventional GWAS, was to conduct single-SNP association analysis using all of the 543,927 GWAS SNPs at the genome-wide level, however after up-weighting the 3,814 apical SNPs via the SFDR control.[16] The second analysis, focusing only on the 3,814 high priority SNPs using a multi-SNP/gene analysis, tested the prioritization hypothesis itself to determine if multiple proteins present on the apical plasma membrane contribute to meconium ileus susceptibility.

As in the conventional GWAS, SNPs from SLC6A14 showed the strongest evidence for association with meconium ileus in the single-SNP GWAS-HD analysis (Supplementary Fig. 5), despite not being in the high priority group, reflecting the robustness of SFDR[17]. In addition, SNPs from SLC26A9, and two additional apical genes, ATP2B2 and SLC9A3, showed association evidence with q values <0.05 (Table 3). A gene-based analysis (Online Methods) of ATP2B2 (P=0.0006) and SLC9A3 (P=0.0001) indicated evidence for allelic heterogeneity after comparing results with single-SNP analysis (Table 3). SLC9A3 was replicated in the French cohort (P=0.017) while ATP2B2 was not (P=0.283) (Supplementary Table 3).

Table 3

Ranked SNPs with q values < 0.05 from GWAS or GWAS-HD

			GWAS			GWAS-HD
SNP	CHR	BP	P value	q valuea	rankb	q valuea	rankb	Apicalc	GENE
rs3788766	X	115480867	1.28×10−12	0	1	0	1	0	SLC6A14
rs4077468	1	204181380	9.88×10−9	0.0018	3	0	2	1	SLC26A9
rs4077469	1	204181508	9.88×10−9	0.0018	2	0	3	1	SLC26A9
rs7512462	1	204166218	2.14×10−8	0.0023	5	0	4	1	SLC26A9
rs7419153	1	204183932	1.01×10−7	0.0091	6	0.0001	5	1	SLC26A9
rs7415921	1	204177506	4.49×10−7	0.0346	7	0.0003	6	1	SLC26A9
rs12047830	1	204183322	3.72×10−6	0.167	12	0.0022	7	1	SLC26A9
rs5905283	X	115479909	1.69×10−8	0.0023	4	0.0045	8	0	SLC6A14
rs1318819d	3	10413649	2.11×10−5	0.6309	18	0.0107	9	1	ATP2B2
rs4684689d	3	10423426	2.76×10−5	0.6828	20	0.0123	10	1	ATP2B2
rs495435d	3	10443723	4.80×10−5	0.7766	30	0.0191	11	1	ATP2B2
rs12741299	1	204181139	1.23×10−4	0.7766	81	0.0426	12	1	SLC26A9
rs1874361	1	204174809	1.31×10−4	0.7865	90	0.0426	13	1	SLC26A9
rs17563161d	5	550624	1.47×10−4	0.796	98	0.0437	14	1	SLC9A3

The q value is a genome-wide adjusted P value that controls the false discovery rate.

The rank indicates the genome-wide ordering of the SNPs based on the original association evidence (GWAS) or SFDR q value after incorporating the apical plasma membrane hypothesis (GWAS-HD). The GWAS and GWAS-HD rank results of all 543,927 SNPs are provided in Supplementary Fig. 5.

Apical indicates whether a gene was (= 1) or was not (= 0) on the generated apical membrane constituent list.

The SNPs from ATP2B2 and SLC9A3 identified by GWAS-HD; results in the French replication sample are P=0.1606, 0.1948, 0.0674, 0.0665 for rs1318819, rs4684689, rs495435, and rs17563161 respectively, where rs17563161 was imputed.

Next, restricting analysis to the 3,814 SNPs annotated to the 155 apical genes (which does not include SLC6A14), we tested the apical prioritization hypothesis as part of the GWAS-HD. Here we observed genome-wide significant evidence for association between meconium ileus and multiple constituents of the apical plasma membrane (permutation P=0.0002, testing all 3,814 SNPs jointly and not subject to multiple hypothesis testing (Fig. 2a and 2b)). Even with the exclusion of SLC26A9 (as well as SLC6A14), the apical hypothesis remained significant (P=0.0058). Thus, GWAS-HD further established the involvement of other genes coding for apical constituents despite insufficient power to detect individual SNPs, even within the context of our prioritized GWAS. For comparison, we also constructed a null hypothesis list of membrane-localized genes. As expected, the 224 GO-annotated nuclear envelope genes tagged by 3,537 GWAS SNPs showed no relationship with meconium ileus (permutation P=0.4639; Fig. 3).

Figure 3

Assessment of the nuclear envelope null hypothesis. A list of 231 genes was generated from the nuclear envelope as defined by GO annotation (GO:0005635) similarly as for the apical membrane list. In total, 3,537 GWAS SNPs are within ±10 kb of the boundaries of 224 tagged genes (NCBI36/hg18). A priori, the nuclear envelope list should not contain genes associated with meconium ileus (under the null of no association).

(a) QQ-plot of the nuclear envelope gene SNPs in the discovery sample.

(b) Statistical significance of the nuclear envelope hypothesis in the discovery sample. Statistical evaluation indicates that genes listed in the nuclear envelope are, as expected, not significantly associated with meconium ileus (permutation P=4639/10,000=0.4639).

The French cohort with genome-wide data provided independent validation of the apical hypothesis (permutation P=0.022; Fig. 2c and 2d; Online Methods). The statistical significance of this gene set (which excludes SLC6A14) in the French cohort remained after further excluding SLC26A9 (P=0.021) and then both SLC26A9 and SLC9A3 (P=0.023). Although analysis of the apical hypothesis in a larger independent cohort should be considered as part of future efforts, the replication in the smaller French cohort supports a (common) mechanism of contributing genes, even though the responsible gene variants across the two datasets may not be the same.

To determine which apical genes were driving the association, the degree of genetic heterogeneity in meconium ileus, and the common contributors across the French and North American samples, we implemented Lasso[18]. Using the North American sample, we jointly analyzed all 3,740 SNPs available in the apical genes (which include SLC26A9 and SLC9A3), and SLC6A14 (Online Methods), that were not in perfect linkage disequilibrium. Forty-eight SNPs spanning 36 different genes were retained by Lasso in the multivariate regression model (Supplementary Table 3). These SNPs jointly explained ~17% of the meconium ileus variation in the North American sample. The percentage explained by the same 48 SNPs in the French sample was 8.1% (Online Methods), presumably due to the smaller sample size and genetic heterogeneity. We then tested the significance of a score for each individual in the French cohort constructed from a weighted sum of the number of risk alleles (defined in the North American sample) across the 48 SNPs[19]. The significant association between meconium ileus and this score (P=0.0137; Online Methods) provided additional complimentary evidence of common contributors between the two cohorts, with SNPs in SLC9A3 and SLC6A14 being two specific examples.

In summary, conventional GWAS identified SNPs in SLC6A14 and SLC26A9 as significantly associated with meconium ileus. GWAS-HD single-SNP analysis identified the same SNPs in SLC6A14 and SLC26A9, as well as SNPs in SLC9A3; and multi-SNP analysis provided evidence that multiple constituents of the apical plasma membrane are collectively associated with meconium ileus, yielding considerable additional information beyond single SNP/gene associations. GWAS-HD can be applied to other Mendelian disorders, or even complex traits provided there is a biologically-based hypothesis and participating relevant genes can be compiled.

Although gene prioritization has been used in other approaches such as pathway or gene enrichment analyses[20-22], GWAS-HD involves key differences. First, in contrast to the previous inclusion/exclusion approaches where all genotyped SNPs/genes are not analyzed simultaneously, GWAS-HD performs parallel single-SNP analysis of all GWAS SNPs, and multi-SNP/gene analysis focusing on the set of SNPs/genes of interest. The prioritized single-SNP analysis interrogates all available SNPs via the SFDR control, yet enables increased statistical power for regions favoured a priori. For example, SLC6A14 would be omitted by inclusion/exclusion approaches, yet it remained the highest ranked gene for association with meconium ileus in GWAS-HD. Second, methods such as interactive pathway analysis can be restrictive because contributing genes/proteins must relate to each other via direct or indirect links, which may be disturbed when a cog component (such as CFTR) is dysfunctional in the disease state. While gene products that participate in maturation or delivery of CFTR may be contributory, consideration of local components of processes that may compensate for the ion and fluid flow disturbance in CF is enabled in the apical hypothesis. Third, distinct from an exhaustive search of all plausible interactive pathways, GWAS-HD focuses on a single biological hypothesis and provides statistical significance for all genes involved jointly, alleviating the multiple testing burden.

It should be noted that the specific statistics or models used in our GWAS-HD application such as SFDR, Lasso, and the sum and score statistics may not be the most powerful ones in any specific setting. For example, the adaptive rank truncated product statistic[23] could be used to identify the common subset of associated apical genes across two samples; and there are alternative weighting and prioritization approaches[24].

The SLC9A3 gene codes for a sodium/hydrogen exchanger that when disrupted has been shown to decrease intestinal obstructions in a CF mouse model[25]. SLC6A14 codes for a sodium- and chloride-dependant neutral and basic amino acid transporter[26,27]. SLC26A9 encodes an anion transporter, likely a chloride channel with multiple modes that include chloride/bicarbonate exchanger and sodium-anion cotransporter capabilities[28]. It has also been reported to physically interact with CFTR[29] and be influenced by CFTR activity, at least in lung-related tissues[30].

It is notable that SLC9A3 has previously been associated with infections and pulmonary function in CF[31]. In the consortium discovery sample, rs6864158 (MAF=0.43) in SLC9A3 was associated with both the lung (P=0.0003, analyzed previously[5]) and meconium ileus (P=0.0001); this provides evidence that some genes may play a role in multiple CF co-morbidities. Both SLC6A14 and SLC26A9 are robustly expressed in human lung epithelia, sweat gland, as well as intestinal epithelia as measured by RT-PCR (not shown). We anticipate that meconium ileus modifier genes may also influence early pathology in other CF-affected organs, and thus could provide important insights into the mechanisms of CF disease severity and co-morbidities.

These findings collectively have important practical implications for CF, where therapeutic strategies should consider pharmacologic modulation of epithelial function, in complementation with paradigms aimed at directly improving function or delivery of the mutated CFTR gene product to the apical membrane[32].