Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Widespread Shortening of 3’ Untranslated Regions and Increased Exon Inclusion Are Evolutionarily Conserved Features of Innate Immune Responses to Infection

Abstract

The contribution of pre-mRNA processing mechanisms to the regulation of immune responses remains poorly studied despite emerging examples of their role as regulators of immune defenses. We sought to investigate the role of mRNA processing in the cellular responses of human macrophages to live bacterial infections. Here, we used mRNA sequencing to quantify gene expression and isoform abundances in primary macrophages from 60 individuals, before and after infection with Listeria monocytogenes and Salmonella typhimurium. In response to both bacteria we identified thousands of genes that significantly change isoform usage in response to infection, characterized by an overall increase in isoform diversity after infection. In response to both bacteria, we found global shifts towards (i) the inclusion of cassette exons and (ii) shorter 3’ UTRs, with near-universal shifts towards usage of more upstream polyadenylation sites. Using complementary data collected in non-human primates, we show that these features are evolutionarily conserved among primates. Following infection, we identify candidate RNA processing factors whose expression is associated with individual-specific variation in isoform abundance. Finally, by profiling microRNA levels, we show that 3’ UTRs with reduced abundance after infection are significantly enriched for target sites for particular miRNAs. These results suggest that the pervasive usage of shorter 3’ UTRs is a mechanism for particular genes to evade repression by immune-activated miRNAs. Collectively, our results suggest that dynamic changes in RNA processing may play key roles in the regulation of innate immune responses.

Author Summary

Changes in gene regulation have long been known to play important roles in both innate and adaptive immune responses. While transcriptional responses to infection have been well-characterized, much less is known about the extent to which co-transcriptional mechanisms of mRNA processing are involved in the regulation of immune defenses. In this study, we sought to investigate the role of mRNA processing in the cellular responses of human macrophages to live bacterial infection. Using primary human macrophages derived from whole blood samples from 60 individuals, we sequenced mRNA both before and after infection with two live bacteria. We show that immune responses to infection are accompanied by pervasive changes in mRNA isoform usage, with systematic shifts towards increased cassette exon inclusion and shortening of Tandem 3’ UTRs post-infection. These patterns are conserved in nonhuman primates, supporting their functional importance across evolutionary time. Complementary microRNA profiling revealed that shortened 3’ UTRs are enriched for target sites of macrophage-expressed miRNAs, many of which are specifically activated after infection to regulate the innate immune response. Our results therefore provide the first genome-wide empirical support for the idea that actively regulated shifts towards shorter 3’ UTRs might allow specific genes to evade repression by immune-activated miRNAs.

Citation: Pai AA, Baharian G, Pagé Sabourin A, Brinkworth JF, Nédélec Y, Foley JW, et al. (2016) Widespread Shortening of 3’ Untranslated Regions and Increased Exon Inclusion Are Evolutionarily Conserved Features of Innate Immune Responses to Infection. PLoS Genet 12(9): e1006338. https://doi.org/10.1371/journal.pgen.1006338

Editor: Christine A. Wells, The University of Queensland, AUSTRALIA

Received: June 24, 2016; Accepted: September 2, 2016; Published: September 30, 2016

Copyright: © 2016 Pai et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: Data generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE73765, which comprises of the mRNA-seq data (GSE73502) and small RNA-seq data (GSE73478).

Funding: This work was supported by grants from the Canadian Institutes of Health Research (232519), the Human Frontiers Science Program (CDA-00025/2012), a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada (RGPIN/435917-2013), and the Canada Research Chairs Program (950-228993) (to LBB). Computational resources at MIT supported by the National Science Foundation (Grant No. 0821391). AAP. was supported by a Jane Coffin Childs postdoctoral fellowship. GB was supported by a postdoctoral fellowship from the Fonds de la Recherche en Santé du Québec. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Innate immune responses depend on robust and coordinated gene expression programs involving the transcriptional regulation of thousands of genes [1–3]. These regulatory cascades start with the detection of microbial-associated products by pattern recognition receptors, which include Toll Like Receptors (TLRs), NOD-like receptors and specific C-type lectins [4,5]. These initial steps are followed by the activation of key transcription factors (e.g., NF-κB and interferon regulatory factors) that orchestrate the inflammatory and/or antiviral response signals involved in pathogen clearance and the subsequent development of appropriate adaptive immune responses [4].

Although much attention has been devoted to characterizing transcriptional changes in response to infectious agents or other immune stimuli, we still know remarkably little about the contribution of co-transcriptional changes–specifically, changes in alternative pre-mRNA processing–to the regulation of the immune system [6–8]. Alternative splicing can impact cellular function by creating distinct mRNA transcripts from the same gene, which may encode unique proteins that have distinct or opposing functions. For instance, in human and mouse cells, several alternatively spliced forms of genes involved in the TLR pathway have been shown to function as negative regulators of TLR signaling in order to prevent uncontrolled inflammation [9–11]. Additionally, several in-depth studies of individual genes suggest that alternative splicing plays an important role in increasing the diversity of transcripts encoding MHC molecules [12] and in modulating intracellular signaling and intercellular interactions through the expression of various isoforms of key cytokines [13,14], cytokine receptors [9,15,16], kinases [17] and adaptor proteins [18–20]. Yet, the extent to which changes in isoform usage are a hallmark of immune responses to infection remains largely unexplored at a genome-wide level [8,21]. Moreover, we still do not know which RNA processing mechanisms contribute most to the regulation of immune responses, or how these disparate mechanisms are coordinated.

To address these questions, we investigated genome-wide changes in transcriptome patterns after independent infections with Listeria monocytogenes or Salmonella typhimurium in primary human macrophages. Because of the distinct molecular composition of these two pathogens and the way they interact with host cells, they activate distinct innate immune pathways after infection [22,23]. Thus, our study design allowed us to evaluate the extent to which changes in isoform usage are pathogen-specific or more generally observed in response to bacterial infection. Furthermore, the large number of individuals in our study allowed us to use natural variation in RNA processing to gain insight into fine-tuned inter-individual regulation of immune responses. Our results provide a comprehensive picture of the role that RNA processing might play in regulating early innate immune responses to bacterial infection in human antigen-presenting cells.

Results

Infection with either Listeria or Salmonella induces dramatic changes in mRNA expression levels. Following 2 hours of infection with each bacteria, we collected RNA-seq data from 60 matched non-infected and infected samples, with an average of 30 million reads sequenced per sample (see Methods; S1 Table). The first principal component of the complete gene expression data set clearly separated infected from non-infected samples, and both PC1 and PC2 clustered infected samples by pathogen (i.e. Listeria or Salmonella; Fig 1A). Accordingly, we observed a large number of differences in gene expression levels between infected and non-infected cells, with 5,809 (39%) and 7,618 (51%) of genes showing evidence for differential gene expression (DGE) after infection with Listeria and Salmonella, respectively (using DESeq2, FDR ≤ 0.1% and |log2(fold change)| > 0.5; S1A Fig, S2 Table). As expected, the sets of genes that responded to either infection were strongly enriched (FDR ≤ 1.8×10−6) for genes involved in immune-related biological processes such as the regulation of cytokine production, inflammatory responses, or T-cell activation (S3 Table).

thumbnail
Download:
Fig 1. Gene expression and isoform proportion differences in response to bacterial infection.

(A) Principal component analysis of gene expression data from all samples (PC1 and PC2 on the x- and y-axis, respectively). (B) IL24, a gene with significant changes in isoform usage upon infection with Listeria and Salmonella. For each IL24 isoform in response to infection, plotted are the average relative isoform usages across samples (left panel, top) with their isoform structures (left panel, bottom) and corresponding fold changes (right panel; log2 scale; with standard error bars). Isoforms are ordered by relative abundance in non-infected samples, and colored circles (right panel) correspond to colors in barplot (left panel). (C) Number of genes showing only DIU, only DGE, and both DIU and DGE upon infection with Listeria and Salmonella, (11,353 genes tested).

https://doi.org/10.1371/journal.pgen.1006338.g001

In order to uncover changes in isoform usage in response to infection using our RNA-sequencing data, we applied two complementary approaches for analyzing RNA processing. We first measured abundances of full isoforms to gain a holistic perspective of the transcriptome landscape. Second, we used the relative abundances of individual exons within a gene to gain a finer understanding of the specific RNA processing mechanisms likely to be involved in the regulation of immune responses.

Pervasive differential isoform usage in response to infection

We initially sought to assess changes in isoform usage after infection using proportional abundances of the different isoforms encoded by the same gene. To do so, we used the transcripts per million (TPM) values reported by the RSEM software [24] and calculated relative proportions by dividing isoform-level expression by the overall expression level of the gene (i.e., TPM summed across all possible isoforms). We then quantified differential isoform usage (DIU) between conditions for genes with at least two annotated isoforms (N = 11,353) using a multivariate generalization of the Welch's t-test that allowed us to test whether, as a group, the proportional abundances of the different isoforms in each gene were significantly different between infected and non-infected cells (Methods). We found that 1,456 (13%; after Listeria infection) and 2,862 (25%; after Salmonella infection) genes showed evidence for significant DIU after infection (FDR ≤ 1%; Fig 1B for an example of an DIU gene, S1B Fig, S4 Table). Across all DIU genes, we observed a mean of 10% (± 8% s.d.) and 13% (± 10% s.d.) change in relative isoform usage upon infection with Listeria and Salmonella, respectively (defined as the maximum change upon infection in relative isoform usage across transcripts of each gene, S1C Fig). Following infection with Listeria and Salmonella, 6.2% and 10.2%, respectively, of genes with a dominant isoform prior to infection (see methods) switch to using an alternative dominant isoform (S1E Fig).

DIU genes were significantly enriched for genes involved in immune responses (FDR ≤ 0.01; S3 Table), including several cytokines (e.g. IL1B, IL7, IL20, and IL24 in 2B), chemokines (CCL15), regulators of inflammatory signals (e.g. ADORA3, MAP3K14) and genes encoding co-stimulatory molecules required for T cell activation and survival (e.g. CD28) (S2 Fig). Notably, 86% of genes with DIU upon infection with Listeria were also classified as having DIU after infection with Salmonella, suggesting that changes in isoform usage are likely to be common across a wide variety of immune triggers. To assess the robustness of our findings, we tested for DIU using isoform-specific expression levels calculated with Kallisto, an alignment-free quantification method [25] (see Supplementary Methods in S1 Text). At an FDR of 5%, we observed a 74% (Listeria) and 78% (Salmonella) agreement with DIU genes identified using RSEM, confirming that identification of DIU genes is largely robust to the method used for isoform quantification.

Considering the large proportion of genes with significant DIU upon infection, we sought to understand whether these changes arose from a shift in the usage of a dominant isoform or increased isoform diversity within a given gene. To do so, we calculated the Shannon diversity index, which measures the evenness of the isoform usage distribution for each gene before and after infection (low values reflect usage of one or few isoforms; high values reflect either the usage of a more diverse set of isoforms or more equal representation of the same set of isoforms). ΔShannon thus quantifies the change in isoform diversity after infection (S5 Table). The majority of genes (63% in Listeria and 68% in Salmonella) show an increase in isoform diversity (ΔShannon > 0) after infection (S1D Fig), indicating a shift away from usage of the primary (pre-infection) isoform after infection and a general increase in the diversity of transcript species following infection.

Previously, it had been reported that DGE and DIU generally act independently to shape the transcriptomes of different mammalian tissues [26]. In contrast, we found that DIU genes were significantly enriched for genes with DGE (both up- and down-regulated genes, Fisher's exact test, P ≤ 7.5x10-6 for both Listeria and Salmonella) (S3A Fig). Despite such enrichment, a substantial proportion of genes (47% (690 genes) in Listeria infected samples and 39% (1,105 genes) in Salmonella infected samples) with significant changes in isoform usage after infection do not exhibit DGE following infection (1C). Even after relaxing the FDR threshold for DGE by 100-fold, 685 (Listeria) and 1096 (Salmonella) DIU genes still exhibited no significant DGE. Thus, a considerable fraction of transcriptome changes in response to infection do occur solely at the level of RNA processing, and independently of changes in mean expression levels.

Directed shifts in 3’ UTR length and alternative splicing following infection

Though useful in documenting the substantial transcriptome changes upon infection, aggregate isoform levels calculated by RSEM do not distinguish between different types of RNA processing changes. Thus, we next examined changes in five distinct categories of alternative RNA processing, allowing us to zoom in on particular molecular mechanisms that underlie differential isoform usage in response to infection. Specifically, we used human transcript annotations from Ensembl to quantify usage of (1) alternative first exons (AFEs), (2) alternative last exons (ALEs), (3) alternative polyadenylation sites, leading to tandem 3’ untranslated regions (TandemUTRs), (4) retained introns (RIs), and (5) skipped exons (SEs). AFEs, ALEs, and TandemUTRs correspond to the alternative usage of terminal exons primarily affecting UTR composition, whereas RIs and SEs correspond to internal splicing events that usually affect the open reading frame.

For each gene or exon within each RNA processing category, we used the MISO software [27] to calculate a "percent spliced in" (PSI or Ψ) value (S6 Table), defined as the proportion of transcripts from a gene that contain the “inclusion” isoform (defined as the longer isoform for RIs, SEs, and TandemUTRs, or use of the exon most distal to the gene for AFEs and ALEs). ΔΨ values thus represent the difference between PSI values calculated for the infected versus non-infected samples. Overall, we observed many significant changes in RNA processing (N ≥ 1,098) across all categories in both bacteria (significance was defined as Bayes Factor > 5 in at least 10% of individuals and |mean ΔΨ| > 0.05), compared to null expectations derived from measuring changes among pairs of non-infected samples (N = 29) (2A, S7 Table). Our criteria for determining significant changes also allows us to choose exons that are more likely to be consistently changing across many individuals; significantly changing exons have lower variance in ΔΨ values across individuals than exons that are not changing after infection (S4 Fig). The greatest proportion of changes after infection occurred among retained intron and TandemUTR events. PSI values for 7.3% and 14% of retained introns and 7.6% and 16.7% of TandemUTRs significantly changed after infection with Listeria or Salmonella, respectively. When we considered the set of genes associated with at least one significant change in RNA processing, we observed an over-representation of Gene Ontology categories involved in immune cell processes and the cellular response to a stimulus (Fig 2B, S8 Table).

thumbnail
Download:
Fig 2. RNA processing changes in response to bacterial infection.

(A) Proportion of events for RNA processing category that are significantly changing after infection with Listeria (left), Salmonella (middle), or variation between non-infected samples as a control (right). Numbers indicate the number of significant changes per category. (B) Significantly Gene Ontology categories for genes with any significant RNA processing change (FDR ≤ 10%). (C) Distribution of ΔΨ values for each RNA processing category. Negative values represent less inclusion, while positive values represent more inclusion, as defined by the schematic exon representations.

https://doi.org/10.1371/journal.pgen.1006338.g002

Genome-wide shifts towards more prevalent usage of inclusion or exclusion isoforms have previously been observed as signatures of cellular stress responses or developmental processes [28–33]. To examine directional shifts in particular RNA processing mechanisms, we used mean ΔΨ values (defined as the mean ΔΨ across all individuals for each gene or exon) and observed a striking global shift towards the usage of upstream polyadenylation sites in tandem 3’ UTR regions, indicating pervasive 3’ UTR shortening following infection with either Listeria or Salmonella (Fig 2C; both P < 2.2×10−16 with Mann-Whitney U test; 94.4% and 98% of significant TandemUTRs changes in Listeria and Salmonella). Although less frequent, we also see a general trend towards usage of upstream polyadenylation sites in alternative last exons. Interestingly, we found a striking correlation in the extent to which a given individual has a global shift towards 3’ UTR shortening and towards usage of an upstream ALE (Pearson R = 0.83, P = 4.97 × 10−16 for Listeria and Pearson R = 0.78, P = 2.65 × 10−13 for Salmonella, S10A Fig). These correlations suggest that the changes in ALEs and TandemUTRs are likely being regulated by similar or coordinated mechanisms. Our findings are consistent with previous reports that suggested that ALEs and TandemUTRs could both be regulated by changes in the activities of cleavage and polyadenylation factors, rather than splicing machinery [34–36]. Additionally, we found a substantial increase in the inclusion of skipped exons after infection (both P < 2.2×10−16 with Mann-Whitney U test; 79.8% and 77.3% of significant SE changes in Listeria and Salmonella, S5 Fig). The consistency of these patterns suggests that a dedicated post-transcriptional regulatory program might underlie genome-wide shifts in RNA processing after bacterial infection.

Quantifying changes at the terminal ends of transcripts is challenging due to known biases in standard RNA-seq protocols [37]. Thus, we validated our observation of 3’ UTR shortening upon infection by sequencing RNA from 6 individuals after infection with Listeria and Salmonella with a 3’RNA-seq protocol, which specifically captures the 3’ most end of transcripts (see Methods). TPM expression values calculated from 3’ RNA-seq data highly correlate with TPMs from RNA-seq data, establishing that this method is able to quantitatively measure of the number of transcripts in a cell (R > 0.82 across all conditions, P < 2.2 × 10−16, S6 Fig). Meta-gene plots around the polyadenylation sites of core and extended regions of TandemUTRs show a notable increase after infection in sequencing coverage around the upstream polyadenylation site relative to the downstream polyadenylation site, specifically for genes with significantly shorter 3’ UTRs upon infection as assessed by the RNA-seq data (Fig 3A) Finally, to look at 3’ UTR shortening in specific genes, we calculate Ψ values for TandemUTRs using the 3’RNA-seq data (see Supplementary Methods in S1 Text) and find again that genes that are significantly changing 3’ UTR usage in the RNA-seq data exhibit a significant shift towards negative ΔΨ values in the 3’RNA-seq data (P < 2.2 × 10−16; Fig 3B).

thumbnail
Download:
Fig 3. 3’RNA sequencing shows increased usage of upstream polyadenylation sites upon infection.

(A) Meta-gene distributions of 3’RNA-seq read densities at the upstream polyA sites (core regions, left) and downstream polyA sites (extended regions, right) of Tandem 3’ UTRs after infection with Listeria or Salmonella (top and bottom, respectively). Shown are the read distributions for non-infected samples across all Tandem 3’ UTRs (black) and infected samples at Tandem 3’ UTRs that significantly change after infection (yellow) or show no change after infection (blue), as called by the RNA-seq data. (B) Distribution of ΔΨ values calculated from 3’RNA-seq data for Tandem 3’ UTRs. We observe significant shifts (P < 2.2 × 10−16 for both Listeria and Salmonella) towards negative ΔΨ values in Tandem UTRs that are identified as significantly changing in RNA-seq data (yellow) relative to Tandem UTRs without any change after infection (blue).

https://doi.org/10.1371/journal.pgen.1006338.g003

RNA processing changes persist 24 hours after infection

To evaluate whether the strong global shifts we observed were the result of a sustained immune response to the bacteria rather than a short-lived cellular stress response, we sequenced RNA from a subset of our individuals (N = 6) after 24 hours (h) of infection with Listeria or Salmonella. Samples sequenced after 24h were paired with matched non-infected controls cultured for 24h to calculate ΔΨ values for alternative RNA processing categories using MISO [27]. Across all categories, we found many more significant changes in RNA processing and larger ΔΨ values after infection for 24h than after 2h (S7 Fig). These results further support a key role for RNA processing in controlling innate immune responses to infection. For genes or exons that changed significantly after 2h of infection in this set of 6 individuals (see Methods), the majority of changes (67%) remained detectable at 24h, usually with similar magnitude and direction (Fig 4A). The one notable exception is that introns that are preferentially retained at 2h in infected cells tended to be spliced out at 24h. Our results thus indicate that the splicing of alternative retained introns during the course of the immune response to infection is temporally dynamic.

thumbnail
Download:
Fig 4. Directed shifts in RNA processing persist across time, stimulus conditions, and closely related species.

(A) Correlations between ΔΨ values after 2 hours of infection (x-axis) and ΔΨ values after 24 hours of infection (y-axis), presented as density plots per RNA processing category in contrasting colors. Plotted are events that are significant after 2 hours of either Listeria or Salmonella infection. (B) Distributions of ΔΨ values for skipped exons (purple) and TandemUTRs (blue) following infection of whole blood cells with lipopolysaccharide, assessed in human (top), chimpanzee (middle) and rhesus macaque (bottom) individuals (N = 6 per species). We observe prominent global shifts in isoform distributions in human and macaque (SEhuman P = 0.003, TandemUTRhuman P = 2.4×10−13; SEmacaque P = 0.05, TandemUTRmacaque P = 8.8×10−6) with a more modest trend observed in chimpanzee, potentially due to poor transcript annotations in the chimpanzee genome (SEchimpanzee P = 0.26, TandemUTRchimpanzee P = 0.002). All p-values were calculated using a Student’s t-test testing deviation from a mean of zero.

https://doi.org/10.1371/journal.pgen.1006338.g004

Shortening of 3’ UTR regions and increased exon inclusion are evolutionarily conserved responses to infection in primates

Since the innate immune response is ancient, one might expect core aspects of this pathway to be evolutionarily conserved among primates. To explore this question, we used high-depth RNA-seq data from human, chimpanzee, and rhesus macaque primary whole blood cells (N = 6 individuals from each species) before and after 4 hours of stimulation with lipopolysaccharide (LPS), the major component of the outer membrane of Gram-negative bacteria. Using an analysis pipeline similar to our analysis of infection after 24 hours of infection with Listeria or Salmonella, we calculated ΔΨ values for alternative RNA processing categories using MISO [27] (see Supplementary Methods in S1 Text, S9 Table). As with macrophages, we observed consistent strong shifts towards increased exon inclusion and shorter 3’ UTRs in human whole blood cells stimulated with LPS. Furthermore, we observed similar trends toward increased exon inclusion and significant 3' UTR shortening in macaques and chimpanzees, supporting a conserved response of RNA processing pathways during innate immune responses in primates (Fig 4B). In addition, these data show that the shortening of 3’ UTRs and increased exon inclusion are not limited to isolated macrophages but also observed as part of the immune response engaged by interacting primary granulocytes and lymphocytes, more closely mimicking in vivo cellular responses to infection.

Increased exon inclusion is associated with increased gene expression levels Since we observed strong connections between overall isoform usage and differential gene expression, we explored the relationship between the occurrence of particular RNA processing changes and the direction of gene expression changes in response to infection. We found that genes with significant differences in skipped exon usage pre- and post-infection (up to 80% of which have increased skipped exon inclusion after infection) were more likely to be up-regulated after infection (T-test, both P ≤ 1×10−4 for Listeria and Salmonella; Fig 5A; S8A Fig). This association remained when gene expression levels were calculated using reads mapping to constitutive exons only (S8A Fig), eliminating concerns about alternative spliced isoforms impacting gene expression estimates. Our observation is consistent with a recent study that reported increased gene expression in tissues with increased inclusion of evolutionarily novel SEs [38]. Genes with other types of splicing changes showed no trend in expression changes.

thumbnail
Download:
Fig 5. Relationship between RNA processing changes and gene expression changes.

(A) Distribution of fold changes in gene expression (y-axis, log2 scale) for genes with significant skipped exon changes after infection (purple) and genes with no change after infection (gray). (B) Distribution of Spearman correlations between ΔΨ and fold change in gene expression across 60 individuals per gene (solid line), for genes with only one annotated alternative event and significantly changed SE usage (purple; nL = 46 genes and nS = 97 genes) or tandem UTR usage (blue; nL = 36 genes and nS = 86 genes). Dotted lines show distribution of the correlation coefficients after permuting ΔΨ values.

https://doi.org/10.1371/journal.pgen.1006338.g005

Genes with significant skipped exon usage were enriched for Gene Ontology categories such as “response to stimulus,” and also for “RNA processing” and “RNA stability” categories (S9A Fig). Included within these latter categories are several members of two major splicing factor families, SR proteins and hnRNPs, which showed increased inclusion of skipped exons and increased gene expression after infection (S9B Fig). The SR protein and hnRNP families are known to function in complex auto- and cross-regulatory networks [39,40], often antagonizing one another’s effects on splicing [41–43]. Consistent with the slightly higher up-regulation of members of the hnRNP family, we observed a consistent enrichment for intronic splicing silencers in upstream intronic regions around significantly excluded skipped exons (S9C Fig). These elements are often bound by hnRNPs to promote exclusion of a cassette exon [42])

Variation across individuals provides insight into putative trans-regulatory factors

A unique feature of our study design is our sampling of 60 individuals after infection with both Listeria and Salmonella. By taking advantage of natural variation in the regulation of RNA processing to infection, we aimed to gain further insight into the connections between RNA processing and changes in gene expression levels in response to infection. For each gene, we calculated the correlation between the inter-individual variation in RNA processing changes (ΔΨ) and the fold changes in gene expression levels upon infection. For most categories, we observed shifts in the distribution of correlations between RNA processing and gene expression levels relative to permuted controls (S8B Fig), with skipped exons and TandemUTRs showing the most consistent patterns (Kolmogorov-Smirnov test, P ≤ 0.005; Fig 5B). Increased skipped exon inclusion correlates with increased gene expression of the gene across individuals in 69% of the genes with significant skipped exons. In addition, preferential expression of shorter alternative 3’ UTR isoforms tends to be correlated with increased up-regulation of the associated genes in response to infection. Our findings thus suggest that RNA processing changes may directly impact gene expression levels, or at least are commonly regulated with changes in transcription or mRNA stability.

Substantial inter-individual variation in the genome-wide average RNA processing patterns also correlated strongly with average gene expression changes upon infection, particularly for changes in 3’ UTR usage. That is, individuals who exhibited more skipped exon inclusion or 3’ UTR shortening also tended to exhibit larger global shifts towards up-regulation of gene expression levels upon infection (S10B Fig).

Next, we sought to exploit our relatively large sample size to identify candidate trans-factors impacting the amount of skipped exon and TandemUTR changes observed across individuals. Specifically, we examined the relationship between mean changes in exon inclusion and fold change in expression for all expressed genes, across the 60 individuals. When we focus on proteins with known RNA binding functions or domains–hypothesizing that these factors are more likely to regulate post-transcriptional mechanisms–we find many significant correlations (FDR ≤ 1%) for known regulators of the corresponding RNA processing category (S11 Fig). For instance, genes whose changes in expression correlated with individual-specific mean 3’ UTR shortening included several factors previously implicated in regulating polyadenylation site usage (such as cleavage and polyadenylation factors [44], hnRNP H proteins [27], and hnRNP F [45]; S11A Fig). Additionally, many known regulators of alternative splicing–including a number of SR proteins and hnRNPs–have increased expression levels in individuals with larger extents of skipped exon inclusion (S11B Fig).

3’ UTR shortening as a means to evade repression by immune-related miRNAs The largest global shift observed in alternative isoform abundance was for 3’ UTR shortening, where as many as 98% of significantly changing TandemUTRs shifted towards usage of an upstream polyA site post-infection. Previously, global 3’ UTR shortening has been associated with proliferating cells, particularly in the context of cellular transitions in development [28], cell differentiation [46], cancerous [29] states, or global 3’ UTR lengthening in embryonic development [33]. However, macrophages do not proliferate, which we confirmed with a BrdU labeling assay in both resting and infected conditions (S12 Fig). This result indicates that the pervasive 3’ UTR shortening observed in response to infection is due to an active cellular process that is independent of cell division and might be mechanistically distinct from that leading to 3’ UTR shortening in proliferating cells.

Previous studies in proliferating cells postulated that 3’ UTR shortening can act as a way to evade regulation by microRNAs (miRNAs), since crucial miRNA target sites are most often found in 3’ UTR regions [28]. To evaluate this hypothesis we performed small-RNA sequencing in 6 individuals after 2h and 24h of infection with Listeria and Salmonella. When focusing specifically on miRNAs expressed in macrophages (S10 Table), we found that the extended UTR regions of infection-sensitive 3’ UTRs have a significantly higher density of miRNA target sites compared to TandemUTRs that do not change in response to infection (Fig 6A). This increased density of miRNA target sites was restricted to the extension region, which is subject to shortening after infection: no difference in the density of target sites was observed in the common “core” 3’ UTR regions of the same genes (S13 Fig).

thumbnail
Download:
Fig 6. Tandem 3’ UTR shortening allows evasion of regulation by miRNAs.

(A) Distribution of frequency of miRNA target sites per nucleotide in the extended regions of Tandem UTRs that either show no change after infection (grey) or significantly change after infection (blue). (B) Significantly enriched miRNA target sites (FDR ≤ 10%, |FC| > 1.5) in the extended regions of significantly changing Tandem UTRs after infection with Listeria-only (top), with Salmonella-only (middle), or both bacteria (bottom). For each bacteria, the barplots in the left panels show the fold enrichment (x-axis, log2 scale) of target sites in the extended regions. White bars represent non-significant enrichments. Panels on the right show the fold change in miRNA expression (x-axis, log2 scale with standard error bars) after either 2 hours of infection (light colors) or 24 hours of infection (dark colors).

https://doi.org/10.1371/journal.pgen.1006338.g006

Next, we tested if the increased density of miRNA target sites was driven by the enrichment of target sites for particular miRNAs expressed in macrophages. For each miRNA that was expressed in either non-infected or infected macrophages, we calculated an enrichment score assessing whether the target sites of that miRNA were significantly enriched in the extended region of significantly shortened 3’ UTRs (with a background distribution matched for sequence composition; S14 Fig; see Supplemental Methods in S1 Text). We found 15 miRNAs with significantly enriched target sites (FDR ≤ 10% and enrichment ≥ 1.5 fold) in either the Listeria or Salmonella conditions, with 6 miRNAs over-represented in both bacterial conditions (Fig 6B). Interestingly, of the miRNAs with the highest enrichment after infection with either Listeria or Salmonella, two (miR-146b and miR-125a) have previously been shown to be important regulators of the innate immune response [47–51]. Accordingly, we found that 4 of the miRNAs with the highest enrichment after infection (miR-125a, miR146b, miR-3661, and miR-151b) are all up-regulated following infection, and strongly so after 24h of infection (Fig 6B). In the majority of cases, greater 3’ UTR shortening was associated with a stronger increase in gene expression across individuals, suggesting that shifts towards shorter 3’ UTRs after infection may allow these transcripts to escape from repression by specific immune-induced miRNAs.

Discussion

Taken together, our results provide strong evidence that RNA processing may play a key role in the regulation of innate immune responses to infection. Despite known differences in the regulatory pathways elicited in macrophages by Listeria and Salmonella, our data revealed striking similarities in the overall patterns of RNA processing induced in response to both pathogens. Given the strong overlaps, we chose throughout this study to focus on consistent patterns across the two bacteria.

Though genes that have differential isoform usage in response to infection are more likely to show changes in gene expression, we found that as many as 47% of genes showing differential isoform usage have no evidence for differences in gene expression upon infection. This observation suggests that a considerable proportion of genes are regulated predominantly by post-transcriptional mechanisms, with minimal changes in transcriptional regulation. Differential isoform usage is coupled to a prominent increase in the diversity of isoforms following infection. Intriguingly, genes showing the strongest increases in diversity upon infection were strongly enriched among down-regulated genes (S3B Fig). This observation raises the possibility that the increased isoform diversity could reflect a shift toward less stable isoforms, e.g., those targeted by nonsense-mediated decay (NMD), thus reducing transcript levels [52,53].

Upon infection of human macrophages with either Listeria or Salmonella, we observed substantial shifts towards increased skipped exon inclusion and shortening of 3’ UTRs. These RNA processing patterns remained prominent when interrogated in human whole blood after an alternative immune stimulus (LPS), despite the fact the number of significant RNA processing changes in whole blood was generally smaller than that observed in purified macrophages infected with live bacteria. This reduced effect is likely explained by the fact that in whole blood, we are interrogating a very heterogeneous cell population within which only a fraction (~3–6%) of the cell types are actually responding to LPS stimulation, such as monocytes/macrophages [54]. Additionally, transcriptional changes induced by live bacteria such as Listeria or Salmonella, which concomitantly activates multiple innate immune pathways [22,23], are larger than transcriptional changes that can induced by a single purified ligand, such as LPS, which activates primarily the TLR4 pathway [55]. Despite this reduced statistical power, we still observed conserved RNA processing signatures of immune stimulation through the primate lineage, supporting the likely functional importance of these cellular responses and the ubiquity of RNA processing changes after immune stimulation across evolution.

The most prominent hallmark of RNA processing changes after infection is a global shortening of 3’ UTRs. Interestingly, this effect resembles the global 3’ UTR shortening seen in proliferating activated T cells and transformed cells [28] [29] [46]. Given that both T-cell and macrophage activation result in the activation of shared signaling cascades (e.g., c-Jun N-terminal kinase (JNK)–a key regulator of apoptosis–and nuclear factor κ-B (NF-κB) activation), it is possible that, at least to some extent, 3’ UTR shortening observed in T-cells and macrophages stems from the activation of these shared pathways. Previous studies hypothesized that 3’ UTR shortening allows mRNAs to escape regulation by dynamically changing miRNAs, suggesting that macrophages may co-opt generally used mechanisms for the same purpose. [28,29]. This shift towards shorter 3’ UTRs could result from increased miRNA-dependent degradation of transcripts with longer 3’ UTRs post-infection. However, genes that shift towards shorter 3’ UTRs were more highly expressed post-infection, making this explanation unlikely. Furthermore, the increase in overall miRNA target site density is reminiscent of a previous observation that genes with regulated 3’ UTR usage across tissues have more target sites for ubiquitously expressed miRNAs in their distal 3’ UTRs [56]. However, in contrast to results across tissues, we observe that the target sites in our distal 3’ UTR regions are directly related to miRNAs whose expression specifically changes after infection.

Specifically, our data support previous hypotheses and find evidence for that 3’ UTR shortening in macrophages is associated with the loss of target sites for a subset of immune-regulated miRNAs, including miR-146b, miR-125a, and miR-151b (Fig 6B). These miRNAs target the extended 3’ UTR regions of several critical activators of the innate immune response–including the IRF5 transcription factor [57] and the MAP kinases MAPKAP1 and MAP2K4 [58]–all of which shift to usage of shorter 3’ UTRs after infection. It is tempting to speculate that without such 3’ UTR shortening, the binding of one or more of miRNAs to the 3’ UTR region of these genes would compromise their up-regulation and the subsequent activation of downstream immune defense pathways. Given that many miRNAs are involved in either (or both) mRNA degradation and translational efficiency, we cannot know for sure which of these processes is prevented by the 3’ UTR shortening of these genes.

While our results point to the importance of 3’ UTR shortening to allow genes to evade post-transcriptional regulation by miRNAs, this is not generalizable to all miRNAs induced in response to infection. Indeed, some of the miRNAs most highly induced after infection, such as miR-155 and miR-29b, are not enriched for binding sites in the longer 3’UTR regions. Thus, while miRNA regulation is crucial to mammalian immune response [59–62], the mechanisms that dictate what specific 3’UTR are shortened or not will require further investigation. Indeed, it is likely that miRNA targets whose 3’ UTRs are shortened in response to infection to evade miRNA regulation are actually key targets of the same miRNAs in other cellular conditions.

The convergence towards similar isoform outcomes across many disparate genes suggests the activation of factors acting in trans to drive global shifts towards inclusion of cassette exons or usage of upstream proximal polyadenyalation sites. Taking advantage of our relatively large sample size, we were able to identify candidate trans-factors whose up-regulation upon infection might have widespread impacts on RNA processing patterns, including some known regulators of alternative splicing or 3’end processing (S11 Fig). Sets of RNA binding factors likely act in combination to influence final transcriptome states [63]. For instance, a general up-regulation of splicing machinery following infection [64], allowing the cell to recognize and use additional splice sites, might contribute to increased exon inclusion and greater isoform diversity observed following infection. Members of both the hnRNP and SR protein splicing factor families show up-regulation of their overall gene expression levels (S9B Fig), likely contributing to specific instances of exon inclusion as well as exclusion. Correspondingly, changes in the regulation of many canonical cleavage and polyadenylation factors were tied to individual-specific shifts in 3’ UTR shortening. The significant correlation between global shifts in ALE usage and TandemUTRs across individuals further support the notion that cleavage and poladenylation pathways are regulating the changing 3’ UTR landscape upon infection (S10A Fig). Intriguingly, fold changes in hnRNP H1 gene expression have high correlations with both skipped exon inclusion and 3’ UTR shortening, consistent with its known splicing function and suggesting that this factor may play multiple roles in shaping the transcriptome in response to bacterial infection.

An alternative mechanism that has been proposed for the global regulation of 3’ UTR shortening in T-cell activation and neuronal differentiation involves a phenomenon called telescripting. This is a phenomenon of long-distance RNA polymerase II transcription due to protection of nascent transcripts from premature cleavage and polyadenylation by binding of the U1 snRNP. In activated or differentiating cells, moderately lower U1 snRNA levels relative to global increases in nascent mRNA production restrict this secondary function of U1 [65]. While we did not directly measure mRNA levels, we observe a moderate but significant increase in total RNA concentration after 2 hours of infection that is consistent with the increase in nascent RNA observed by Berg et al. during neuronal differentiation [65] (across 60 samples, paired t-test; P = 0.006 for Listeria and P = 4.57 × 10−6 for Salmonella), combined with no changes in the levels of U1 snRNA (relative to 5S rRNA, Supplemental Methods in S1 Text, S15 Fig). Thus it is possible that a delayed response of U1 snRNA could promote the usage of upstream polyadenylation sites in both TandemUTRs and ALEs. Interestingly, Berg et al. also observe increased up-regulation of individual internal exons upon moderate U1 knockdown in neurons [65], which is consistent with our observation of increased skipped exon inclusion.

Substantial inter-individual variation in the genome-wide average RNA processing patterns also correlated strongly with average gene expression changes upon infection, particularly for changes in 3’ UTRs and skipped exons. That is, individuals who exhibited more skipped exon inclusion or usage of upstream polyadenylation sites also tended to exhibit larger global shifts towards up-regulation of gene expression levels upon infection (S10B Fig). Though these patterns might result from stochastic variation in levels of one or several trans-factors, they are also likely to be influenced by genetic differences among individuals. It will be interesting to ask whether global isoform differences might potentially be reflective of variation in an individual’s susceptibility to infection.

Materials and Methods

Complete details of the experimental and statistical procedures can be found in SI Materials and Methods. Briefly, blood samples from 60 healthy donors were obtained from Indiana Blood Center with informed consent and ethics approval from the Research Ethics Board at the CHU Sainte-Justine (protocol #4022). All individuals recruited in this study were healthy males between the ages of 18 and 55 y old. Blood mononuclear cells from each donor were isolated by Ficoll-Paque centrifugation and blood monocytes were purified from peripheral blood mononuclear cells (PBMCs) by positive selection with magnetic CD14 MicroBeads (Miltenyi Biotec). In order to derive macrophages, monocytes were then cultured for 7 days in RPMI-1640 (Fisher) supplemented with 10% heat-inactivated FBS (FBS premium, US origin, Wisent), L-glutamine (Fisher) and M-CSF (20ng/mL; R&D systems). After validating the differentiation/activation status of the monocyte-derived macrophages we infected them at a multiplicity of infection (MOI) of 10:1 for Salmonella typhimurium and an MOI of 5:1 for Listeria monocytogenes for 2- and 24-hours. For primate comparisons, whole blood was drawn from 6 healthy adult humans (CHU Sainte-Justine), 6 common chimpanzees, (Texas Biomedical Research Institute), and 6 rhesus macaques (Yerkes Regional Primate Center). Human samples were acquired with informed consent and ethics approval from the Research Ethics Board (CHU Sainte-Justine, #3557). Non-human primate samples were acquired in accordance with individual institutional IACUCC requirements. For all species, 1ml of blood was drawn into a media-containing tube (Trueculture tube, Myriad, US) spiked with ultrapure LPS (Invitrogen, USA) or endotoxin-free water (control). Samples were stimulated with 1ug/ml LPS, at 37C for 4 hours before total blood leukocytes were isolated and RNA collected. Genome-wide gene expression profiles of untreated and infected/treated samples were obtained by RNA-sequencing for both mRNA transcripts and small RNAs. After a series of quality checks (SI Materials and Methods), mRNA transcript abundances were estimated using RSEM [66] and Kallisto [25] and microRNA expression estimates were obtained as previously described[67]. To detect genes with differential isoform usage between two groups of non-infected and infected samples we used a multivariate generalization of the Welch's t-test. Shannon entropy Hsh (also known as Shannon index) was applied to measure the diversity of isoforms for each target gene before and after infection. Changes across individual RNA processing events were quantified using the MISO software package (v0.4.9) [27] using default settings and hg19 version 1 annotations. Events were considered to be significantly altered post-infection if at least 10% (n > = 6) of individuals had a BF > = 5 and the |mean ΔΨ| ≥ 0.05 (S6 Table). To conduct a targeting sequencing of 3’ ends of transcripts, We used Smart-3SEQ (Foley et al, manuscript in preparation), a modified version of the 3SEQ method [68]. We used a custom gene ontology script to test for enrichment of functional annotations among genes that significantly changed isoform usage in response to infection (SI Materials and Methods). All predicted miRNA target sites within annotated TandemUTR regions were obtained using TargetScan (v6.2) [69].

Ethics statement

Institutional ethics (IRB) approval was obtained by the CHU Sainte-Justine Institutional Review Board (approved project #3557) and all subjects gave written informed consent prior to participation.

Data access

Data generated in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.hlm.nih.gov/geo/) under the accession number GSE73765, which comprises of the mRNA-seq data (GSE73502) and small RNA-seq data (GSE73478).

Supporting Information

S1 Text. Supplementary Methods and Materials.

https://doi.org/10.1371/journal.pgen.1006338.s001

(PDF)

S1 Fig. Differential isoform usage after infection.

(A) Volcano plots of differential expressionafter infection with Listeria and Salmonella in the left and right panels, respectively.–log10 p-values (y-axis) testing for differential expression are plotted against average log2 fold changes in expression levels (x-axis) for genes that are not differentially expressed (black) and genes that with significant differential expression after infection (FDR ≤ 0.1% and |log2(fold change)| ≥ 0.5; blue). (B) Distribution of p-values for the differential isoform usage test upon infection with Listeria and Salmonella. Expected distribution of p-values under the null hypothesis of no significant difference between the mean relative isoform usage is shown in grey. (C) A comparison of DIU effect sizes between DIU genes at 1% FDR (following infection with Listeria in dark pink and Salmonella in gold) and the background of all genes (lighter colors). Effect sizes are defined as maximum change in relative isoform abundances per gene upon infection. (D) Distributions of ΔShannon diversity index (ΔSDI) after infection with Listeria and Salmonella. Null distribution (black dotted line) was generated by permuting samples across conditions. (E) The percentage of genes with a dominant isoform before infection (left) and the fraction of these genes where the dominant isoform changes after infection (right).

https://doi.org/10.1371/journal.pgen.1006338.s002

(PDF)

S2 Fig. Representative examples of immune-related genes with significant DIU after infection.

Heatmaps below each example represent the variation in isoform usage across the 60 individuals, where each vertical bar represents one individual. The dominant isoform in non-infected samples is represented in grey, while the predominant isoform in infected samples is in a color. Darker bars represent increased relative usage, while lighter bars represent decreased relative usage.

https://doi.org/10.1371/journal.pgen.1006338.s003

(PDF)

S3 Fig. Gene expression and differential isoform diversity after infection.

(A) Proportion of genes with differential isoform usage (y-axis) among all genes (grey), differentially expressed genes that are down-regulated after infection (orange), and differentially expressed genes that are up-regulated after infection (green). (B) log2 fold changes (with standard error bars) for odds of low ΔShannon and high ΔShannon among the set of all the differentially expressed genes, up-regulated differentially expressed genes, and down-regulated differentially expressed genes, following infection with Listeria (dark pink) and Salmonella (gold).

https://doi.org/10.1371/journal.pgen.1006338.s004

(PDF)

S4 Fig. Variance in RNA processing across 60 individuals.

Distribution of coefficient of variation (x-axis) in ΔΨ values for each isoform in a given RNA processing category.

https://doi.org/10.1371/journal.pgen.1006338.s005

(PDF)

S5 Fig. Two representative examples of genes with significant skipped exon changes after infection.

Sashimi plots for two genes, PTK2B and STK40, that have significant changes in skipped exon usage both after 2 hours and 24 hours of infection with either bacteria.

https://doi.org/10.1371/journal.pgen.1006338.s006

(PDF)

S6 Fig. Relationship between 3'RNA-seq and RNA-seq data.

Correlations between TPMs from 3' RNA-seq data (y-axis) and TPMs from RNA-seq data (x-axis).

https://doi.org/10.1371/journal.pgen.1006338.s007

(PDF)

S7 Fig. RNA processing changes 24hr after infection.

(A) Proportion of significantly changing events (x-axis) after 2 hours of infection (dark colors) and 24 hours of infection (light colors) with either Listeria (top) or Salmonella (bottom). Numbers indicate the significant events at corresponding timepoints using only the 6 individuals used for these cross-timepoint analyses. (B) Distribution of ΔΨ values (x-axis) for significantly changing events in each event type after either 2 hours (top) or 24 hours (bottom) of infection.

https://doi.org/10.1371/journal.pgen.1006338.s008

(PDF)

S8 Fig. Relationships between RNA processing and gene expression changes after infection.

(A) Distributions of overall fold changes in gene expression (y-axis, log2 scale) for genes that have significant splicing changes in each event type (colored boxplots) and genes with no splicing changes after infection (grey). Gene expression values are calculated using either full transcript models (top) or only constitutively included exons (bottom). (B) Distribution of spearman correlations between the ΔΨ of an event and the fold change in gene expression, across individuals, for events that are not significantly changing after infection (left) and significantly changing events (right). Solid lines represent the observed values and dotted lines represent a distribution of correlations after permuting the correspondence between ΔΨs and fold changes in gene expression.

https://doi.org/10.1371/journal.pgen.1006338.s009

(PDF)

S9 Fig. Characteristics of significantly changing skipped exon events.

(A) Significantly enriched gene ontology categories for genes with significant skipped exon changes after infection. Color indicates the fold enrichment of the number of observed genes relative to the number of genes expected to be in that category. (B) Fold changes in gene expression (y-axis, log2 scale) after infection for 2 families of splicing factors (hnRNPs and SR proteins) relative to a background distribution of fold changes in all genes (C) Fold enrichments of splicing regulatory elements (SREs) in exonic regions and surrounding intronic regions (±100bp). SREs assessed included exonic splicing enhancers (green in exons), exonic splicing silencer (red in exons), intronic splicing enhancers (green in introns), and intronic splicing silencers (green in introns). Enrichments were calculated separately for significantly included skipped exons (top) and significantly excluded skipped exons (bottom). Shading of the bars indicates the significance of the enrichment.

https://doi.org/10.1371/journal.pgen.1006338.s010

(PDF)

S10 Fig. Inter-individual variation of global shifts in alternative splicing.

(A) Correlations between the mean ΔΨ values per individual for ALEs (x-axis) and TandemUTRs (y-axis). (B) Spearman correlations between the mean ΔΨ value per individual and the mean fold change of gene expression per individual for corresponding genes.

https://doi.org/10.1371/journal.pgen.1006338.s011

(PDF)

S11 Fig. Identifying RNA binding factors that might be involved in trans-regulation of RNA processing after infection.

For both Tandem UTRs (A) and Skipped exons (B), the top panel is a scatter plot of Spearman correlations between the individual-specific mean ΔΨ values and individual-specific fold change of gene expression values for all expressed genes (grey) in both Listeria (x-axis) and Salmonella (y-axis). Genes with significant correlations (FDR ≤ 1%) in both Listeria and Salmonella conditions are plotted in black, and those factors with known RNA-binding properties are colored by their functional category. The bottom panels show distributions of the average Spearman correlation values for each of the RNA-binding functional categories with significant correlations.

https://doi.org/10.1371/journal.pgen.1006338.s012

(PDF)

S12 Fig. Cellular proliferation after bacterial infection.

BrdU cell proliferation assay in macrophages (top panel) and LCLs (bottom panel) in non-infected cells and in cells infected with Listeria or Salmonella for both 2 and 24 hours. BrdU incorporates into newly synthesized DNA and therefore the quantity of BrdU incorporated into cells (x-axis) is a direct indication of cell proliferation. No evidence for cellular proliferation was observed in macrophages, in contrast to the high rates of proliferating cells observed in our positive control population of LCLs.

https://doi.org/10.1371/journal.pgen.1006338.s013

(PDF)

S13 Fig. Density of miRNA target sites in 3’ UTR regions.

Barplots in grey indicate Tandem 3’ UTRs that are not changing after infection, while colored barplots indicate Tandem 3’ UTRs that are significantly changing after either Listeria (pink) or Salmonella (yellow) infections.

https://doi.org/10.1371/journal.pgen.1006338.s014

(PDF)

S14 Fig. Nucleotide composition of Tandem 3’ UTR regions.

The distribution of GC content in core (brown) and extended (blue) regions of Tandem 3' UTRs that are either not changing or significantly changing after infection (left panel). While there 3' UTRs that are significantly changing generally have greater overall GC content, this is true for both the core and extended regions of the 3’ UTRs. Thus, the distributions of relative GC content when comparing the regions is distributed around 1 for both significantly changing (yellow) and not changing (blue) Tandem 3' UTRs (right panel).

https://doi.org/10.1371/journal.pgen.1006338.s015

(PDF)

S15 Fig. Investigating effects of telescripting.

Distributions of fold changes in 5s rRNA (grey, n = 20), U1 (green, n = 20), and total RNA (yellow, n = 60) concentrations after infection. 5s rRNA concentrations were calculated by taking the ΔCT value from qPCR measurements across 20 samples. Relative U1 snRNA concentrations were calculated by taking the ΔΔCT values from qPCR measurements across 20 samples, where U1 snRNA is measured relative to 5s rRNA concentrations in the same samples. Total RNA concentrations were estimated from Nanodrop measurements of RNA extraction yields across all 60 samples in our study. All samples (both non-infected and infected) were plated at exactly the same macrophage cellular density at the start of the experiment and we have confirmed that these macrophages do not proliferate either before or after infection (S12 Fig). There are no significant shifts in the distribution of 5s rRNA or U1 snRNA fold changes after infection (t-test for null of μ = 1, both P > 0.01 for 5s rRNA and both P > 0.2 for U1 snRNA), while the distribution of total RNA fold changes are significantly increased after infection (t-test for null of μ = 1, P < 0.001 for both Listeria and Salmonella samples).

https://doi.org/10.1371/journal.pgen.1006338.s016

(PDF)

S1 Table. mRNA-seq statistics.

Details of RNA-seq data for samples after infection at 2 hours or 24 hours. Libraries were made as described in Supplementary Section 1.2.1 in S1 Text.

https://doi.org/10.1371/journal.pgen.1006338.s017

(XLSX)

S2 Table. Gene Level Differential Gene Expression.

Estimates of differential gene expression for the 14,954 genes tested by DESeq2, as described in Supplementary Section 1.3.2 in S1 Text. Included are the log2 fold change estimates, p-values for differential expression, and Benjamini-Hochberg adjusted p-values for differential expression between non-infected samples and both Listeria and Salmonella samples.

https://doi.org/10.1371/journal.pgen.1006338.s018

(XLSX)

S3 Table. Gene Ontology Enrichments for Differential Expression.

Enriched Gene Ontology categories for genes that are differentially expressed, or show differential isoform usage after infection. Gene Ontology analyses were performed using GSEA, as described in Supplementary Section 1.3.8 in S1 Text.

https://doi.org/10.1371/journal.pgen.1006338.s019

(XLSX)

S4 Table. Differential Isoform Usage.

Estimates of differential isoform usage for the 11,353 genes tested as described in Supplementary Section 1.3.3 in S1 Text. Included are the changes in relative isoform usage for the isoform with the greatest relative change after infection, the p-value and Benjamini-Hochberg FDR for DIU using isoform expression levels calculated using RSEM, and the FDR for DIU using isoform expression levels calculated using Kallisto.

https://doi.org/10.1371/journal.pgen.1006338.s020

(XLSX)

S5 Table. Change in Isoform Diversity after Infection.

Estimates of Shannon diversity indicies for 11,353 genes both before and after infection with both Listeria and Salmonella, as described in Supplementary Section 1.3.4 in S1 Text.

https://doi.org/10.1371/journal.pgen.1006338.s021

(XLSX)

S6 Table. Differential isoform usage across 5 specific RNA processing categories.

Estimates of Ψ and ΔΨ values calculated by MISO (Supplementary Section1.3.5 in S1 Text) for each individual for each for the 5 different RNA processing categories assess in this study: alternative first exons (AFE), alternative last exons (ALE), retained introns (RI), skipped exons (SE), and Tandem 3’ UTRs (TandemUTR). Also included are Bayes Factors testing for differential isoform usage per event for each individual, mean ΔΨ values across all individuals and individuals with significant changes, and a column indicating whether an event was called significant based on the criteria defined in Supplementary Section 1.3.5 in S1 Text.

https://doi.org/10.1371/journal.pgen.1006338.s022

(XLSX)

S7 Table. Summaries of Significant RNA processing categories.

A summary of the proportion of annotated events that were detected and determined to be significant (according to the criteria defined in Supplementary Section 1.3.5 in S1 Text) for each condition and RNA processing category. Included are tables for metrics after both 2 hours of infection and 24 hours of infection.

https://doi.org/10.1371/journal.pgen.1006338.s023

(XLSX)

S8 Table. Gene Ontology Enrichments for significant RNA processing changes.

Top gene ontology categories for genes that have significant RNA processing changes after infection with either Listeria or Salmonella. All tests were done using an iterative gene ontology analysis as described in Supplementary Section 1.3.5 in S1 Text and gene names refer to UniProt gene accessions.

https://doi.org/10.1371/journal.pgen.1006338.s024

(XLSX)

S9 Table. Summary of Significant RNA processing categories following infection with LPS, across primate species.

A summary of the proportion of annotated evenst that were detected and determined to be significant (according to the criteria defined in Supplementary Section 1.3.5 in S1 Text) for each species and RNA processing category following infection of primary whole blood cells with LPS.

https://doi.org/10.1371/journal.pgen.1006338.s025

(XLSX)

S10 Table. miRNA mapping and differential expression statistics.

Details of short RNA-seq data for 6 samples sequenced after both 2 hours and 24 hours of infection with either Listeria or Salmonella, as described in Supplementary Sections 1.2.3 and 1.3.10 in S1 Text. Also included are the log2 fold change estimates and p-values for differential expression of miRNAs after both 2 hours and 24 hours of infection, tested by DESeq2 as described in Supplementary Section 1.3.10 in S1 Text.

https://doi.org/10.1371/journal.pgen.1006338.s026

(XLS)

Acknowledgments

We thank P Freese for a categorized list of RNA binding proteins, G McVicker for an initial script to perform iterative gene ontology analyses and P Sudmant for help parsing primate transcript annotations. We thank Y Katz, JM Taliaferro, P Sudmant, J Tung, and members of the Barreiro and Burge labs for helpful discussions and comments on the manuscript. We thank Calcul Québec and Compute Canada for managing and providing access to the supercomputer Briaree from the University of Montreal, used to do many computations reported in the manuscript.

Author Contributions

  1. Conceived and designed the experiments: AAP LBB.
  2. Performed the experiments: APS JFB AD VY.
  3. Analyzed the data: AAP GB YN JCG.
  4. Contributed reagents/materials/analysis tools: JWF KJS ZPJ REL CBB.
  5. Wrote the paper: AAP GB LBB.

References

  1. 1. Huang Q, Liu D, Majewski P, Schulte LC, Korn JM, Young RA, et al. The Plasticity of Dendritic Cell Responses to Pathogens and Their Components. Science. 2001;294: 870–875. pmid:11679675
  2. 2. Smale ST. Selective Transcription in Response to an Inflammatory Stimulus. Cell. 2010;140: 833–844. http://dx.doi.org/10.1016/j.cell.2010.01.037 pmid:20303874
  3. 3. Medzhitov R, Horng T. Transcriptional control of the inflammatory response. Nature Reviews Immunology. Nature Publishing Group; 2009;9: 692–703. pmid:19859064
  4. 4. Medzhitov R, Janeway CA Jr. Innate immune recognition and control of adaptive immune responses. Seminars in Immunology. 1998;10: 351–353. http://dx.doi.org/10.1006/smim.1998.0136 pmid:9799709
  5. 5. Kawai T, Akira S. The role of pattern-recognition receptors in innate immunity: update on Toll-like receptors. Nat Immunol. Nature Publishing Group; 2010;11: 373–384. pmid:20404851
  6. 6. Lynch KW. Consequences of regulated pre-mRNA splicing in the immune system. Nature Reviews Immunology. 2004;4: 931–940. pmid:15573128
  7. 7. Martinez NM, Lynch KW. Control of alternative splicing in immune responses: many regulators, many predictions, much still to learn. Immunological Reviews. Wiley Online Library; 2013;253: 216–236. pmid:23550649
  8. 8. Alasoo K, Estrada FM, Hale C, Gordon S, Powrie F. Transcriptional profiling of macrophages derived from monocytes and iPS cells identifies a conserved response to LPS and novel alternative transcription. bioRxiv. 2015.
    • 9. Rao N, Nguyen S, Ngo K, Fung-Leung W-P. A novel splice variant of interleukin-1 receptor (IL-1R)-associated kinase 1 plays a negative regulatory role in Toll/IL-1R-induced inflammatory signaling. Molecular and Cellular Biology. Am Soc Microbiol; 2005;25: 6521–6532. pmid:16024789
    • 10. Wells CA, Chalk AM, Forrest A, Taylor D, Waddell N, Schroder K, et al. Alternate transcription of the Toll-like receptor signaling cascade. Genome Biology. BioMed Central Ltd; 2006;7: R10.
    • 11. O’Connor BP, Danhorn T, De Arras L, Flatley BR, Marcus RA, Farias-Hesson E, et al. Regulation of Toll-like Receptor Signaling by the SF3a mRNA Splicing Complex. Lin X, editor. PLoS Genet. Public Library of Science; 2015;11: e1004932. pmid:25658809
    • 12. Ishitani A, Geraghty DE. Alternative splicing of HLA-G transcripts yields proteins with primary structures resembling both class I and class II antigens. Proceedings of the National Academy of Sciences. National Acad Sciences; 1992;89: 3947–3951.
    • 13. Bihl MP, Heinimann K, Rudiger JJ, Eickelberg O, Perruchoud AP, Tamm M, et al. Identification of a novel IL-6 isoform binding to the endogenous IL-6 receptor. American journal of respiratory cell and molecular biology. Am Thoracic Soc; 2002;27: 48–56. pmid:12091245
    • 14. Nishimura H, Yajima T, Naiki Y, Tsunobuchi H, Umemura M, Itano K, et al. Differential roles of interleukin 15 mRNA isoforms generated by alternative splicing in immune responses in vivo. The Journal of experimental medicine. Rockefeller Univ Press; 2000;191: 157–170. pmid:10620614
    • 15. Goodwin RG, Friend D, Ziegler SF, Jerzy R, Falk BA, Gimpel S, et al. Cloning of the human and murine interleukin-7 receptors: demonstration of a soluble form and homology to a new receptor superfamily. Cell. Elsevier; 1990;60: 941–951. pmid:2317865
    • 16. Koskinen LL, Einarsdottir E, Dukes E, Heap GA, Dubois P, Korponay-Szabo IR, et al. Association study of the IL18RAP locus in three European populations with coeliac disease. Hum Mol Genet. Oxford Univ Press; 2009;18: 1148–1155. pmid:19103669
    • 17. Jensen LE, Whitehead AS. IRAK1b, a novel alternative splice variant of interleukin-1 receptor-associated kinase (IRAK), mediates interleukin-1 signaling and has prolonged stability. Journal of Biological Chemistry. ASBMB; 2001;276: 29037–29044. pmid:11397809
    • 18. Gray P, Michelsen KS, Sirois CM, Lowe E, Shimada K, Crother TR, et al. Identification of a novel human MD-2 splice variant that negatively regulates Lipopolysaccharide-induced TLR4 signaling. The Journal of Immunology. Am Assoc Immnol; 2010;184: 6359–6366. pmid:20435923
    • 19. Ohta S, Bahrun U, Tanaka M, Kimoto M. Identification of a novel isoform of MD-2 that downregulates lipopolysaccharide signaling. Biochem Biophys Res Commun. Elsevier; 2004;323: 1103–1108. pmid:15381113
    • 20. Janssens S, Burns K, Tschopp J, Beyaert R. Regulation of interleukin-1-and lipopolysaccharide-induced NF-κB activation by alternative splicing of MyD88. Current Biology. Elsevier; 2002;12: 467–471. pmid:11909531
    • 21. Rodrigues R, Grosso AR, Moita L. Genome-wide analysis of alternative splicing during dendritic cell response to a bacterial challenge. PLoS One. 2013.
    • 22. Haraga A, Ohlson MB, Miller SI. Salmonellae interplay with host cells. Nature Reviews Microbiology. Nature Publishing Group; 2008;6: 53–66. pmid:18026123
    • 23. Pamer EG. Immune responses to Listeria monocytogenes. Nature Reviews Immunology. Nature Publishing Group; 2004;4: 812–823. pmid:15459672
    • 24. Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN. RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics. 2010;26: 493–500. pmid:20022975
    • 25. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology. Nature Publishing Group; 2016.
    • 26. Barbosa-Morais NL, Irimia M, Pan Q, Xiong HY, Gueroussov S, Lee LJ, et al. The evolutionary landscape of alternative splicing in vertebrate species. Science. American Association for the Advancement of Science; 2012;338: 1587–1593. pmid:23258890
    • 27. Katz Y, Wang ET, Airoldi EM, Burge CB. Analysis and design of RNA sequencing experiments for identifying isoform regulation. Nat Meth. 2010;7: 1009–1015.
    • 28. Sandberg R, Neilson JR, Sarma A, Sharp PA, Burge CB. Proliferating Cells Express mRNAs with Shortened 3' Untranslated Regions and Fewer MicroRNA Target Sites. Science. 2008;320: 1643–1647. pmid:18566288
    • 29. Mayr C, Bartel DP. Widespread Shortening of 3´UTRs by Alternative Cleavage and Polyadenylation Activates Oncogenes in Cancer Cells. Cell. 2009;138: 673–684. pmid:19703394
    • 30. Wong JJL, Ritchie W, Ebner OA, Selbach M, Wong JWH, Huang Y, et al. Orchestrated Intron Retention Regulates Normal Granulocyte Differentiation. Cell. Elsevier Inc; 2013;154: 583–595. pmid:23911323
    • 31. Lackford B, Yao C, Charles GM, Weng L, Zheng X, Choi EA, et al. Fip1 regulates mRNA alternative polyadenylation to promote stem cell self-renewal. The EMBO Journal. 2014;33: 878–889. pmid:24596251
    • 32. Shalgi R, Hurt JA, Lindquist S, Burge CB. Widespread inhibition of posttranscriptional splicing shapes the cellular transcriptome following heat shock. Cell Reports. Elsevier; 2014;7: 1362–1370. pmid:24857664
    • 33. Ji Z, Lee JY, Pan Z, Jiang B, Bin Tian. Progressive lengthening of 3' untranslated regions of mRNAs by alternative polyadenylation during mouse embryonic development. Valcarcel J, editor. PNAS. Proceedings of the …; 2009;106: 7028–7033. pmid:19372383
    • 34. Takagaki Y, Manley JL. Levels of Polyadenylation Factor CstF-64 Control IgM Heavy Chain mRNA Accumulation and Other Events Associated with B Cell Differentiation. Molecular Cell. 1998;2: 761–771. pmid:9885564
    • 35. Di Giammartino DC, Nishida K, Manley JL. Mechanisms and Consequences of Alternative Polyadenylation. Molecular Cell. 2011;43: 853–866. pmid:21925375
    • 36. Taliaferro JM, Vidaki M, Oliveira R, Olson S, Zhan L, Saxena T, et al. Distal Alternative Last Exons Localize mRNAs to Neural Projections. Molecular Cell. 2016;61: 821–833. pmid:26907613
    • 37. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10: 57–63. pmid:19015660
    • 38. Merkin JJ, Chen P, Alexis MS, Hautaniemi SK. Origins and impacts of new mammalian exons. Cell Reports. 2015;10: 1992–2005. pmid:25801031
    • 39. Huelga SC, Vu AQ, Arnold JD, Liang TY, Liu PP, Yan BY, et al. Integrative Genome-wide Analysis Reveals Cooperative Regulation of Alternative Splicing by hnRNP Proteins. Cell Reports. 2012;1: 167–178. pmid:22574288
    • 40. Pandit S, Zhou Y, Shiue L, Coutinho-Mansfield G, Li H, Qiu J, et al. Genome-wide analysis reveals SR protein cooperation and competition in regulated splicing. Molecular Cell. Elsevier; 2013;50: 223–235.
    • 41. House AE, Lynch KW. Regulation of Alternative Splicing: More than Just the ABCs. Journal of Biological Chemistry. 2008;283: 1217–1221. pmid:18024429
    • 42. Wang Z, Burge CB. Splicing regulation: From a parts list of regulatory elements to an integrated splicing code. RNA. 2008;14: 802–813. pmid:18369186
    • 43. Busch A, Hertel KJ. Evolution of SR protein and hnRNP splicing regulatory factors. WIREs RNA. Wiley Online Library; 2012;3: 1–12. pmid:21898828
    • 44. Zheng D, Tian B. RNA-binding proteins in regulation of alternative cleavage and polyadenylation. Systems Biology of RNA Binding Proteins. Springer; 2014. pp. 97–127.
    • 45. Veraldi KL, Arhin GK, Martincic K, Chung-Ganster LH, Wilusz J, Milcarek C. hnRNP F Influences Binding of a 64-Kilodalton Subunit of Cleavage Stimulation Factor to mRNA Precursors in Mouse B Cells. Molecular and Cellular Biology. 2001;21: 1228–1238. pmid:11158309
    • 46. Lackford B, Yao C, Charles GM, Weng L, Zheng X, Choi EA, et al. Fip1 regulates mRNA alternative polyadenylation to promote stem cell self-renewal. The EMBO Journal. 2014;33: 878–889. pmid:24596251
    • 47. Taganov KD, Boldin MP, Chang KJ. NF-κB-dependent induction of microRNA miR-146, an inhibitor targeted to signaling proteins of innate immune responses. 2006.
      • 48. Perry MM, Williams AE, Tsitsiou E, Larner-Svensson HM, Lindsay MA. Divergent intracellular pathways regulate interleukin-1β-induced miR-146a and miR-146b expression and chemokine release in human alveolar epithelial cells. FEBS Letters. Federation of European Biochemical Societies; 2009;583: 3349–3355.
      • 49. Boldin MP, Baltimore D. MicroRNAs, new effectors and regulators of NF-κB. Immunological Reviews. 2012;246: 205–220. pmid:22435557
      • 50. Kim SW, Ramasamy K, Bouamar H. MicroRNAs miR-125a and miR-125b constitutively activate the NF-κB pathway by targeting the tumor necrosis factor alpha-induced protein 3 (TNFAIP3, A20). 2012. pmid:22550173
      • 51. Banerjee S, Cui H, Xie N, Tan Z, Yang S, Icyuz M, et al. miR-125a-5p Regulates Differential Activation of Macrophages and Inflammation. Journal of Biological Chemistry. 2013;288: 35428–35436. pmid:24151079
      • 52. Lewis BP, Green RE, Brenner SE. Evidence for the widespread coupling of alternative splicing and nonsense-mediated mRNA decay in humans. PNAS. 2003;100: 189–192. pmid:12502788
      • 53. Green RE, Lewis BP, Hillman RT, Blanchette M, Lareau LF, Garnett AT, et al. Widespread predicted nonsense-mediated mRNA decay of alternatively-spliced transcripts of human normal and disease genes. Bioinformatics. Oxford Univ Press; 2003;19: i118–i121. pmid:12855447
      • 54. Zarember KA, Godowski PJ. Tissue Expression of Human Toll-Like Receptors and Differential Regulation of Toll-Like Receptor mRNAs in Leukocytes in Response to Microbes, Their Products, and Cytokines. The Journal of Immunology. American Association of Immunologists; 2002;168: 554–561. pmid:11777946
      • 55. Poltorak A, He X, Smirnova I, Liu M-Y, Van Huffel C, Du X, et al. Defective LPS Signaling in C3H/HeJ and C57BL/10ScCr Mice: Mutations in Tlr4 Gene. Science. American Association for the Advancement of Science; 1998;282: 2085–2088. pmid:9851930
      • 56. Lianoglou S, Garg V, Yang JL, Leslie CS, Mayr C. Ubiquitously transcribed genes use alternative polyadenylation to achieve tissue-specific expression. Genes & Development. Cold Spring Harbor Lab; 2013;27: 2380–2396.
      • 57. Lawrence T, Natoli G. Transcriptional regulation of macrophage polarization: enabling diversity with identity. Nature Reviews Immunology. Nature Publishing Group; 2011;11: 750–761. pmid:22025054
      • 58. Arthur JSC, Ley SC. Mitogen-activated protein kinases in innate immunity. Nature Reviews Immunology. Nature Publishing Group; 2013;13: 679–692. pmid:23954936
      • 59. Xiao C, Rajewsky K. MicroRNA Control in the Immune System: Basic Principles. Cell. 2009;136: 26–36. pmid:19135886
      • 60. Ma F, Xu S, Liu X, Zhang Q, Xu X, Liu M, et al. The microRNA miR-29 controls innate and adaptive immune responses to intracellular bacterial infection by targeting interferon-[gamma]. Nat Immunol. Nature Research; 2011;12: 861–869.
      • 61. Vigorito E, Kohlhaas S, Lu D, Leyland R. miR‐155: an ancient regulator of the immune system. Immunological Reviews. 2013;253: 146–157. pmid:23550644
      • 62. Luo X, Ranade K, Talker R, Jallal B, Shen N, Yao Y. microRNA-mediated regulation of innate immune response in rheumatic diseases. Arthritis Research & Therapy 2013 15:2. BioMed Central; 2013;15: 210.
      • 63. Brooks AN, Duff MO, May G, Yang L, Bolisetty M, Landolin J, et al. Regulation of alternative splicing in Drosophila by 56 RNA binding proteins. Genome Research. Cold Spring Harbor Lab; 2015;: gr.192518.115.
      • 64. Munding EM, Shiue L, Katzman S, Donohue JP, Ares M Jr. Competition between Pre-mRNAs for the Splicing MachineryDrives Global Regulation of Splicing. Molecular Cell. Elsevier Inc; 2013;51: 338–348.
      • 65. Berg MG, Singh LN, Younis I, Liu Q, Pinto AM, Kaida D, et al. U1 snRNP Determines mRNA Length and Regulates Isoform Expression. Cell. 2012;150: 53–64. pmid:22770214
      • 66. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. BioMed Central Ltd; 2011. p. 323. pmid:21816040
      • 67. Siddle KJ, Deschamps M, Tailleux L. A genomic portrait of the genetic architecture and regulatory impact of microRNA expression in response to infection. Genome Research. 2014.
      • 68. Beck AH, Weng Z, Witten DM, Zhu S, Foley JW, Lacroute P, et al. 3′-End Sequencing for Expression Quantification (3SEQ) from Archival Tumor Samples. Ng IOL, editor. PLoS One. Public Library of Science; 2010;5: e8768. pmid:20098735
      • 69. Friedman RC, Farh KK-H, Burge CB, Bartel DP. Most mammalian mRNAs are conserved targets of microRNAs. Genome Research. Cold Spring Harbor Lab; 2009;19: 92–105. https://doi.org/10.1101/gr.082701.108 pmid:18955434