Dynamic Regulation of Tandem 3′ Untranslated Regions in Zebrafish Spleen Cells during Immune Response

Poly(A) tails are found on the 3′ end of nearly all mature eukaryotic mRNA and are important for the stability, translation, and transport of mRNA (1, 2). The poly(A) tail is acquired through the process of polyadenylation, and pre-mRNAs may be subjected to alternative polyadenylation (APA) due to multiple polyadenylation signals, resulting in different mRNA isoforms. Investigation of existing mRNA databases shows that more than half of human genes exhibit two or more 3′ untranslated region (UTR) isoforms (3, 4), and high-throughput sequencing-based analyses of human, mouse, and zerbrafish transcriptomes have revealed that approximately half of all expressed genes undergo APA (5–7). In recent studies, APA has been suggested to play an important role in tumorigenesis, development, and cell differentiation (8–10). APA may alter the coding sequence when its signals are located in upstream introns or exons. APA involving multiple poly(A) signals in the terminal exon has no impact on the protein sequence, but instead it generates mRNA isoforms with tandem 3′ UTRs of differing lengths (3, 11). Generation of tandem 3′ UTRs is the most common consequence and accounts for more than half of all APA events (7, 12). Shortened tandem 3′ UTRs have been shown to be preferentially expressed in proliferating cells (13, 14), whereas lengthened tandem 3′ UTRs were more commonly observed in differentiating cells (7, 15). After activation, murine CD4⁺ T lymphocytes express mRNAs with shortened 3′ UTRs, and globally decreased expression of distal 3′ UTRs has been observed in human T cells, B cells, and monocytes upon stimulation with a variety of activating factors. The shortened 3′ UTRs were reported to possess fewer microRNA (miRNA) target sites (13). This investigation was based on in vitro experiments and only focused on only a limited number of genes with two known APA sites. To our knowledge, so far there are no studies addressing APA-mediated immune regulation, and its role in immune response is largely unknown.

The zebrafish is a valuable model for immunological studies owing to the high similarities of its innate and adaptive immune systems to those of mammals. Nearly all immune cells and molecules of mammals have equivalents in zebrafish (16–21). With its fully developed immune system, through the use of genomic sequencing, genetic screening, and real-time visualization, zebrafish have become an ideal model for studying in vivo host immune responses to infection (22).

Experimental infection models have been established to study both human and fish pathogens, including bacteria, viruses, and parasites (19, 23–25). The acute phase response of zebrafish upon infection with Staphylococcus aureus has been shown to exhibit close similarity to that of mammals, not only in molecular constitution, but also in regulation of gene expression and functional strategies (24). These previous observations make zebrafish an ideal model for comprehensive study of the dynamic regulation of the complexities of immune response to infections in vivo.

To understand the involvement of APA in immune responses in vivo, the sequencing APA sites (SAPAS) strategy (26) was used to create a genome-wide profile of the APA-derived 3′ UTRs in the acute immune response. Our analysis revealed genes involved in immune signaling pathways to be globally activated, with the TLR signaling, apoptosis, and JAK-STAT signaling pathways displaying the greatest number of APA genes. In the APA-regulated immune signaling pathways, most APA genes were cytokines, cytokine receptors, and signal transducer genes. Our comprehensive analyses of APA genes indicate that APA may play an important role in regulation of the immune response against infection in vivo.

Wild-type zebrafish were reared and propagated in recirculating systems at 28°C. Adult zebrafish (∼5 mo old) were transferred to an isolated flow-through system and allowed to acclimate for 2 d prior to infection.

Overnight cultures of S. aureus were diluted 1:50 in Luria–Bertani broth, grown at 37°C with shaking, and harvested in the logarithmic phase of growth. The cells were washed twice in sterile PBS and resuspended and diluted to the appropriate concentration in sterile PBS. The final concentration of bacteria was determined by plating serial dilutions on solid media.

Adult zebrafish were anesthetized in 160 μg/ml tricaine and injected i.p. with 4 × 10⁶ clonal formation units S. aureus. Fifteen hours later, spleens from 20 bacterially challenged fish or 20 naive fish were harvested and processed to prepare bacterially challenged (BC) and unchallenged (UC) cDNA libraries.

The SAPAS libraries were prepared as described previously (26). Briefly, total RNA was extracted from samples using RNeasy Mini kits (Qiagen). Approximately 10 μg total RNA was fragmented by heating, and the first-strand cDNA was synthesized by reverse transcription with anchored oligodeoxythymine primers and a 5′ template switching adaptor. PCR was then performed to amplify the cDNA, and the number of cycles was optimized to ensure that the double-stranded cDNA remained in the exponential phase of amplification. After the PCR, fragments of 200–300 bp in size were excised from a PAGE gel and purified. The final pooled fragments were quantified and sequenced from the 3′ end with Solexa.

Analyses of the SAPAS-generated database were performed as previously described (6). Briefly, the SAPAS-generated raw reads were trimmed, filtered, and mapped to the zebrafish genome (Zv9 assembly, GCA_000002035.2, http://www.ensembl.org/Danio_rerio/Info/Index). After internal priming filtering, the resulting uniquely mapped reads (called effective reads) were clustered to define poly(A) sites. The total read counts from each sample were normalized to the libraries with fewer effective reads (e.g., the BC library), and poly(A) sites with two or more normalized reads were used for further analysis. A poly(A) site that overlapped with an annotated 3′ UTR from Ensembl was defined as a tandem poly(A) site. We used a linear model to test APA switch events (26), and a gene with an absolute value of test-resulted tandem 3′ UTR switching index of >0.1, and a p value smaller than the threshold corresponding to a BH-sense false discovery rate of 0.05 was defined as an APA gene. To analyze the relative gene expression level, the effective SAPAS reads were normalized to reads per million reads. For analysis of expression level changes between two samples, we calculated the logarithm (to the base 2) of the ratio of normalized reads for each gene. For cases in which a gene had no read in a particular sample, we assigned one read to that gene as its expression level for that particular sample. The common 3′ UTR was defined as the shared portion of the 3′ UTR, and the extended 3′ UTR was defined as the ‘‘altered’’ portion of the 3′ UTR that differs between two samples. For APA genes with exactly two tandem poly(A) sites, the common 3′ UTR is defined as the 3′ UTR between the stop codon and the proximal poly(A) site, and the extended 3′ UTR as the 3′ UTR between the proximal poly(A) site and the distal site (26). For searching motif, we considered all possible 4- to 7-mers and searched separately the common and extended regions of each APA gene. WWTT(ATTTA)TTWW motif allowing two mismatches in either side in the common and extended 3′ UTR of APA genes, AU-rich element (ARE) scores were calculated with the Web-based version of AREScore (http://arescore.dkfz.de/arescore.pl) (27). The p values were obtained using a Fisher exact test with a Bonferroni correction. A Gene Ontology enrichment analysis was performed on the DAVID Web site (28).

Thirty novel APA sites with relatively high expression levels were chosen, and 3′ RACE was performed with gene-specific primers using the 3′-Full RACE Core Set version 2.0 (Takara) according to the manufacturer’s instructions. The PCR products were sequenced with an Applied Biosystems 3730 machine. Twelve genes with significant length differences in their 3′ UTRs between the UC and the BC libraries were chosen for quantitative real-time RT-PCR (qRT-PCR). The poly(A) sites of these genes were divided into two supersites (the proximal and distal sites), and the region upstream of the supersites was targeted for qRT-PCR. Double-stranded cDNA was synthesized from total RNA using the PrimeScript RT reagent kit (Takara). qRT-PCR was performed on a LightCycler 480 real-time PCR machine (Roche) using the SYBR Premix Ex Taq II kit (Takara) in a 10 μl reaction system. All samples were analyzed in triplicate.

3′ RACE PCR was performed to validate the 3′ UTR isoforms of CRFB1, CRFB6, IFNφ1, Baxa, MyD88, and NIK, and the PCR products were sequenced with an ABI Prism 3730 DNA analyzer (Applied Biosystems). Eight 3′ UTR isoforms were obtained and inserted into psiCHECK-2 (Promega), downstream of the Renilla luciferase gene. To isolate the effects of the 3′ UTR on protein expression from efficiency of polyadenylation, both the proximal and distal poly(A) signal sequences were deleted and a synthetic poly(A) was used to terminate each transcript. All ATTTA pentamers in the longer 3′ UTR isoforms were mutated into GCCCG to construct ΔATTTA mutational plasmids. EK293T cells were plated in 48-well plates (at 5 × 10⁴ cells/well) and transfected 24 h later with recombinant plasmids DNA (50 ng/well). Luciferase activities were measured by the Dual-Luciferase reporter assay system (Promega). Luciferase activities were normalized to firefly luciferase activities. Each experiment was performed at least in triplicate and repeated at least twice.

Unless indicated otherwise, all data were plotted as means ± SD. Significant differences between groups were determined by a two-tailed Student t test.

To study whether the global patterns of RNA polymerase II transcription termination is affected during the immune response in zebrafish, the newly developed SAPAS strategy was used to profile poly(A) sites in BC and UC zebrafish spleen. In total, 27,174,773 raw reads were obtained from Solexa sequencing (Illumina Genome Analyzer II). After mapping and filtering as previously described, 18,293,287 reads were directly used to infer the transcript cleavage sites (Table I). Among the filtered reads, 82.6% were mapped to annotated 3′ UTRs, and 4.8% were mapped to 1 kg nucleotide downstream regions. Those remaining (12.6%) were mapped to the coding regions of the annotated protein-coding genes (1.8%), noncoding genes (1.2%), and intergenic regions (9.6%) (Fig. 1A). After clustering the filtered reads as described (3), 39,139 poly(A) sites were identified (Table I), of which 23.6% were mapped to Ensembl transcription termination sites (TTSs; Fig. 1B). The remaining 76.4% were novel. Based on our dataset, 12,724 Ensembl canonical genes were identified by at least one poly(A) site (Table I), representing 51.3% of all canonical genes. A total of 10,032 genes were identified by at least one SAPAS-sampled poly(A) site in the last exon, among which 3933 genes had two or more tandem poly(A) sites and 1588 genes harbored three or more tandem APA sites (Fig. 1C, 1D).

Sequence analysis of poly(A) sites from different genomic locations showed all novel poly(A) site categories except CDS exhibit a similar nucleotide distribution pattern to that of Ensembl TTS, characterized by the canonical poly(A) signal hexamer AAUAAA/AUUAAA and the downstream GU/U-rich element (Fig. 1E, 1F). A 3′ RACE analysis of 30 randomly selected sites, five from each of the six categories, revealed that 27 (90%) were authentic TTSs (Fig. 1G). These results confirmed SAPAS as a powerful and reliable tool for poly(A) sites profiling.

Zebrafish spleen is enriched with immune cells and plays important roles in the immune response (18, 25, 29–31). To test whether the immune system in zebrafish spleen was activated by bacterial stimulation, we compared the read mapping results from the UC and BC and found that ∼6.5% of UC reads and 9.0% of BC reads were mapped to immune response genes. When comparing total genes and immune response genes, the BC library detected more genes than did the UC library, indicating that both global and immune response gene expression profiles are skewed to gene induction rather than gene suppression after bacterial infection (Fig. 2A, 2B). A list of candidate immune response genes was manually collected with reference to KEGG pathways of zebrafish (http://www.genome.jp/kegg/). For pathways that have not been annotated in zebrafish, the corresponding human pathway was referenced. The gene list included 1480 genes from >50 gene families. We tentatively classified these genes as: pattern recognition, effectors, cytokines and signal transducers, and kinases and transcription factors. This classification is inevitably subjective, as some genes can be assigned to different categories, and the functions of many genes were inferred from sequence similarity. Most immune gene families were upregulated, and genes that were categorized as the pattern recognition proteins and effectors showed higher induction than did those that were categorized as cytokines and signal transducers and kinases and transcription factors. Ten of 26 gene families of pattern recognition proteins and effectors were upregulated more than 10-fold, whereas only 3 of 28 gene families of cytokines and signal transducers and kinases and transcription factors were upregulated >10-fold (Fig. 2C).

Having demonstrated that the SAPAS strategy could be effectively used to profile poly(A) sites of zebrafish spleen and that the immune response was globally activated, we performed the test of linear trend alternative to independence to compare the tandem 3′ UTR lengths of the BC and UC libraries as previously described (6, 26). We found 848 genes (false discovery rate of 0.01) with tandem 3′ UTR lengths that showed significant differences between the BC and UC libraries, with 407 genes switched to longer tandem 3′ UTRs and 441 genes switched to shorter tandem 3′ UTRs after bacterial infection (Fig. 3A). However, for genes with switched APA sites, we did not observe a negative correlation between 3′ UTR length and gene expression, and proximal to distal switching and distal to proximal switching of APA sites were similarly represented (Fig. 3A), suggesting that the process of APA during an in vivo immune response is more complex than that observed for in vitro responses (13).

To validate the results generated by the SAPAS method, we performed qRT-PCR on 12 genes that showed significant 3′ UTR length differences in the UC and BC libraries. Ten of 12 genes were confirmed, with failed PCR on two genes (Fig. 3B).

According to the list of candidate immune response genes, 86 immune response genes were identified as APA genes after bacterial infection and, of these, 45 genes switched to expression of longer 3′ UTRs, and 41 genes switched to express shorter 3′ UTRs. Most (60 of 86) of these immune response APA genes were classified as cytokines and signal transducers, whereas four were pattern recognition genes, six were effector genes, and 16 were kinase and transcription factor genes (Fig. 3C, Supplemental Table I). These observations indicate that APA preferentially impacts the intermediate signal transduction cascade as opposed to signal initiation or termination during the immune response. TLR signaling, apoptosis, and JAK-STAT signaling pathways were enriched with APA genes (Fig. 3D, Supplemental Table I).

TLRs play a critical role in the early innate immune response to invading pathogens by recognizing pathogen-associated molecular patterns, and the activation of TLR signaling pathways leads to the release of cytokines that result in innate and adaptive immune responses (32, 33). We found that the genes involved in this pathway were globally upregulated (Fig. 4), and 12 genes in this pathway switched APA sites upon stimulation. Eight of the 12 genes served as signal transducer genes (MyD88, SIGIRR, TAB1, TAB2, TRAF7, CYLD, IκBβ, and IκBε), 3 genes were kinase genes (IKKβ, IKKγ, and MKK4b), and 1 gene was a transcription factor (IRF2) (Fig. 4, Supplemental Table I).

The apoptosis pathway plays an important role in the immune response and during inflammation by facilitating intercellular communication (34). Seventeen genes from the apoptosis pathway switched poly(A) sites in the BC library: three ligand genes (CD40, TNFSF12, and TNFSF13b), one receptor gene (TNFRSF1A), one exonuclease gene (TREX1), one effector caspase gene (CASP 3b), and two antiapoptotic genes (BCL2L1 and BCL2L13) switched to a longer 3′ UTR, whereas one initiator caspase gene (CASP 8L), one exonuclease gene (ENDO-G), three negative regulator genes (DFF45, cFLIP, and BIRC6), two positive regulator genes (BAXa and Baxb), and two kinase genes (NIK and PRKACab) switched to a shorter 3′ UTR (Fig. 4, Supplemental Table I). The tandem 3′ UTR switching of these genes may contribute to the activation of the apoptosis pathway after microbial pathogen stimulation.

The JAK-STAT pathway is a major signaling alternative to the second messenger system, which is activated by IFNs, ILs, growth factors, as well as other chemical messengers, and it regulates the proliferation, differentiation, and apoptosis of cells, as well as the immune response (35, 36). Most of the detected genes involved in JAK-STAT pathway were upregulated, indicating that this pathway was fully activated after bacterial infection (Fig. 4, Supplemental Table I). Thirteen genes from the JAK-STAT pathway switched poly(A) sites after bacterial challenge. One IFN gene (IFNφ1) and four receptor genes (IL-13Ra1, CRFB1, CRFB2, and CRFB7) switched to a longer 3′ UTR, whereas one IL gene (IL16), two receptor genes (CRFB4 and CRFB6), two STAT genes (STAT1a and STAT6), and one negative regulator gene (PIAS4a) switched to a shorter 3′ UTR (Fig. 4, Supplemental Table I). The tandem 3′ UTR switching of these genes may contribute to fine-tuning the regulation of this pivotal multifunctional pathway after pathogenic microbe stimulation.

miRNAs are small noncoding RNAs that bind to the 3′ UTR of mRNAs and function in posttranscriptional regulation of the expression of their target genes (37, 38). We found more miRNA target sites in the extended 3′ UTR than in the common region (Fig. 5A). Mammalian miR-146 family members (miR-146a and miR-146b) are inflammation-inducible miRNAs involved in negative feedback regulation of the TLR signaling pathway (39), which are also induced by bacterial infection in zebrafish in a MyD88- and TRAF6-dependent manner (40). The qRT-PCR assay revealed that miR-146a and miR-146b were also significantly upregulated after bacterial infection (Fig. 5B). Several miR-146 target sites were identified in the 3′ UTR of two cytokine receptor family members (CRFB1 and CRFB6), which are located on the extended 3′ UTR (Fig. 5C). To investigate whether the expression of CRFB1 and CRFB6 is regulated by miR-146 family members, the long and short 3′ UTR isoforms of CRFB1 and CRFB6 were cloned and inserted into luciferase reporter vectors, and these recombinant vectors were transfected into HEK293T cells with zebrafish miR-146a or miR-146b expression vectors. Luciferase assays showed that miR-146a and miR-146b significantly suppressed the expression of the Renilla luciferase gene inserted with the long 3′ UTRs that contains miR-146 target sites, but had no impact on the expression of the Renilla luciferase gene inserted with short 3′ UTRs or long 3′ UTRs with mutated miR-146 target sites (Fig. 5D). These findings indicate that miRNAs may play an important role in regulating the expression of APA genes, especially for immune response APA genes.

AREs are well-characterized cis-acting regulatory sequences located in the 3′ UTR of mRNAs that often consist of one or several AUUUA pentamers in an adenosine- and uridine-rich region (41). It is estimated that ∼5–8% of the human transcriptome contains AREs (27, 42). When searching the ATTTA pentamers in the common and extended 3′ UTRs of APA genes with shortened or lengthened 3′ UTR, we found that the extended 3′ UTRs contained more ATTTA pentamers than did the common 3′ UTRs (Fig. 6A). Additionally, searching of the WWTTATTTATTWW motifs in the common and extended 3′ UTRs of APA genes also showed that the extended 3′ UTRs contained more WWTTATTTATTWW motifs than did the common 3′ UTRs (Fig. 6B). We sought to investigate other cis-regulatory motifs in the switched 3′ UTR and observed five motifs in extended regions that were enriched compared with the common regions. The canonical poly(A) signal AATAAA and its variants were significantly enriched in the extended regions of genes with both shortened or lengthened 3′ UTR, and some U-rich motifs (TTGT, TTTT, and TGTTT) were also enriched in the extended regions of genes with lengthened 3′ UTR, but not in genes with shortened 3′ UTR (Fig. 6C).

A primary function of the AREs is to target specific mRNAs for rapid decay through the interaction with ARE-binding proteins (ARE-BPs) (43, 44). A large number of mRNAs encoding immune-regulatory proteins have AREs in the 3′ UTR, and most mRNAs regulated through ARE-mediated mRNA decay encode cytokines and chemokines, signal transducers that mediate intercellular signaling and coordinate immune responses (41, 45–47). Sixteen ARE-BPs have been identified in mammals (48), and 11 of them were found in zebrafish. Eight ARE-BPs were detected in our UC and BC libraries, all of which were upregulated, except zinc finger protein 36 C3H type-like 1a and 1b (ZFP36L1a and ZFP36L1b) (Fig. 6D, Supplemental Table II). Human AgR (HuR) is one of the well-characterized ARE-BPs, which is implicated in the stabilization of many AU-rich mRNAs, including those that participate in immunity, such as TNFα, IL-3,and IL-8 (49–51). Several tandem 3′ UTR isoforms have been described for both human and mouse HuR mRNAs. Mouse HuR autoregulates its expression by promoting alternative polyadenylation (52). HuR is the most abundant ARE-BP gene in zebrafish spleen and underwent a switch to shorter 3′ UTRs upon pathogenic microbe stimulation (Supplemental Table II). Among ARE-BPs that promote mRNA decay, TTP destabilizes several proinflammatory cytokine mRNAs, including TNFα, GM-CSF, and IL-2 (44, 53). A high mRNA level of this gene was detected in zebrafish spleen (Supplemental Table II).

We used the Web-based version of AREScore (http://arescore.dkfz.de/arescore.pl) to calculate ARE scores in the common and extended 3′ UTRs of immune response APA genes and found the ARE scores of extended 3′ UTRs to be significantly higher than those of the common 3′ UTRs. Of these genes, the extended 3′ UTRs also contained more WWTTATTTATTWW motifs than did the common 3′ UTRs (Fig. 6E). To investigate whether AREs affect the expression of target genes, 3′ UTR isoforms of four immune response APA genes with different ARE scores (Fig. 6F) were cloned and inserted into luciferase reporter vectors, and these recombinant vectors were transfected into HEK293T cells with zebrafish HuR or TTP expression vectors. Luciferase assays showed that zebrafish HuR significantly increased the expression of the Renilla luciferase gene inserted with the 3′ UTRs with greater numbers of ATTTA pentamers and higher ARE scores, that is, the full-length 3′ UTRs of IFN-φ1, Baxa, and MyD88, whereas TTP exerted the opposite effect. Both HuR and TTP had little effect on the expression of the Renilla luciferase gene inserted with 3′ UTRs containing fewer ATTTA pentamers and exhibiting lower ARE scores. Deletion of ATTTA pentamers abolished the effect of HuR and TTP on luciferase activities (Fig. 6G). These findings indicate that AREs and ARE-BPs may play important roles in regulating the expression of APA genes during the immune response, especially of genes with lengthened 3′ UTRs.

Polyadenylation is dynamically and elegantly regulated by a macromolecular complex called the mRNA 3′ processing complex (54, 55). Several trans-acting factors are required for 3′ processing, including the cleavage and polyadenylation specificity factor (CPSF), cleavage stimulation factor (CstF), cleavage factor I (CF Im), cleavage factor II, poly(A) polymerase, symplekin, and RNA polymerase II. CPSF specifically recognizes the A(A/U)UAAA element, at least in part, through CPSF-160. CstF binds to the downstream G/U-rich sequence via CstF-64. Recognition of the poly(A) signal in the absence of the canonical A(A/U)UAAA element primarily depends on the presence of an upstream UGUA sequence motif, which functions in association with CF Im (56, 57). The function of cleavage factor II in mRNA 3′ processing remains poorly understood. The mentioned factors assemble onto mRNA polyadenylation signals to form the 3′ processing complex (Fig. 7A). Regulated processing at the pre-mRNA 3′ ends can be induced by the differential expression of constitutive polyadenylation factors. In most cases, the consequence of this regulatory mechanism is the selection of alternative poly(A) signals that recruit the polyadenylation system in an inefficient manner due to the presence of suboptimal cis-acting elements. In the CPSF-30 mutant of Arabidopsis, poly(A) sites of most genes are altered (58). In our study, CPSF-30 was induced up to 2.9-fold upon bacterial infection, and another subunit, CPSF-73, was induced 3.0-fold after pathogen challenge (Fig. 7B, Supplemental Table II). Tryptophan–aspartic acid (W-D) dipeptide repeat domain 33 (WDR33) is involved in the polyadenylation system by interacting with CPSF complex (54). Zebrafish WDR33 was induced up to 10-fold upon bacterial infection in this study.

Increased expression of CstF-64 in human cell lines has been associated with APA of several pre-mRNAs under various biological conditions (59–62), and knockdown of CstF-64 and its paralog CstF-64τ by RNA interference exerts a high impact on the global APA pattern (63). Although the expression level of CstF-64 was not significantly changed in our study, another subunit of CstF, CstF-50, was induced up to 22.5-fold upon bacterial infection, which may also influence alternative polyadenylation by promoting the assembly of CstF complex (Fig. 7B, Supplemental Table II).

The differential expression of CF Im25 and CF Im68 in mouse male germ cells has been associated with the use of alternative proximal poly(A) signals in a number of transcripts (64), and knockdown of CF Im25 results in alternative poly(A) signal selection (65). CF Im25 and CF Im68 were induced up to 21.5- and 1.5-fold, respectively, upon bacterial infection in our study (Fig. 7B, Supplemental Table II).

A number of RNA binding proteins have been implicated in APA control (reviewed by Di Giammartino et al. in Ref. 66); for example, several splicing factors are also known to influence 3′ processing (reviewed by Millevoi and Vagner in Ref. 55). Polypyrimidine tract–binding protein (PTB), which is a regulator of alternative splicing and internal ribosome entry site–driven translation (reviewed by Valcárcel and Gebauer in Ref. 67), can interfere with CstF binding to the downstream sequence element (68) or can stimulate polyadenylation by increasing the binding of heterogeneous nuclear ribonucleoprotein H (hnRNP H, previously known as a splicing factor) to the upstream sequence element, which in turn recruits CstF and poly(A) polymerase and stimulates cleavage (69, 70). The homologs of PTB in zebrafish, PTB1a and PTB1b, were both upregulated after pathogen stimulation, along with as hnRNP H1. U1 small nuclear ribonucleoprotein A (U1A) plays a key role in inhibiting cleavage and polyadenylation at the secretory poly(A) site of IgM (71), and a direct interaction between U1A and the CPSF-160 increased the efficacy of cleavage and polyadenylation (72). U1A was significantly upregulated following pathogen stimulation in our study (Fig. 7B, Supplemental Table II).

Some polyadenylation-related genes also underwent APA site switching. CPSF-30 switched to a shorter 3′ UTR postinfection, whereas Pcf11, Symplekin, PP1Aa, PP1B, and ELL2 switched to a longer 3′ UTR after stimulation (Supplemental Table II). The APA site switching of these genes may impact their protein output by changing mRNA stability or translation efficiency, thus influencing polyadenylation complex formation.

mRNA isoforms with different 3′ UTRs produced by APA are distinct in their stability, translation efficiency, and cellular location (11). A comprehensive genome-wide assessment of APA genes is required for understanding APA involvement in biological processes. However, traditional quantification methods, such as quantitative PCR and microarray, cannot meet these requirements. By simultaneously measuring the expression and position of poly(A) sites, the SAPAS strategy provides a powerful and reliable method for this kind of transcriptome profiling (6, 26). Using this method, we present the genome-wide profiling of APA sites in an immune organ in vivo. For the genes with switched APA sites identified in our study, we did not observe a clear correlation between the length of the 3′ UTR and gene expression level. Moreover, proximal to distal switching and distal to proximal switching of APA sites were observed in similar numbers, suggesting that the APA functioning during an in vivo immune response may be different from those active in in vitro studies (13), and that they may be regulated by a complex mechanism. The global APA profile of a cell is apparently associated with its proliferation and differentiation status (reviewed by Di Giammartino et al. in Ref. 66). Proliferating cells tend to produce mRNAs with shorter 3′ UTRs by using proximal poly(A) sites, whereas differentiated cells predominately access distal poly(A) sites, resulting in mRNAs with longer 3′ UTRs (7, 13–15, 73, 74). The spleen contains several types of immune cells, including T and B cells, macrophages, neutrophil, dendritic cells, and NK cells. These cells exhibit different sensitivities and accessibilities to pathogens and are activated in a different temporal and spatial manner, forming a complex, codependent regulation network (75). This may contribute to the complexity of immune responses in the spleen via the regulatory network of APA. Previous analyses have studied only genes with limited known APA sites (13), and their observations may not be conclusive. Several recently developed strategies based on next-generation sequencing technology have revealed a large number of previously unidentified poly(A) sites. For example, 69.3% of poly(A) sites identified by the SAPAS strategy in human breast cancer cell lines were novel (26), 56.1% of DRS reads overlapped with previously annotated poly(A) sites (76), and only 12% of poly(A) sites identified in zebrafish embryonic development samples were previously reported (6). In the present study, fewer than a quarter of the identified poly(A) sites were annotated in Ensembl. These newly identified poly(A) sites may dramatically change the analysis of APA and provide a more comprehensive and authentic profile of 3′ UTRs. Such profiles may ultimately help explain how variation in 3′ UTRs is regulated in a given biological process, such as the immune response.

Although no negative correlation was observed between the 3′ UTR length and gene expression, and the proximal to distal switching and distal to proximal switching of APA sites were similarly represented, in contrast to previous reports (13, 14), we found that the TLR signaling, apoptosis, and JAK-STAT signaling pathways, which play central roles in the antibacterial immune response (35, 77), were globally activated and enriched in APA genes, and most of these APA genes were upregulated. This suggests that APA is extensively involved in the expression regulation of genes involved in immune signaling in vivo, which in turn modulates the immune response to bacterial infection.

The loss of miRNA target sites in the shorter 3′ UTRs has been shown to play an important role in oncogenesis and activation of T lymphocytes (13, 14). This mechanism may also contribute to the regulation of APA gene expression in the immune response. A greater number of miRNA target sites were consistently identified in the extended region of 3′ UTR than in the common region. miR-146 family members were induced by bacterial infection and may regulate the immune response synergistically with APA by targeting CRFB1 and CRFB6. ARE was also more enriched in the extended region than in the common region of 3′ UTR. Moreover, ARE-BP genes were generally upregulated upon bacterial infection, indicating that ARE and its binding protein may play an important role in controlling the expression of APA genes, especially genes with lengthened 3′ UTR. These cis-acting elements, miRNA target sites, and ARE may also influence the relative stability of mRNA isoforms, which in turn contribute to poly(A) site switching.

The differential expression of core polyadenylation factors exerts a global impact on the selection of alternative poly(A) sites (10, 66). CstF-64 is a well-characterized APA regulator of IgM and NFATc in lymphocyte activation in a splice-dependent manner (59, 60, 78), but the expression of CstF-64 was slightly downregulated upon bacterial infection, whereas CstF-50, CF Im25, and WDR33 were induced >10-fold in our study, suggesting that different strategies can be adopted to modulate 3′ processing. Some transcription and splicing factor genes involved in APA were also differentially expressed in the immune response, and, notably, some genes switched to a longer or shorter 3′ UTR postinfection. The variation in expression level and poly(A) site switching of these genes may contribute to APA during the immune response. APA-mediated regulation of mRNA transcripts with varying 3′ UTRs may be a new paradigm for regulation of immune response, as it is for development (6).

The authors have no financial conflicts of interest.