Primary macrophages and J774 cells respond differently to infection with Mycobacterium tuberculosis

Macrophages play an essential role in the early immune response to Mycobacterium tuberculosis and are the cell type preferentially infected in vivo. Primary macrophages and macrophage-like cell lines are commonly used as infection models, although the physiological relevance of cell lines, particularly for host-pathogen interaction studies, is debatable. Here we use high-throughput RNA-sequencing to analyse transcriptome dynamics of two macrophage models in response to M. tuberculosis infection. Specifically, we study the early response of bone marrow-derived mouse macrophages and cell line J774 to infection with live and γ-irradiated (killed) M. tuberculosis. We show that infection with live bacilli specifically alters the expression of host genes such as Rsad2, Ifit1/2/3 and Rig-I, whose potential roles in resistance to M. tuberculosis infection have not yet been investigated. In addition, the response of primary macrophages is faster and more intense than that of J774 cells in terms of number of differentially expressed genes and magnitude of induction/repression. Our results point to potentially novel processes leading to immune containment early during M. tuberculosis infection, and support the idea that important differences exist between primary macrophages and cell lines, which should be taken into account when choosing a macrophage model to study host-pathogen interactions.

Scientific RepoRts | 7:42225 | DOI: 10.1038/srep42225 to be genetically unstable, with hybrid phenotypes and possibly atypical signalling mechanisms. Therefore, the question remains on whether cell lines are physiologically relevant systems, particularly for host-pathogen interaction studies.
Previous comparisons of macrophage infection models for M. tuberculosis have focused on certain phenotypic traits such as intracellular growth, cytokine production, or phagosome characteristics 5,9,10 . Jordao and collaborators found that Mycobacterium bovis grows more efficiently in hMDMs than in J774 over a 7-day time course, even though both types of macrophages can restrict growth during the first 24 hours of infection through a mechanism involving phagosome acidification and nitric oxide production 9 . Mendoza and collaborators compared M. tuberculosis, M. bovis and M. bovis BCG infection of hMDMs, THP-1 and U937, and found that U937 has a lower phagocytic capacity than the other two. They also compared production of nitrite, TNF-α and IL-12 in resting and IFN-γ activated macrophages, and found that THP-1 resembles more closely hMDMs 10 .
Here we focus on two mouse macrophage models that are widely used to study M. tuberculosis infection: BMDMs (primary macrophages) and J774A.1 (a cell line, hereafter referred to as "J774" for simplicity). We conduct an in depth analysis of their transcriptomic response to infection with live M. tuberculosis and stimulation with γ -irradiated (killed) M. tuberculosis. We find that BMDMs respond very strongly to infection whereas the response of J774 is delayed and more discreet both in terms of number of differentially expressed genes and magnitude of induction/repression. Furthermore, we identify macrophage genes whose possible function in resistance to M. tuberculosis infection has not yet been investigated.

Results
Global transcriptome profiles. BMDMs and J774 macrophages were either infected with live M. tuberculosis H37Rv or stimulated with γ -irradiated (killed) M. tuberculosis H37Rv for 4 hours. Total RNA was isolated at 4 and 24 hours post-infection (hpi), and CFUs were determined by plating dilutions of lysed macrophages. At 24 hpi, a 0.34 log and a 0.2 log decrease in CFUs was observed in BMDMs and J774, respectively (see Supplementary Fig. S1). Total RNA was used to construct individually barcoded, strand-specific RNA-seq libraries after depletion of ribosomal and mitochondrial RNA, and the libraries were sequenced on an Illumina ® HiSeq ™ 2000 apparatus.
To investigate general trends in the data, principal component analysis (PCA) was carried out (see Supplementary Fig. S2). PCA confirmed no outlier samples and showed separate grouping of uninfected macrophages (controls), and macrophages infected with live bacteria or stimulated with dead bacteria. The transcriptome profiles of BMDMs infected with live bacteria or stimulated with dead bacteria were similar to each other at 4 hpi and became more distinct (from each other and especially from the 4-hpi samples) at 24 hpi. For the J774 cells, PCA showed grouping of the transcriptomes of cells infected with live bacteria or stimulated with dead bacteria at 4 hpi and those of cells stimulated with dead bacteria at 24 hpi; that is, in the case of J774 cells, the 24-hpi transcriptome of cells stimulated with dead bacteria was not clearly different from the 4-hpi transcriptomes.
In addition, the transcriptome of BMDMs uninfected controls (but not that of J774 uninfected cells) changed considerably from 4 hpi to 24 hpi ( Supplementary Fig. S2, and Supplementary Data S1). That is, considerable changes occurred in the transcriptome of uninfected BMDMs while they were maintained in culture, emphasising the need to process a matched non-infected control for each time point when using primary macrophages. These differences in the transcriptome of in vitro grown primary macrophages over time had been previously observed by others both with murine macrophages 11 and human macrophages 12 .

Differences between the transcriptional responses of BMDMs and J774 cells.
Differential expression analysis demonstrated a clear response of M. tuberculosis-infected BMDMs as early as 4 hpi. A total of 2,494 significantly differentially expressed genes (DEGs) with an absolute fold change (FC) greater than 2 (|FC|> 2 or |Log 2 FC|> 1) were detected at 4 hpi in BMDMs infected with live bacteria compared to uninfected BMDMs (FDR < 0.05) (Fig. 1a, and Supplementary Data S2). By contrast, only 234 DEGs were detected for J774 when Of these 234 DEGs, 183 were also differentially expressed in BMDMs ( Fig. 2a and b, and Supplementary Data S3), representing 78% and 7.4% of the DEGs in J774 and BMDMs, respectively. Furthermore, 182 of the 183 common DEGs were differentially expressed in the same direction, with 122 upregulated (top right quadrant of the scatter plot in Fig. 2a), and 60 downregulated (bottom left quadrant, Fig. 2a). However, the absolute fold change in expression of 136 (74.3%) of the 183 common DEGs was higher in BMDMs than in J774 (Fig. 2c). This was a general trend, with a Log 2 FC range between − 5.39 and + 9.76 in BMDMs, whereas in J774 it was narrower, between − 2.33 and + 5.54.
The number of DEGs in J774 (but not in BMDMs) increased considerably from 4 to 24 hpi, but was still lower than the number of DEGs in BMDMs (Fig. 1, and Supplementary Data S2). 918 genes with a |FC|> 2 were differentially expressed at 24 hpi in both BMDMs and J774 cells ( Fig. 2d and e, and Supplementary Data S3). These represent 41% and 54% of the total DEGs at 24 hpi, respectively. Of these 918 common DEGs, 903 were differentially expressed in the same direction, with 398 upregulated (top right quadrant of the scatter plot), and 505  Fig. 2d). At 24 hpi, the level of response was more similar between the two types of cells, with only 55% of the common DEGs presenting a Log 2 FC significantly different (Fig. 2f). Overall, the global level of induction/repression was similar at 24 hpi in both types of macrophages ranging from − 6.00 to + 9.77 in BMDMs, and − 5.41 to + 9.84 in J774. Of note, 104 of the common DEGs detected at 4 hpi were also differentially expressed in both cell types at 24 hpi, including most of the top differentially expressed ones.
Core response of murine macrophages to infection with live M. tuberculosis. A more detailed examination of the genes differentially expressed in both J774 and BMDMs highlighted many genes previously described as important during M. tuberculosis infection (Supplementary Data S3). These include receptors for the recognition and phagocytosis of M. tuberculosis, such as the C-type lectin receptor Mincle (Clec4e) and the class A scavenger receptor MARCO, both of which were induced; the mannose receptor (Mcr1), that was repressed; and the cytosolic receptors Nod1 and Nod2, and Tlr2, induced. The products of other shared induced genes are involved in the activation and recruitment of antimicrobial mechanisms, like the CD40 receptor that activates antimicrobial mechanisms in infected macrophages upon interaction with ligand CD40L in T-cells 13 ; several IFN-inducible GTPases such as Gbp6, Gbp7, Irgm1 (Lrg-47) and Irgm2, which confer cell-autonomous immunity to mycobacterial infection 14,15 ; Isg15, which induces IFNγ -production in lymphocytes 16 ; Irg1, with bactericidal properties against M. tuberculosis through various mechanisms including inhibition of its isocitrate lyase [17][18][19] ; Il1b, which stimulates antimycobacterial immunity in macrophages 20,21 , and caspase-1 (Casp1), required for processing of pro-IL-1β into its active form 22 . Mmp9 and Ccl5 (Rantes), associated with recruitment of monocytes and macrophages during granuloma formation 23,24 , and the neutrophil and lymphocyte-recruiting chemokines Cxcl2, Cxcl10 and Cxcl11 were also induced; as were the pro-apoptosis genes Ptgs2 (Cox2) and Ptges (Pges), while the pro-necrosis gene Alox5 was repressed.
Several interferon-stimulated genes (ISGs) were also highly induced, such as Rsad2 (Viperin), Mx1, Mx2, Ifit1, Ifit2, Ifit3, Ifi205, Apol9a, and Apol9b, all with known antimicrobial functions 25-27 but without a known function during M. tuberculosis infection. Interestingly, the gene encoding RIG-I (Ddx58) was also induced. RIG-I is a pattern recognition receptor involved in viral RNA recognition, critical for the host antiviral response. The top downregulated genes included Rtn4rl1; the apoptosis activator Bmg; the receptor Epha2, whose absence is associated with a reduced bacterial load during the chronic phase of M. tuberculosis infection in mice 28 ; Gpr34, whose deficiency in mice has been associated with an altered immune response 29,30 ; the chemokine Ccl24; and the scavenger receptor Fcrls.

Pathway analysis of DEGs in macrophages infected with live M. tuberculosis. To gain insight into the cellular and molecular functions of the M. tuberculosis-induced host genes, DEGs in macrophages infected
with live bacteria were subjected to pathway analysis using QIAGEN's Ingenuity ® Pathway Analysis (IPA ® ).
A full list of all the enriched pathways for each condition can be found online in the Supplementary Data S4. Comparison analysis was used to visualise the results across all four conditions. The top enriched pathways are shown in Fig. 3. Pathways are ranked according to the enrichment score (Fisher's exact test P-value, Fig. 3a), and to the Z-score which predicts activation or inhibition of the pathways by comparing observed and predicted regulation patterns (Fig. 3b). The same pathways were enriched in BMDMs and J774 although with different scores. For most pathways, the enrichment score and the Z-score was higher in BMDMs than in J774, particularly at 4 hpi, reflecting the stronger response in this type of macrophages early after the infection (Fig. 3).
TREM1 (triggering receptor expressed on myeloid cells 1) was one of the most significantly enriched pathways, with a high activation score in all conditions. TREM1 is a DAP12-associated receptor that plays an essential role in innate immunity by fine-tuning the inflammatory response 31 . In our dataset, many of the genes encoding cytokines and chemokines whose expression and secretion is induced by TREM1 activation, as well as several genes encoding cell surface proteins, were upregulated (see Supplementary Fig. S3). In addition to regulating the innate immune response, TREM1 modulates the adaptive immune response by inhibiting the regulators SIGIRR, which downregulates the Th1 response, and ST2, which enhances the Th2 response 32 . Both SIGIRR and ST2 were downregulated in our dataset.
Several pathways involved in bacterial recognition by host cells were also activated such as "Activation of IRF by cytosolic pattern recognition receptors", "Role of RIG-1 like receptors in antiviral innate immunity", "Toll-like receptor signalling", or the more general pathway "Role of pattern recognition receptors in recognition of bacteria and virus" which includes TLR, NLR, RIG and complement receptors (Fig. 3).
Only a few pathways had a negative Z-score (inhibition), including PPAR signalling, LXR/RXR activation, and antioxidant action of vitamin C (Fig. 3b). Of these, PPAR signalling had the highest enrichment score. The PPAR (proliferator-activated receptor) family are ligand-activated transcriptional regulators activated by fatty acids and their derivatives. They regulate the expression of genes involved in lipid metabolism, and their function is supressed by cytokines such as IL-1 and TNF-α .
Cell type-specific response to M. tuberculosis infection. A full list of the genes differentially expressed exclusively in one of the two macrophage models at a particular time point can be found online (Supplementary Data S3). Some of the genes induced exclusively in BMDMs included the chemokines Cxcl1, Cxcl3 and Cxcl9, the cytokines Il1a, Il6, and Il12b; Rasgrp1, which participates in T cell activation; Mmp14, which regulates monocyte migration and collagen destruction in tuberculosis 33 ; Serpinb2, also known as plasminogen activator inhibitor type 2 (Pai-2), a major product of activated monocytes/macrophages which has been linked to regulation of Th1 responses 34 ; and Fpr2, a recently described pathogen recognition receptor that recognizes bacterial signal peptides and triggers classical innate immune responses, such as intracellular Ca 2+ mobilization, generation of reactive oxygen species, release of metallopeptidase, and chemotactic cell migration 35  a Notch target gene that suppress TLR-triggered pro-inflammatory cytokine production 36 , was one of the most downregulated genes in BMDMs.
Among the top induced genes solely in J774 there was Nos2, highly induced at 24 hpi; Rab15, involved in endocytic recycling 37,38 ; the gene encoding for the laminin subunit gamma-2 (Lamc2); Kdm6b also known as Jmjd3, a lysine-specific demethylase that modulates pro-inflammatory responses in macrophages 39 and has been recently implicated in foamy macrophage formation during mycobacterial infection infection 40 ; the IFN-inducible GTPase Gbp1; and Il23a, required for long term control of M. tuberculosis 41 (Supplementary Data S3). Tlr5 and Tlr8 were some of the downregulated genes exclusively in J774, together with Il10 that was repressed at 24 hpi in J774 while it was induced in BMDMs at 4 hpi.
Pathway analysis did not reveal any cell-specific pathways, although at 4 hpi some pathways were enriched exclusively in BMDMs, which is not unexpected given the modest response exhibited by J774 cells at this early time point (Fig. 3). Some of the pathways induced in BMDMs at 4 hpi but not in J774 until 24 hpi were "Interferon signalling" and "Role of RIG1-like receptors in antiviral innate immunity".
Immune response of macrophages to stimulation with γ-irradiated M. tuberculosish. To study if the differences detected between both cell types were due to differences in the response to either phagocytosis or infection, the transcriptome of macrophages stimulated with dead (γ -irradiated) M. tuberculosis was also analysed by RNA-seq. Similar to what was observed when infecting with live bacteria, differential expression analysis of macrophages stimulated with dead bacteria compared to uninfected macrophages at 4 hpi showed a stronger response in BMDMs (1,764 DEGs) than in J774 (211 DEGs) (Fig. 1). At 24 hpi the number of DEGs slightly decreased in BMDMs (1,338 DEGs) and increased in J774 (337 DEGs). Moreover, the magnitude of regulation in macrophages infected with live bacteria and stimulated with dead bacteria was similar in BMDMs at both time points, and in J774 at 4 hpi, but lower in J774 stimulated with dead bacteria at 24 hpi. Specifically, the Log 2 FC in Pairwise comparisons of DEGs in macrophages stimulated with dead bacteria or infected with live bacteria showed a marked overlap, particularly in BMDM. 79% and 82% of the DEGs in BMDMs stimulated with dead bacteria were also differentially expressed in BMDMs infected with live bacteria at 4 hpi and 24 hpi, respectively ( Fig. 4a and b). These correspond to 58% and 48% of the DEGs in BMDMs infected with live bacteria, respectively. Similarly, 64% and 65% of the DEGs in J774 stimulated with dead bacteria were also differentially expressed in J774 infected with live bacteria at 4 hpi and 24 hpi, respectively ( Fig. 4c and d). These correspond to 57% of the DEGs in J774 infected with live bacteria at 4 hpi, but only 13% of the DEGs in J774 infected with live bacteria at 24 hpi ( Fig. 4c and d). This reflects the substantial difference in the total number of DEGs between the two conditions in J774 at 24 hpi (Fig. 1b). Nevertheless, in all four comparisons virtually all the common genes were differentially expressed in the same direction and magnitude (Fig. 4). Consequently, pathway analysis results were similar in macrophages stimulated with dead bacteria and macrophages infected with live bacteria, and all the pathways had the same direction of activation/inhibition prediction (Fig. 5). As expected, the degree of activation/inhibition (Z-score absolute value) was similar in macrophages stimulated with dead bacteria or infected with live bacteria at 4 hpi but lower in macrophages stimulated with dead bacteria at 24 hpi. Overall, these results indicate that the differences detected in the transcriptomic responses between the two cell types occur even just upon phagocytosis of dead M. tuberculosis.

Discussion
Macrophages play an essential role in the early immune response against M. tuberculosis and are the cell type preferentially infected in vivo. The failure of an appropriate macrophage response is critical for the establishment of infection and progression of the disease. M. tuberculosis has evolved multiple evasion mechanisms in order to survive and grow inside macrophages. To study the antimicrobial response of macrophages and the evasion mechanisms of M. tuberculosis both primary macrophages and cell lines are commonly used. BMDMs are a popular model because they are easy to obtain, and can be sourced from a variety of specific genetic backgrounds. Likewise, the cell line J774 is frequently used since it is readily available, and no animals are involved, which results in a reduction in the number of animals used in research. In the present study, we compare the transcriptomic response of BMDMs and the macrophage-like cell line J774 to infection with M. tuberculosis H37Rv by using high-throughput RNA-seq.
We have found that BMDMs respond very strongly to infection whereas the response of J774 to M. tuberculosis infection is weaker both in terms of number of DEGs and magnitude of induction/repression, at least during the early time points of the infection. A strong host response similar to that obtained in BMDMs has been previously observed using microarrays in murine BMDMs infected with live M. tuberculosis strains CDC1551 and HN878 by Koo et al. 42 . A high number of DEGs were detected early in the infection (6 hpi) followed by a decrease in the number of DEGs at 24 hpi 42 . Other studies using human or bovine MDMs have reported a more modest response in the first hours of the infection followed by a marked increase in the amount of DEGs detected at around 18-24 hpi, similar to what we have observed in the J774 cell line [43][44][45] . Interestingly, previous work has shown that infection with Yersinia enterocolitica results in a higher number of upregulated genes in BMDMs than in J774 cells 46,47 . This seems to be due to a higher susceptibility of J774 to inhibition by specific anti-host effector proteins from Y. enterocolitica that interfere with host signalling pathways 47 . Whether a similar scenario occurs upon infection with M. tuberculosis would be difficult to establish considering the multifactorial nature of its virulence. An additional explanation would be that J774, which derive from a lymphoma, are already in a relatively activated state with enhanced glycolysis, which may have contributed to the reduced response observed at 4 hpi compared to the more naïve BMDM.
It is important to note that the BMDMs used in the present study were obtained from C57BL/6 mice, an extensively studied mouse strain from which most knockout mutants derive, whereas the J774 cell line derives from Balb/c mice. Some differences have been described in the immune response of these two mouse strains to M. tuberculosis infection and BCG vaccination in vivo, entailing mainly a greater early innate response (measured by production of IL-12, IFN-γ , TNF-α , MCP-1 and neutrophil influx) but a suppressed Th1 response in Balb/c compared to C57BL/6 48,49 . However, both strains control the infection similarly and have been classified as resistant to M. tuberculosis infection in comparison to other more susceptible mouse strains 50 . Besides, they both mount a comparable protective immunity upon vaccination with BCG 51 . Additionally, their BMDMs respond very similarly to M. tuberculosis infection in vitro 52,53 . However, we cannot rule out that the differences we have observed are the result of both the different origin of the cells and of them being primary cells and a cell line.
Except for the more modest response of J774 cells, the type of response of both macrophage types to M. tuberculosis infection was quite similar in terms of DEGs and pathways enriched, particularly at 24 hpi. Importantly, several previously described antimycobacterial mechanisms were induced in our dataset such as Irgm1 (Lrg-47) and Irgm2, or Isg15. Many other ISGs with known antimicrobial function were also induced, such as Rsad2 (Viperin) or Ifit1/2/3. To the best of our knowledge, a specific role for many of these genes in M. tuberculosis infection has not been described. Therefore, the study of their involvement in the defence against M. tuberculosis is worth exploring.
An interesting pathway that was highly induced in both macrophages was TREM1. TREM1 has been shown to regulate the host response at an early stage during infection of various pathogens such as Streptococcus pneumoniae 54 , Klebsiella pneumoniae 55 , and Streptococcus suis 56 , and it is induced in humans during active tuberculosis 57 . Although the ligand of TREM1 is unknown, it has been speculated that pathogen-associated molecular patterns could directly activate it 31 , therefore it would be interesting to study if any of the M. tuberculosis cell wall components is a ligand for this receptor.
The PPAR signalling pathway was repressed. Current evidence indicates that mycobacterial infection causes a time-dependent increase in PPARγ expression through a macrophage mannose receptor-dependent pathway 58 . This induction results in increased lipid droplet formation and down-modulation of the macrophage response, suggesting that PPARγ expression might aid the mycobacteria in circumventing the host response. In fact, knockdown of PPARγ in human macrophages has been shown to result in an increase in TNF production and a decrease in M. tuberculosis growth at 24 hpi 58 . The inhibition of this pathway together with the antimicrobial mechanisms previously mentioned could explain the control of M. tuberculosis growth by the macrophage at this early stage of the infection 9,59-61 . PPAR is a well-known pathway associated with diabetes 62 . Interestingly there is now a recognised connection between diabetes and tuberculosis 63 , although the mechanisms leading to this interaction are still unclear 64,65 . However, with the PPARγ pathway common between these diseases, it could be worthwhile to explore this connection in more detail.
It has been recently recognised that M. tuberculosis derived-products reach the cytosol early during infection via the formation of pores in the phagosomal membrane in an ESX1-dependent process 66 . Both M. tuberculosis DNA and cyclic-di-adenosine monophosphate (c-di-AMP) are recognized by cytosolic receptors triggering production of type I IFN, ISGs and pro-inflammatory cytokines such as IL-1β through the (cGAS)-STING-TBK1-IRF3 and AIM2-NLRP3-IL-1β signalling pathways 67,68 . Surprisingly, the STING/TBK1 axis was not differentially expressed in our dataset but the genes for other DNA sensors such as Dai (Zbp) and Ifi204, and for RNA sensors such as RIG-I (Ddx58), LGP2 (Dhx58), and MDA5 (Ifih1) were induced, as were the downstream transcription factors NF-κ B, IRF7, IRF9, STAT1 and STAT2, the effectors IFNβ and IL-1β , and many of the ISGs, particularly in BMDMs (at both time points, and either upon infection with live bacteria or stimulation with dead bacteria), but only marginally in J774 infected with live bacteria at 24 hpi. The importance of cytosolic nucleic acid receptors in the antimicrobial defence against bacteria has been increasingly recognized 69,70 . The relevance of the activation of these signalling pathways during M. tuberculosis infection and the origin of the DNA/RNA that they would sense deserves further investigation.
The similarity in the host response in macrophages infected with live bacteria or stimulated with γ -irradiated bacteria at this early stage of the infection is not surprising if we consider that, contrary to what happens with heat-killing, inactivation by γ -irradiation preserves cell integrity 71 . Therefore, the same pathogen associated molecular patterns (PAMPs) are recognized by the cell's pattern recognition receptors (PRRs) leading to the same downstream signalling pathways being activated 72 . Furthermore, several studies have shown that cells inactivated by γ -irradiation are still metabolically active and de novo transcription and translation still occurs 71,73 . Consequently, it would be safe to speculate that during macrophage stimulation with γ -irradiated M. tuberculosis the early host-pathogen interactions are conserved. This would most likely include activation of the cytosolic surveillance pathway through ESX-1-mediated cytosolic access, in the same way that γ -irradiated Francisella tularensis is able to escape the endosome into the cytoplasm 74 . In fact, many of the hallmark genes of the cytosolic surveillance pathway were highly induced in BMDMs stimulated with dead bacteria. The induction was particularly high at 4 hpi and it decreased at 24 hpi (Log 2 FC ≈ 1-2) while it was still high in macrophages infected with live bacteria (Log 2 FC ≈ 3-4). It is likely that, at later time points, macrophages stimulated with inactivated bacteria incapable of replication can digest the bacteria and slowly shutdown their response. This would also explain the low number of DEGs at 24 hpi in macrophages stimulated with dead bacteria compared to macrophages infected with live bacteria.
In conclusion, comparison of the response of bone marrow-derived primary macrophages and the macrophage-like cell line J774 has highlighted differences in the timing and strength of the response. These differences should be taken into account when using a macrophage model for the study of host-pathogen interactions. Furthermore, the use of high-throughput RNA-seq has allowed delineation of the macrophage transcriptomic Scientific RepoRts | 7:42225 | DOI: 10.1038/srep42225 response to M. tuberculosis infection in great depth, identifying genes whose function in resistance to M. tuberculosis infection has not yet been investigated. These findings point to potentially novel processes leading to immune containment early in the infection.

Infection of BMDMs and J774 cells with
Strand-specific RNA-seq library preparation and sequencing. Approximately 500 ng of RNA from each sample was used to prepare individually barcoded strand-specific RNA-seq libraries. rRNA and mitochondrial RNA was removed using Ribo-Zero Gold Removal Magnetic Kit for human/mouse/rat (Epicentre, now Illumina) following the manufacturer's instructions. mRNA was purified using Agencourt RNAClean XP kit.
Libraries were constructed using ScriptSeq ™ Complete Kit (Epicentre, now Illumina) as per manufacturer's instructions. Quality and quantity of the libraries were determined using a Bioanalyzer (Agilent) and Qubit ® (Invitrogen), respectively. Libraries were then quantified by qPCR in an Applied Biosystems ® 7500 with the Kapa Biosystems kit for library quantitation (Kapa Biosystems). Next-generation sequencing was performed using an Illumina HiSeq 2000 flow cell with one 150-bp end run. The NCBI Gene Expression Omnibus (GEO) accession number for the RNA-seq data reported in this paper is GSE88801.

RNA-seq data analysis.
Reads were aligned versus the Ensembl mouse genome version mm10 (GRCm38) using STAR 75 Figure S5). Differential expression was analysed in SeqMonk after quantification of raw counts, strand-specific, merged isoforms, using the R (version 3.2.2) package DESeq2. A gene was classified as up-or downregulated using a cut-off value of more or less than 2-fold expression difference, with an adjusted P-value < 0.05, following multiple testing correction. DEGs were further analysed using QIAGEN's Ingenuity ® Pathway Analysis (IPA ® , QIAGEN) to identify over-represented canonical pathways.