Genetic analysis of the Trichuris muris-induced model of colitis reveals QTL overlap and a novel gene cluster for establishing colonic inflammation

Systemic and colonic phenotyping of chronic T. muriscolitis in an F2 population of resistant and susceptible mice

An F2 inter-cross of resistant (BALB/c) and susceptible (AKR) mice was phenotyped 35 days post-infection, a time-point when chronic inflammation is established in susceptible mice (Figure 1). Colonic worm burden was not normally distributed; the majority of animals were resistant, with the largest worm burdens harboured by a small number of individuals (Figure 1A). This pattern of worm load distribution is indicative of an out-bred cohort [10]. Serum parasite-specific IgG1 (T_H2 specific) and IgG2a (T_H1 specific) was measured. Many individuals had a combination of both serotypes. To determine the predominant phenotype expressed a serum IgG1:2a ratio was calculated. A highly significant difference was observed between the mean of IgG1:2a ratio in resistant and susceptible mice (Mann Whitney U test, p < 0.0001, Figure 1B). At day 35 post-infection, 83% (110/133) of individuals with persistent worm burden demonstrated serum antibody titre IgG1 < IgG2a, indicative of a polarised T_H1 immune response. A dominant T_H2 immune response (IgG1 > IgG2a) correlated with worm expulsion. Furthermore, females demonstrated more resistance compared to males; females had significantly fewer worms at D35 post infection (Additional file 1: Figure S1A) and significantly higher IgG1:2a ratio (p < 0.0001, Additional file 1: Figure S1B) indicative of a dominant Th2 immune response.

Colonic histological assessment demonstrated persistent T. muris infection and large bowel inflammation. Mild-to-moderate inflammatory changes included: transmural tissue oedema and associated leukocytic infiltration (lymphocytes, macrophages, neutrophils); prominent mucosal and submucosal reactive lymphoid aggregates; colonic crypt hyperplasia and hypertrophy (Figure 1C-E). Significant correlation between histological parameters of inflammation (e.g. crypt length), immune response phenotype (Figure 1F: Spearman’s Rs = -0.54) and worm burden (Figure 1G: Spearman’s Rs = 0.84), were demonstrated. 98.5% of mice with persistent helminthosis demonstrated colonic inflammatory changes.

Whole genome Linkage analysis

Of particular interest, Tm3 (92.4-118.3 Mbp, chromosome 3) demonstrated complete overlap with a susceptibility locus identified in three unrelated murine models of spontaneous experimental colitis: G-protein alpha inhibitory 2 chain knock-out (Gnai2^-/-) mice (Gpdc1 locus) [11]; C3H/HeJBir IL10-deficient mice (Cdcs1 locus)[12]; T-bet^-/-Rag2^-/- double-deficient mice that resemble ulcerative colitis (TRUC) (Cdcs1 locus) [13]. The Cdcs1 region has been shown to contain at least three distinct regions [14, 15]. Here, we show complete overlap of Tm3 with Cdcs1.1 (Figure 2), a region shown to contribute strongly to the severity of colitis in C3H/HeJBir mice [14, 15].

Prioritization of QTL candidate genes via pathway-driven workflow analysis

Schematic representation and key stage data from each of the qtl_to_pathway [16] and refseq_ids_to_pathways [17] workflows are shown (Figure 3). In total, 1419 genes were identified within the 7 T. muris QTL. Genes attributed to one or more KEGG biological pathways were determined (Figure 3.1). Simultaneously, of 5476 genes with significant transcriptional differences during T. muris infection [9], 2504 genes displayed either an up-regulated or down-regulated change in expression of ≥ 1.4 fold (i.e. a 40% or greater change in expression over naïve controls, Figure 3.2). Biological KEGG pathways associated with significant gene expression data were determined. In total, 1158 of these genes (46%) were involved in 204 separate KEGG pathways (Figure 3.3).

The cross correlation of functional pathways containing QTL genes and genes demonstrating significant expression were identified, linking genotype and phenotype trait interactions (Figure 3.3). Finally, polymorphic genes between parental AKR and BALB/c mice were identified within each locus (Figure 3.4).

As an example, 344 Ensembl gene ID’s were detected within Tm3. Of these, 97 (28%) were designated as functionally important within molecular interaction networks, as assigned by the KEGG pathway database. Significantly expressed microarray genes were similarly assigned biological pathways. For Tm3, the cross correlation of common pathway data and the exclusion of any gene which lacked SNPs between parental strains, identified 17 Quantitative Trait genes (Figure 3.4, Column D). In comparison, 61 KEGG-assigned polymorphic genes did not demonstrate any change in transcriptional activity (Column C). Of the genes yet to be allocated a KEGG pathway, 16 of 191 genes displayed significant transcription (Column B). The same process was undertaken for all 7 QTL (Figure 3).

Chromosome 3 candidates

Candidate gene validation

Quantitative PCR analysis was undertaken independently in infected parental strains (days 0, 7, 14, 21, and 35 post-infection) to validate microarray data (Tables 2 and 3).

For gene candidates, Vav3, Ptpn22, FcgR1 and S100a10, qPCR corroborated up-regulated expression found in susceptible AKR on day 35 post infection (Figure 4). Likewise, the down-regulation of Hmgcs2 was confirmed by qPCR. Microarray analysis of Ctss and RORc demonstrated down regulation of colonic gene expression in chronically affected individuals, which could not be replicated (data not shown for RORc).

Animals

Mice were housed with free access to food and water under specific pathogen free conditions. All experiments were performed under regulations of The UK Home Office Animals (Scientific Procedures) Act of 1986.

For QTL analysis, AKR/OlaHsd (susceptible, hereafter referred to AKR) and BALB/cOlaHsd (resistant, hereafter referred to BALB/c) mice (Harlan Olac, UK) were interbred. To generate an F1 population of mice, equal numbers of AKR males vs BALB/c females (F1a offspring), and AKR females vs BALB/c males (F1b offspring) were mated. At least 50 breeding-pairs of F1-mice were then interbred. All F1 vs F1 breeding was performed over the same time-period. To maintain genetic balance, F1a males were bred with F1a and F1b females; and, F1b males with F1a and F1b females. A single generation of 307 F2 mice (male and female) was created for study. All F2 mice were infected at the same time with T. muris ova at 6-8 weeks of age.

Parasites

Trichuris muris parasites were harvested and ova collected and maintained as previously described [18]. All infected mice received 300 T. muris ova in distilled water (200 μl) by oral gavage.

QTL phenotyping

Phenotypic analysis was performed for all 307 F2 mice. Day 35 post-infection, serum samples and intestines were taken at autopsy. Resistance (0 worm load) and susceptibility (>0 worms) were defined. All worm counts were performed by a single investigator over 1 week from caeca frozen at autopsy. This method of storage and counting is used routinely for large experiments and does not affect quantification. Parasite-specific antibody ELISA was performed as described previously [44], using in-house T. muris excretory-secretory (ES) protein. T. muris specific IgG1 (T_H2 specific, driven by IL4) and IgG2a (T_H1 specific, driven by IFNγ) optical density (OD) was measured simultaneously for all samples. All 307 serum ELISAs were performed in one run. For histology, 0.5cm of whole colonic segments from the proximal ascending colon was snap-frozen, thawed in 4% paraformaldehyde, paraffin embedded, and 5μm transverse tissue sections stained with Haematoxylin and Eosin (H&E) simultaneously. 50 randomly assigned colonic specimens were assessed. Proximal colonic specimens were scored according to colonic crypt length (μm), immune cell infiltration and tissue inflammation by a single investigator. Colonic crypt length per individual was taken as a mean across at least 20 crypt units and 3 separate sections. Crypt units were measured using Image-J software [45]. Spearman’s rank correlation coefficient was performed to measure the statistical dependence of (a) worm count and colonic crypt length variables, and (b) IgG1:IgG2a ratio and crypt length variables.

DNA extraction

DNA was isolated (Promega Wizard DNA isolation kit) from tail snips digested in proteinase K digestion buffer (20 mg/ml). DNA concentration was determined by Nanodrop spectrophotometer and then stored at -80°C until analysis.

Linkage Map

165 polymorphic murine microsatellite markers distinguishing between AKR and BALB/c were selected [46]. Whole genome coverage was 85% and median inter-marker distance 12.3 cM. Conversion of marker positions from recombination fraction (cM) to physical position (Mb) was achieved using the Ensembl database [47].

Microsatellite amplification and genotype analysis

Forward polymerase chain reaction (PCR) primers were fluorescently labelled with 6-FAM, HEX or NED (MWG Biotec, Applied Biosystems). 25 ng of genomic DNA was used for each marker. Semi-automated analysis of genotypes on pooled panels of PCR products was performed using an Applied Biosystems 3100 Capillary sequencer with Genescan analysis and Genotyper software.

QTL analysis

IgG1:IgG2a ratios were log₁₀ transformed to achieve parametric distribution. Median, mean and kurtosis values were calculated using QStat, Windows QTL Cartographer 2.5. Normalised data were analysed using multiple interval mapping to optimize and refine QTL positions. A genome-wide permutation test (1000 repeats) determined thresholds for significance; a logarithm of odds (LOD) score of 4.0 or a p value of <5.2×10^-5 was considered significant. A LOD score of 2.5 or p value of 1.6×10^-3 was considered suggestive of linkage according to published guidelines [48]. All significant LOD scores were confirmed by 1-way ANOVA with pairwise comparison, using the Bonferroni correction method. Kruskal-Wallis analysis was used for worm burden and IgG2a data, and converted to an LOD score [49].

Genome-wide colonic transcriptional activity of parental murine strains

Naïve and infected 6-to-8 week old male AKR and BALB/c mice (Harlan Olac, UK) were monitored through to day 35 post-infection (n = 6, 3 experimental replicates for each conditional cohort) as described previously [9]. 3 replicate pooled samples of colonic RNA (ascending colon) were generated for each experimental group. Whole transcriptome microarray expression analysis (Affymetrix Genechip Mouse Exon 1.0 ST Array®) and bioinformatic analysis was performed. The entire genome-wide expression dataset was used for subsequent analysis [9, 50].

The in-silicoprioritization of QTL candidate genes

The use of workflows in the analysis of large-scale genomic data provides a systematic and un-biased mechanism for hypothesis generation [51]. Previously constructed workflows were re-used for the analysis of QTL and gene expression data, to identify biological pathways which correlated with Trichuris muris infection. The identification of candidate genes underlying each QTL was carried out by firstly determining the precise co-ordinates of each genetic marker (Mbp) (Table 1). Each QTL was subsequently entered into the workflow qtl_to_pathway (Additional file 1: Figure S2 [16]). Genes located within each QTL were annotated with additional accession number identifiers (including UniProt ID and Entrez Gene IDs), in order to cross-reference Ensembl database identifiers to KEGG (Kyoto Encyclopaedia of Genes and Genomes) [52] pathway identifiers. As a result, annotated biological pathways were extracted from the KEGG database for inclusion in further analysis.

In parallel, differentially expressed genes identified from the T. muris microarray study [9] were analysed using the refseq_ids_to_pathways workflow (Additional file 1: Figure S3 [17]). This workflow required preliminary analysis of the gene expression data [9] (Partek Genomic Solution version 6.5, 2009, Partek, USA) and conversion of Affymetrix probe-set identification markers to their recognised NCBI RefSeq identification code (refseq ids). An identical process to that of the qtl_to_pathway workflow for gene annotation was then carried out.

The mapping of gene expression data to KEGG highlighted biological pathway activity in the pathogenesis of colonic disease. All genes with significant transcriptional differences between resistant and susceptible strains, in naïve and infected states (ANOVA, factor interaction, p <0.05), were included for analysis. To identify cis-QTL genes of biological relevance to phenotype, those genes with a higher degree of over/under expression (Fold Change ≥ +/-1.4 over naïve levels) during chronic T. muris intestinal inflammation, were used in the workflow analysis (see Figure 3).

The workflow common_pathways (Additional file 1: Figure S4 [53]) was used to identify candidate pathways containing differentially expressed genes within a QTL, in order to obtain an overall view of the mechanisms which may be influencing the expression of the phenotype.

Additional text mining was used to prevent potential candidate genes which lacked KEGG pathway annotation from being discarded. Transcribed QTL genes were analysed using a text mining workflow (Additional file 1: Figure S5 [54]). Briefly, published abstracts were identified from a PubMed search using the term “(“Colitis” AND “Inflammation”) AND (“Human” OR “Mouse”)”. All scientifically relevant keywords contained within individual abstracts were extracted, constructing a phenotype concept profile and allowing the calculation of inverse document frequency (IDF) scores ie a score relating the number of resulting documents which contained the keywords in question. In parallel, abstracts pertaining to selected genes were similarly recorded. The identification of phenotype keywords within individual gene abstracts allowed for the generation of a cosine vector score for each gene ranging from +1 to -1 (+1 = causation of phenotype; 0 = unknown association with phenotype; -1 = preventative of phenotype). Ranked by their cosine vector score, the association with phenotype of a particular gene was displayed. Similarly, individual phenotype keywords were also ranked according to the IDF scores, identifying possible correlations between each gene and the phenotype. All data regarding text mining and workflow approaches are published online [55].

Only QTL genes known to possess SNP variation between parental AKR and BALB/c [56] were subject to further analysis.

Independent replication of candidate gene expression by qPCR (Tm3)

Infected parental strains AKR and BALB/c (Harlan Olac, UK) received 300 T. muris ova by oral gavage. Mice were culled days 0 (naive), 7, 14, 21 and 35 post-infection for analysis (n = 3 for each cohort). mRNA was extracted from 0.5 cm of whole colonic tissue segments, from the ascending colon, according to manufacturer’s instruction (TRIZOL®, Invitrogen). cDNA was synthesised. A full list of gene primers (Eurofins-MWG-Operon, Germany) and their sequences are provided (Additional file 1: Table S1). Samples were quantitatively analysed using KAPA SYBR FAST qPCR Master Mix (Kapa Biosystems Inc., USA) and a Bio-Rad MyIQ™ PCR detection system (Bio-Rad IQ5 optical system software, version 2; Bio-Rad Laboratories Inc.,©). Three replicate cDNA samples were run at a 1:20, a 1:100, and a 1:500 dilutions for each time-point. Threshold cycles were calculated; gene detection within the three serially diluted samples was standardized, and then normalized against housekeeping gene beta-actin (Act-b). Relative fold change in gene quantity was calculated using naïve resistant mice as a reference.