Variable exon usage of differentially-expressed genes associated with resistance of sheep to Teladorsagia circumcincta

Graphical abstract


Introduction
The abomasal nematode, Teladorsagia circumcincta is one of the most common parasites of sheep in cool temperate regions (Urquhart et al., 1987) and is a major drain on the economics of sheep production (Nieuwhof and Bishop, 2005). Several studies have described genetic variation for parasite susceptibility within and between sheep breeds (Dominik, 2005;Mugambi et al., 1997;Terefe et al., 2007), and that faecal egg count (FEC) is a heritable characteristic (Davies et al., 2005) and can be used as a parameter for selective breeding for parasite resistance (Davies et al., 2005;Windon, 1990;Woolaston, 1997). Sheep selected for low FEC control infection by the acquisition of anti-parasite antibody (Stear et al., 1999) driven byTh2 cytokines (Craig et al., 2014), and there is a significant genetic relationship between IgA antibody levels and the ability to control infection Strain et al., 2002). In contrast, sheep selected for high FEC produce little IgA antibody but instead generate an inflammatory Th1/Th17 response that fails to control infection and egg production (Gossner et al., 2012).
FEC and antibody levels have both been used as selectable markers of resistance (McRae et al., 2014;Sayers and Sweeney, 2005;Stear et al., 2006) but the identification of genetic markers has the potential for marker assisted selection (Davies et al., 2006). Both quantitative trait loci mapping and genome-wide association studies have been used to identify selectable genetic markers (Davies et al., 2005;Dominik, 2005;Riggio et al., 2013). More recently high throughput microarrays (Andronicos et al., 2010;Diez-Tascon et al., 2005) and digital gene expression  have been used to identify potential candidate genes for resistance to sheep gastrointestinal nematodes including Haemonchus contortus (Rowe et al., 2008) and T. circumcincta (Gossner et al., 2013;Knight et al., 2011). The aim of these studies was to identify genes that were differentially-expressed in pre-exposed (immune) vs. naïve or resistant vs. susceptible sheep and to try and explain the molecular basis of differential polarization of the immune response associated with each phenotype.
Microarray analysis of abomasal lymph nodes from infected sheep, selected on the basis of FEC, highlighted genes and physiological pathways related to the polarization of T cells and control of inflammation as being associated with resistance (Gossner et al., 2013). These data were obtained using the Affymetrix Ovine Gene 1.1 ST whole-genome array consisting of probes for all exons of each annotated transcript in the sheep Oar v2 genome assembly. The data were originally analysed at the gene-level, with exons summarized to genes, but this analysis also revealed that individual exons for some genes were not equally expressed.
Alternative splicing is a post-transcriptional mechanism for regulating gene expression; and within a gene the variable usage of (or part of) an exon or intron creates multiple distinct transcripts (Blencowe, 2006). Alternative splicing can modify coding sequence, but within the 5 and 3 UTRs it can also alter miRNA-binding and regulatory motifs and influence transcript expression levels or translation, and therefore can have a profound effect on the identity or expression levels of a protein (Martinez and Lynch, 2013).
The focus of this current study was to examine the differential usage of exons of selected genes in the resistant and susceptible sheep, as identified by microarray analysis (Gossner et al., 2013). The eight genes were chosen on the basis that they were differentially-expressed in resistant and susceptible sheep and have a role in T cell biology and inflammation and included arachidonate 15-lipoxygenase (ALOX15) and interleukin 13 (IL13), both of which were significantly increased in resistant animals and had variable exon usage. The microarray data were first validated by fold-change RT-qPCR for two exons of each selected gene; copy number RT-qPCR was then used to test the hypotheses that resistance/susceptibility to T. circumcincta is associated with differential exon usage within the selected genes and that individual exon usage correlates with selected quantitative parameters of resistance.

Animals and tissues
Details of animals, infection protocols, phenotypic parameters and population genetic analyses have been described previously (Beraldi et al., 2008). Briefly, ∼3 months old Blackface female lambs originated from a flock used for quantitative trait analyses (Davies et al., 2006); all animals were housed in worm-free conditions. Ten sheep were sham infected controls and 45 sheep were infected experimentally with ∼2300 L3 T. circumcincta larvae three times a week for three months and sacrificed two days after the last infection. At post mortem ten of the infected group had no adult worms in the abomasal contents, while the other infected animals had adult worm count (AWC) within the abomasum at post mortem ranging from 80 to 11300 (Table S1). Animals were ranked (1-45) according to AWC and FEC (Beraldi et al., 2008). The relative concentration of serum IgA anti-T. circumcincta antibody was also measured for all these animals; and it was determined that both AWC and FEC were significantly negatively correlated with antibody levels and with body weight (BW) (Beraldi et al., 2008). The 7 most resistant (R) sheep (infection rank 1-7) had no detectable AWC or FEC, high IgA antibody levels and high body weight (BW). The 7 most susceptible (S) sheep (infection rank 39-45) were those with the highest AWC (mean 6000, maximum 11300), high FEC (mean 414, maximum 950), low IgA antibody levels and low body weight. Animal experiments were approved by University of Edinburgh Ethical Review Committee and conducted under an Animals (Scientific Procedures) Act 1986 Project Licence. Abomasal (gastric) lymph nodes (ALN) were removed at post mortem and stored at −80 • C in RNAlater (Ambion, Huntingdon, UK).

RNA extraction and cDNA synthesis
Total RNA was isolated from 0.02 g of tissue using the Ribopure Kit (Ambion, UK) according to the manufacturers' instructions. Contaminating DNA was removed by On-column PureLink ® DNase I treatment (Ambion). The quantity, quality and integrity of the RNA samples were determined using a NanoDrop ND-1000 spectrophotometer (Labtech International Ltd.) and Agilent 2200 TapeStation system (Agilent Technologies); all had an RNA Integrity Number of >7.3. cDNA synthesis from 1 g RNA was by SuperScript TM II Reverse Transcription Kit (Invitrogen, UK) using oligo-dT( 15 ) primer and RNaseOUT (Invitrogen).

Cloning of ovine exons
The sequences of differentially-expressed probes of the Affymetrix Ovine Gene 1.1 ST Array (http://www.affymetrix.com/ analysis/index.affx) were used to obtain full exon sequences using NCBI-BLAST against the sheep genome assembly, Oar v3.1 (http://www.livestockgenomics.csiro.au/sheep/oar3.1.php/). Exon-specific primers were selected using Primer-BLAST (http:// www.ncbi.nlm.nih.gov/tools/primer-blast/) and reanalysed using Net Primer (www.premierbiosoft.com/netprimer/). Primers used for RT-PCR and RT-qPCR are shown in Table S2A. The location of the primers and their relationship to the locations of the Affymetrix Ovine Gene 1.1 ST whole-genome array probe sets is shown in Table  S3. RT-PCR used the FastStart Taq DNA Polymerase Kit (Roche) following the manufacturer's protocol. PCR products were analyzed by agarose gel electrophoresis, visualized by GelRed/UV transillumination, purified using MinElute PCR Purification Kit (Qiagen), ligated into pGEM-T Easy vector (Promega) and transformed into JM109High Efficiency Competent Cells (Promega). A random selection of color-screened clones were sequenced (Edinburgh Genomics; https://genomics.ed.ac.uk/) with SP6 and T7 primers using the BigDye ® Terminator v3.1Cycle Sequencing Kit (Applied Biosystems, UK).

Quantitative real-time PCR analysis
qPCR was performed using 1 l template cDNA or linearized plasmid DNA (for the copy number analysis), 7.5 l FastStart Universal SYBRgreen Master (Rox) 2× concentration (Roche), 0.1--0.3 l primers and nuclease-free water to 15 l final volume. All reactions were prepared using a CAS-1200 TM Precision Liquid Handling System and performed on the Rotor-Gene Q (Qiagen). Cycle conditions were: 95 • C for 10 min and then 40 cycles of 95 • C for 10 s, annealing (Table S2A) for 15 s and 72 • C for 30 s followed by melt curve analysis. Optimised RT-qPCR primers had an efficiency between 95 and 105% and R 2 of 0.98 and 0.99. Relative analysis was performed on the 7 most resistant sheep (ranked 1-7) and the 7 most susceptible sheep (ranked 39-45) (Table S1). Absolute expression (copy number) analysis was performed on all 45 infected animals. Transcripts were quantified in cDNA from three separate RT-qPCR reactions for each biological sample; each cDNA sample was assayed in triplicate with GAPDH housekeeping and no-template controls included in all runs.
Relative gene expression levels were calculated in GenEx 5.3.4.157 (MultiD Analyses AB, Sweden) using the comparative 2 − ( Cq) method and normalized to the geometric mean of GAPDH and SDHA. Fold changes were calculated from Cq values using GenEx. To derive the copy number of the target sequence in all 45 sheep, a standard curve (linearized sheep ALOX15 or IL13 plasmid DNA) was used with a dynamic range that spanned at least five orders of magnitude. Copy numbers were calculated from Cq values using the following formula: molecules per g = ((1 × 10 −6 )/(M g/mol)× [6.03 × 10 −23 molecules/mol]) M = size of plasmid × 660 g/mol per bp. The expression levels were normalized by dividing the copy number derived from the standard curve by the calculated normalization factor for each individual sample.

Sequencing 5 intronic regions
Genomic DNA (gDNA) was extracted from the ALN of 6 resistant (infection rank 1, 3, 4, 5, 6, 7) and 6 susceptible (infection rank 39, 41, 42, 43, 44, 45) sheep using the Wizard ® SV Genomic DNA Purification System (Promega). 10-20 mg of ALN tissue was incubated in the digestion solution at 55 • C with 300 rpm shaking overnight and then processed following the manufacturer's protocol; 200 ng gDNA per PCR reaction was used as a template for cloning and sequencing of ALOX15 and IL13 promoter regions (primers and parameters Table S2B).

Statistical analysis
Relative gene expression levels were analyzed statistically in GenEx using an unpaired, 2-tailed t-test to determine the difference between groups. ALOX15 and IL13 copy number results were analyzed in GraphPad Prism v 5 (Graph Pad Software, USA). The data were grouped into resistant, susceptible and intermediate (n = 15 per group) and a Kruskal-Wallis test was performed to determine overall significance, with Dunn's multiple comparison test within the Kruskal-Wallis to determine significance between groups. The correlations between transcript levels and quantitative phenotypes were analyzed with Spearman's rank correlation coefficient (r s ). p-values less than 0.05 were considered statistically significant.

RT-qPCR validation of array analysis
Transcriptome analysis of the ALN has been described previously (Gossner et al., 2013). In brief, cDNA from the abomasal lymph node of 7 resistant, 7 susceptible and 7 uninfected sheep was hybridized to Affymetrix Ovine Gene 1.1 ST Array and analysed, with exons summarized to genes, using the mean expression of all the exons of a gene. These data and protocols are available at ArrayExpress accession number E-MTAB-1580. In this current study the same dataset was used but without summarizing to genes. Of the 165740 exons interrogated, 1196 were significantly differentially-expressed within their gene in the R vs. S comparison, with a fold change ≥2 and p-value ≤0.05 (one-way between-subject ANOVA, unpaired); 930 exons were increased and 266 exons were repressed. Eight genes were selected, based on each being significantly differentially-expressed in the R vs. S comparison, their role in T cell biology and inflammation as well as significant differential exon usage within the gene (Table S4). RT-qPCR was developed for two exons of each gene including one with significant differentialexpression.
Comparison of the array and RT-qPCR fold change analysis for the exons of the selected genes is shown in Table 1. Significant differential exon expression was validated, by RT-qPCR, only for ALOX15 and IL13. Exons 9 and 14 of ALOX15 were 22.02 and 7.44 fold increased respectively, in the resistant compared to susceptible sheep as assessed by array, and 5.3 (p = 0.001) and 1.5 fold (p = 0.03) by RT-qPCR. Exons 1 and 4 of IL13 were 3.6 and 6.91 fold increased by array and 3.77 (p = 0.0002) and 5.16 fold (p = 5.2E − 6) by RT-qPCR. When the data were summarized to the whole gene, ALOX15 was 5.81 fold higher, and IL13 was 3.17 fold higher in the resistant sheep (Gossner et al., 2013). Differential exon expression of the other six genes was not confirmed by RT-qPCR, although all but CD109 demonstrated significant differential expression for at least one exon in the R vs. S comparison.

Association of exon copy number and quantitative phenotype
ALOX15 and IL13 were chosen for copy number measurement in the ALN in all 45 infected animals. The expression levels of ALOX15 exon 9 was significantly higher than exon 14 (mean 11313 ± 4547 and 7729 ± 2360 copies per g total RNA respectively, p = 0.0003). Similarly IL13 exon 4 was expressed significantly higher that exon 1 (152.377 ± 81.222 and 72.406 ± 43.858, p < 0.0001). The data were further analysed by dividing the 45 animals into three groups; the fifteen most resistant sheep with mean AWC of 59 and mean FEC of 1.67; a susceptible group of fifteen sheep with mean AWC of 5167 and mean FEC of 288; and the fifteen intermediate sheep with mean AWC of 1508 and mean FEC of 82 (Table S1).

Sequence of 5 intronic or UTR regions of ALOX15 and IL13
Polymorphisms in intronic promoter or 5 UTR regulatory elements (Black, 2003), possibly associated with differential exon usage, were examined by sequencing approximately 1000 bp 5 to the translation start site of both ALOX15 (LN864492) and IL13  (LN864491) from 6 R and 6 S sheep. Ten single base variations were identified in the ALOX15 region (g.66T>C; g.227C>G; g.246C>G; g.606T>C; g.654G>A; g.656G>C; g.697T>C; g.836T>C; g.960G>A; g.1036A>G) but none segregated with phenotype. An initiator codon variant (A/G) was also identified in two animals, which has already been identified in other sheep breeds (NCBI dbSNP; rs422045752). The only variation found with IL13 was the insertion of three nucleotides (g.674 675insGAA) 222 nucleotides 5 to the translation start site in one animal (lamb 110, rank 7), which was heterozygous at this locus.

Discussion
The resistance of sheep to the gastrointestinal nematode parasite T. circumcincta has a major host genetic component (Murphy et al., 2010;Riggio et al., 2013) and differential susceptibility to this parasite is clearly associated with T cell activation and inflammation (Gossner et al., 2013;Gossner et al., 2012). In addition, mouse models of nematode infections show that variations of genes associated with Th1, Th2 and Th17 T cell maturation are involved in the genetics of resistance (Filbey et al., 2014;Finkelman et al., 1997). Most studies in sheep have focussed on the important role for transcription regulation and primary sequence differences of these genes (Gossner et al., 2013;Riggio et al., 2013); however this current study concentrates on differential exon usage as a source of gene variants. Alternative splicing is a significant cause of protein variation (Martinez and Lynch, 2013) especially in genes expressed in the immune system (Lynch, 2004) including cytokines and receptors involved in T cell polarization (Gaudreau et al., 2012;Martinez et al., 2012); and the original selection of the eight genes was made on the basis that they were differentially-expressed in resistant and susceptible sheep (Gossner et al., 2013) as well as being involved in T cell biology and inflammation.
Fold-change RT-qPCR resulted in the validation of the microarray data of two of the eight selected genes, ALOX15 and IL13. IL13 plays a major role in resistance to gastrointestinal nematodes. In T. circumcincta infected sheep it was the most significant up-regulated gene in the R vs. S comparison by both microarray and RT-qPCR array (Gossner et al., 2013). It is also the dominant cytokine in the control of murine Nippostrongylus brasiliensis infection  and seems to be equally involved with IL4 in the development of protective immunity to Trichuris muris and Heligomosoides polygyrus . IL13 is largely a product of antigen-activated Th2 cells but is also synthesized by mast cells and eosinophils (Gessner et al., 2005). IL13 shares many biological functions with IL4, including the regulation of Th2 development and consequently antibody production and heavy chain switch; it also inhibits inflammatory cytokine production and promotes mast cell proliferation (Minty et al., 1993) and tissue remodelling after parasite-induced injury (Wynn, 2003). IL13 is encoded by four exons on the minus strand of ovine chromosome 5 (NC 019462.1); exon 4 encodes 108 bp at the 3 end of the coding region (19.262.871-19262764) and the 862 bp 3 UTR (19.262.763-19261902). The consistent differential expression of both exons between the R and S groups, and the highly significant correlation of exon 4 levels and sheep rank phenotypes (AWC and FEC) as well as BW, argues against a simple technical explanation for differential exon usage. Expression levels of exon 4 are likely to represent full length and functional IL13 transcripts, which may explain why quantitative levels of IL13 exon 4 are significantly positively correlated BW and consequently negatively correlated with AWC and FEC. The lack of correlation between IL13 exon 4 and IgA (r 2 0.08, p = 0.31) may be due to the fact that the IL13 levels are transcripts at a single time point (12 weeks post-infection), but that serum IgA had accumulated over time.
Another major function of IL13 (and IL4) is the promotion of ALOX15 expression (Nassar et al., 1994) and the consistent upregulation of all the ALOX15 exons in resistant sheep could be explained by the high levels of these cytokines in those sheep. ALOX15 is one of a range of at least six isoenzymes that oxygenate arachidonic acid leading to the production of the lipoxins (Samuelsson et al., 1987). These function to inhibit leukotrienes and therefore are broadly anti-inflammatory. In addition ALOX15 stimulates the clearance of apoptotic cells and promotes tissue healing (McMahon et al., 2001), which is seen in resistant sheep (Gossner et al., 2012). ALOX15 has not been annotated in the Oar v3.1 genome assembly but analysis using the Bos taurus ALOX15 sequence (NM 174501.2) identifies that the sheep gene is encoded by 14 exons on the plus strand of chromosome 11 (NC 019468.1). Exon 9 encodes 90 bp within the coding region (1247-1337 from the translation start site, 26326237-26326327) and exon 14 encodes the terminal 180 bp of the coding region and the complete 3 UTR (26327916-26328935). As with IL13, it is likely that the ALOX15 transcript that includes exon 9 encodes the full length and functional enzyme, and therefore there is a highly significant relationship between exon 9 expression levels and phenotypic parameters of parasite resistance. One transcript and three predicted transcript variants of ALOX15 have been described in cattle, with the X2 variant (XM 005220190.2) having a deleted exon 7. However, no functional consequences of these variants have been described.
Control of exon usage and differential splicing can be associated with intronic regulatory elements (Black, 2003), often in the proximal 1000 bp 5 to the translation start site. These regions were sequenced for both ALOX15 and IL13 from 6 resistant (rank 1-6) and 6 susceptible (rank 40-45) sheep; and although variants were identified for both genes, none segregated with resistance and susceptibility phenotype.
In conclusion, we have validated differential exon usage of two genes that were expressed in the ALN of sheep showing differential resistance to the abomasal nematode T. circumcincta. Both ALOX15 exon 9 and IL13 exon 4 were significantly increased in resistant animals and expression levels of these exons were negatively correlated with quantitative phenotypic traits, including AWC and FEC. Consequently, they represent potential markers for selection for resistance to a common and economically important gastrointestinal parasite.