Dynamic changes in host gene expression associated with H5N8 avian influenza virus infection in mice

Emerging outbreaks of newly found, highly pathogenic avian influenza (HPAI) A(H5N8) viruses have been reported globally. Previous studies have indicated that H5N8 pathogenicity in mice is relatively moderate compared with H5N1 pathogenicity. However, detailed mechanisms underlying avian influenza pathogenicity are still undetermined. We used a high-throughput RNA-seq method to analyse host and pathogen transcriptomes in the lungs of mice infected with A/MD/Korea/W452/2014 (H5N8) and A/EM/Korea/W149/2006 (H5N1) viruses. Sequenced numbers of viral transcripts and expression levels of host immune-related genes at 1 day post infection (dpi) were higher in H5N8-infected than H5N1-infected mice. Dual sequencing of viral transcripts revealed that in contrast to the observations at 1 dpi, higher number of H5N1 genes than H5N8 genes was sequenced at 3 and 7 dpi, which is consistent with higher viral titres and virulence observed in infected lungs in vivo. Ingenuity pathway analysis revealed a more significant upregulation of death receptor signalling, driven by H5N1 than with H5N8 infection at 3 and 7 dpi. Early induction of immune response-related genes may elicit protection in H5N8-infected mice, which correlates with moderate pathogenicity in vivo. Collectively, our data provide new insight into the underlying mechanisms of the differential pathogenicity of avian influenza viruses.

Scientific RepoRts | 5:16512 | DOI: 10.1038/srep16512 to become transmissible between species and to directly infect humans, more detailed studies on the characterisation of the pathogenesis of H5N8 virus infection are necessary.
Our previous work demonstrated that the MDk/W452(H5N8) strain is less pathogenic than is the EM/W149(H5N1) strain in BALB/c mice 7 . In addition, MDk/W452(H5N8) was able to replicate in human nasal epithelium and lung tissues, although its replication efficiency was less than that of EM/ W149(H5N1). Analysis of mouse bronchial alveolar fluids indicated that MDk/W452(H5N8) infection moderately upregulated the expression of cytokines/chemokines, whereas the HPAI A(H5N1) viruses, particularly EM/W149(H5N1), induced robust expression of all the pro-inflammatory cytokines/chemokines, including tumour necrosis factor-α (TNF-α ), interleukin-1β (IL-1β ), IL-6, regulated on activation, normal T cell expressed and secreted (RANTES), and granulocyte macrophage colony-stimulating factor (GM-CSF). Consistent with this finding, persistent activation and recruitment of immune cells (e.g. T cells, macrophages, and monocytes), the hallmark of severely lethal HPAI A(H5N1) virus infection, was not observed with A(H5N8). Thus, in comparison with H5N1-infected mice, H5N8-infected mice showed moderate pathogenicity, although detailed mechanisms of pathogenicity are yet to be determined.
A better understanding of the global gene changes underlying the multi-step progression of pathogenicity during infection could help develop potential therapeutic strategies for emerging infectious diseases 8 . Very recently, high-throughput RNA sequencing (RNA-seq) technology, which is a powerful way to profile the transcriptome with great efficiency and higher accuracy, has been employed in various viral infections and diseases [9][10][11][12][13][14] . RNA-seq technology has the potential to reveal the alteration dynamics of the pathogen genome itself and the systemic change in host gene expression in the process of infection by pathogens, which could help uncover the pathogenesis of the infection and interaction mechanism of pathogens.
In the present study, we used the RNA-seq technology to comprehensively annotate dual host/pathogen transcriptomes following infection of mice with H5N8 and H5N1. Our work provides a global view of virus type-specific and post-infection time-specific mRNA profiles, which illustrate the role and mechanisms of viral replication efficiency and host response during avian influenza infection of mice.

Results
Analysis of the high-throughput sequencing RNA-seq transcriptome. To gain global and dynamic gene expression profiles of A/MD/Korea/W452/2014-(452, H5N8)-and A/EM/Korea/ W149/2006-(149, H5N1)-infected BALB/c mice, RNA-seq was performed to explore the transcriptomes from mouse lungs derived at 1, 3, and 7 day post infection (dpi). H5N1 was used to compare the pathogenicity with recently isolated H5N8 virus. More than 37 million 100-bp paired-end reads were generated. The base quality of RNA-seq reads was checked and analysed using an Agilent RNA 6000 Nano Kit ( Supplementary Fig. 1a). Each sample had relatively high sequencing coverage range, as shown by pie charts (Supplementary Fig. 1b). The uniquely mapped reads covered over 80% of the total genomic bases ( Supplementary Fig. 1c). The read depth seemed to be distributed relatively evenly along the whole body of the genes, reflecting no introduction of obvious bias during randomly primed reverse transcription and subsequent RNA sequencing.
Using the 41,388 annotated genes in the mouse genome database, gene expression was quantified and compared between the control and the virus-infected groups, and differentially expressed genes (DEGs) with a q-value < 0.05 were identified. Figure 1A shows total numbers of both up-and downregulated DEGs for H5N8-and H5N1-infected mice (with + /− 2-fold changes). Interestingly, the number of upregulated genes was larger than the numbers of downregulated genes in both groups of virus-infected mice at both time points, except H5N1 at 7 dpi (Fig. 1A). Furthermore, Venn diagrams were generated to examine the overlapping mRNA profiles for H5N8-and H5N1-infected mice. mRNAs that were up and down regulated at least 2 fold at each time point were used to generate the Venn diagrams. The mRNA differential expression levels in H5N8-and H5N1-infected mice are depicted as three overlapping circles, as shown in Fig. 1A. The numbers indicate the mRNA counts in the indicated area. At 1 dpi, only 3 genes were commonly overlapping upregulated in H5N8-and H5N1-infected mice (Fig. 1B). While higher numbers of genes were upregulated in the lungs of H5N8-infected mice than in those of H5N1-infected mice at 1 dpi, the opposite pattern was observed at 3 and 7 dpi. In H5N1-infected mice, 384 and 760 genes were upregulated more than 2-fold at 3 and 7 dpi, respectively, whereas only 57 and 129 genes were upregulated at 3 and 7 dpi, respectively, in H5N8-infected mice. Overall, there were 245 genes from H5N1 and 492 genes from H5N8 overlapping at 3 and 7 dpi, respectively.
To illustrate that the general expression pattern of transcripts was similarly distributed across all individual biological samples within the group, a volcano plot was generated (Fig. 1C). The ratio of the differential expression (fold change, illustrated in the abscissa of the volcano plot) shows that there was very good correlation between the fold change differences and p-values (i.e. genes with a large fold change difference also had a low p-value in the group-wise comparison). Additionally, further analysis of DEGs were visualized by hierarchical clustering analysis, scatter plot and MA plot ( Supplementary  Fig. 2a-c). These data indicate a good clustering of samples according to levels of similarities in the gene expression pattern and there is a clear distinction between these H5N8-and H5N1-infected mice, and the gene expression profiles were able to discriminate between the two main groups of virus-infected mice. To evaluate the degree of dissimilarity among samples in regard to biological variations and dimensions, we then compared H5N8-and H5N1-infected mouse transcriptomes obtained from RNA-seq Scientific RepoRts | 5:16512 | DOI: 10.1038/srep16512 Figure 1. Global overview of the RNA-seq data of H5N8-and H5N1-infected mice lungs. A global overview of RNA-seq data of mice infected with mock, H5N8 virus, or H5N1 virus is shown. Lungs from three mice per group were harvested at 1, 3, and 7 days post infection (dpi). RNA-seq was performed for 21 samples. (A) Numbers of more than 2 fold up or down-regulated differentially expressed genes (DEGs) identified from the comparison among mock and virus-infected groups (DEGs were identified based on a false discovery rate (FDR) q-value threshold of less than 0.05). (B) Venn diagrams of overlapping DEG profiles for H5N8 and H5N1. DEGs are those displaying a change of more than two-fold with a p-value of less than or equal to 0.05. The mRNA differential expressions in H5N8-and H5N1-infected mice are depicted in three overlapping circles for 2-fold up and down-regulation at 1, 3, and 7 dpi. The numbers indicate the mRNA counts in the indicated area. (C) Volcano plot showing DEGs for H5N1-and H5N8infected mice. The x-axis represents the log 2 values of the fold change observed for each mRNA transcript, and the y-axis represents the log 10 values of the p-values of the significance tests between replicates for each transcript. Data for genes that were not classified as differentially expressed are plotted in black. by performing three-dimensional multidimensional scaling (3D-MDS) ( Supplementary Fig. 3). Overall, H5N1 and H5N8 samples at 1 dpi were distinctly distributed from H5N1 and H5N8 samples at 3 and 7 dpi in the 3D-MDS plot, suggesting that these samples are closely clustered according to virus types and post infection time points. Viral transcript levels correlate with changes in host gene expression. RNA-seq enabled dual analysis of both host and viral transcripts within the same sample. Paired-end reads from influenza-infected samples were mapped to the influenza virus genome. Interestingly, viral transcript counts of hemagglutinin (HA), neuraminidase (NA), and nonstructural protein 1 (NS1) of H5N8 were found to be significantly higher than those of H5N1 at 1 dpi. Given that host transcriptome data also showed more upregulated DEGs for H5N8 in comparison with H5N1 at 1 dpi, these data imply a strong correlation between sequenced numbers of viral transcripts and expression levels of host immune-related genes. Overall, the pattern of influenza virus gene expression was different between H5N8-and H5N1-infected mice in a time-dependent manner. The average number of reads was approximately 10 to 100 times higher in H5N1-infected mice than in H5N8-infected mice at 3 and 7 dpi (Fig. 2). Significantly higher viral expression was observed in the H5N1 group compared with the H5N8 group for every viral gene at 3 dpi (p < 0.05). At 7 dpi, the transcript counts for HA, M, NP, PA, and PB1 were found to be significantly different between the H5N1-and H5N8-infected mice. These results suggest that viral gene expression levels also vary for these avian influenza viruses in a temporal dependent manner.
Distinct and dynamic changes in host DEGs in H5N8-vs H5N1-infected mice. To further analyse the characteristics of DEGs, we focused on the 10 DEGs showing the most marked up or down-regulation from H5N8-and H5N1-infected mice ( Table 1 and 2). At 1 dpi, there was no significantly upregulated immune-related gene in H5N1-infected mice. However, elevated upregulation of a key immune modulator, chemokine (C-X-C motif) ligand 3 (CXCL-3) was detected in H5N8-infected mice. Interestingly, in contrast to that observed at 1 dpi, the number of DEGs of immune-related function increased dramatically at 3 and 7 dpi. Chemokine pathway-associated genes, such as C-X-C motif chemokine 10 (CXCL-10), CXCL-9, chemokine (C-C motif) ligand 7 (CCL-7), and CXCL-2, were among the 10 most upregulated genes for both H5N8-and H5N1-infected mice at 3 dpi. Furthermore, at 7 dpi, CXCL-11 was found to be one of the most upregulated genes in the H5N8-and H5N1-infected mice, although its expression was higher in H5N1-infected mice than in H5N8-infected mice (10.4 vs 8.76 log 2 -fold change (fc) values). We also determined DEGs that were upregulated in virus-infected groups only but not in mock-infected groups. Among those genes, IL-24, which is known to be an important regulator of tumorigenesis 15 , was found to be highly upregulated in the H5N8-(4.87 log 2 fc values; Supplementary Table 1) and H5N1-infected mice (5.98 log 2 fc values; Supplementary Table 2). Furthermore, IL-21 and CCL1 were also commonly found to be upregulated in the H5N8-and H5N1-infected mice, but not in the mock-infected group.
To examine the biological roles of the DEGs, a gene ontology (GO) enrichment analysis was applied to the up-regulated genes. GO term statistical analysis of genes in this cluster confirmed that defence response to virus were significantly over-represented (Supplementary Table 3 and 4). Among the biological process category of GO, we focused primarily on immune system processes because many genes related to a function of anti-viral signalling are included in this GO category ( Supplementary Fig. 4). Interestingly, in comparison with H5N8-infected mice, higher numbers of genes were found to be included in the GO term of immune system processes in the H5N1-infected mice at all times.

Comprehensive analysis of functional enrichment and networks in H5N8-vs H5N1-infected mice.
To investigate possible biological interactions of DEGs and identify important functional networks, datasets representing genes with altered expression profiles derived from RNA-seq analyses were imported into the Ingenuity pathway analysis (IPA) tool. The highest activated networks (high z-score) were identified using IPA. Figure 3A represents a hierarchical clustering heatmap showing the list of activated canonical pathways after influenza virus infection. Notably, triggering receptor expressed on myeloid cells 1 (TREM1) signalling was most activated at 3 and 7 dpi in both H5N8-and H5N1-infected mice. TREM1 is known to activate neutrophils and monocytes/macrophages by signalling through the 12-kDa adapter protein DNAX activation protein (DAP12), and it activates toll-like receptor (TLR) signalling pathways via secretion of proinflammatory cytokines and chemokines in response to microbial infection 16 . Similarly, canonical pathway analysis identified several key innate immune receptor pathways for virus sensing, such as death receptor signalling, toll-like receptor signalling, roles of pattern recognition receptors in virus sensing, retinoic acid-mediated apoptosis signalling, and activation of interferon regulatory factor (IRF) by cytosolic pattern receptors. The heat map also indicates higher activation of death receptor signalling in H5N1-infected mice than in H5N8-infected mice at 3 dpi (z-scores were 4.73 and 3.30 for H5N1 and H5N8, respectively). Figure 3B also shows the heatmap of the transcriptomic expression values, restricted to the lists of the top 20 up and downregulated DEGs, found as specific for each viral strain. For each set of specific transcripts, hierarchical clustering have been performed and represented using dendrograms. As with the top 10 upregulated gene lists in Table 1, heatmap also identified the chemokine signalling pathway as a key Scientific RepoRts | 5:16512 | DOI: 10.1038/srep16512 pathway involving the majority of DEGs (CXCL-11, CCL7, CXCL-10, CXCL-2, and CXCL-9) across all datasets, based on the number of network interactions. In addition to virus-sensing pathways, key players in innate immunity and inflammation, such as p38 MAPK signalling, and iNOS signalling, NF-κ B signalling, IL-8 signalling, IL-1 signalling, and leukocyte extravasation signalling, were also found to be key pathways involving the majority of DEGs in both H5N8-and H5N1-infected mice at 3 and 7 dpi, although upregulation levels were more enhanced in H5N1-infected mice than in H5N8-infected mice. It is likely that the reason why inflammation pathways were found to be more elevated in H5N1-infected mice than in H5N8-infected mice was because there were higher numbers of statiscially significant DEGs found from H5N1-infected mice than those from H5N8-infected mice at 3 and 7 dpi, which is consistent with higher viral titres and virulence observed in vivo 7 .
Next, the expression of different pathways was analysed with IPA, and several gene interaction networks were created. DEGs due to viral infection were indicated as either red (upregulated) or green (downregulated) (Figs 4 and 5, Supplementary Figs 5 and 6). We focused our analyses primarily on pathways and sets of genes with a previous association with interferon signalling at 3 dpi, as shown in Fig. 4. Notable differences included higher upregulation of key downstream signalling molecules such as Janus kinase/signal transducers and activators of transcription (JAK/STAT), interferon-induced protein with tetratricopeptide repeats (IFIT3), and 2′ -5′ -oligoadenylate synthetase 1 (OAS1) in H5N1-infected mice than in H5N8-infected mice. Similarly, a network map of pattern recognition receptor-recognizing viruses was created to illustrate receptor connections with many chemokines, chemokine receptors, and interferon-inducible genes for H5N8-and H5N1-infected mice at 3 dpi ( Supplementary Fig. 5).
We further generated a network map of death receptor signalling for H5N8-and H5N1-infected mice at 3 dpi (Fig. 5). FasL, TNF-α , and Apo2L (TRAIL)-mediated downstream signalling leading to cell apoptosis were activated to a greater extent in H5N1-infected mice than in H5N8-infected mice. In addition, acute phase signalling was also activated more in H5N1-infected mice than in H5N8-infected mice at 7 dpi ( Supplementary Fig. 6) (z-scores were 5.81 and 3.47 for H5N1 and H5N8, respectively), which also correlates with more severe pathogenicity observed in H5N1-infected mice. H5N1 viruses were very lethal and replicated systemically in various tissues in mice; A significant upregulation of death receptor signalling by H5N1 may have been responsible for the severe pathogenicity observed in the H5N1-infected mice in vivo. Consequently, these data collectively suggest that higher upregulation of DEGs from both viral and host transcriptomes induced by H5N1 infection at 3 and 7 dpi may have  Our data highlight distinct and dynamic changes in expression patterns of host genes following MDk/ W452(H5N8) virus infection in mouse lungs compared with those observed following H5N1 virus infection. Previously, we have demonstrated that intranasal inoculation of mice with 10 6 log 10 EID 50 /mL (EID 50 is 50 percent embryo infectious dose) of MDk/W452(H5N8) induced only 6% reduction in body weight and 40% lethality within the 14-day observation period, whereas EM/W149(H5N1) infection was highly lethal and replicated systemically in various tissues 7 . We observed a more modest secretion of inflammatory molecules in the organs of H5N8-infected mice than in those of H5N1-infected mice. Additionally, quantitative realtime PCR was performed to validate RNA-seq data and the transcript levels of cytokines/chemokines, such as CXCL10, CXCL11, IFN-γ , TNF-α and IL-24, were all up-regulated upon H5N1 and H5N8 infection (data not shown). These data correlate well with the RNA-seq data presented in this study, in which significantly upregulated DEGs were mostly associated with chemokine pathways, and H5N1-infected mice expressed higher levels of chemokines and interferon-stimulated genes (ISGs) than did H5N8-infected mice at 3 and 7 dpi (Fig. 3). A robust activation of immune system activation pathways such as acute phase response, p38 MAPK, and iNOS signalling in H5N1-infected mice correlates with high viral replication in these mice, and could play an important role in pathogenesis. Furthermore, our data demonstrating a more marked activation of the death receptor signalling pathway in H5N1-infected mice also correlates with the high pathogenicity observed with this virus in  (Fig. 5). Notably, when we analysed DEGs found only in the infected groups and not sequenced in the mock-infected group, few interesting novel genes and genes with unknown function in response to influenza virus infection were found. One of the most upregulated genes was IL-24, which is a new member of the IL-10 family of cytokines. Its role in virus-mediated immune responses was unknown. It may be worthwhile to carry out further studies to investigate the role of IL-24 in viral replication control.
Previously, our data have shown that virus titres from H5N8-infected mouse lungs were lower than those from H5N1-infected mouse lungs, suggesting a moderate viral replication pattern of H5N8 in comparison to that of H5N1 7 . Consistent with this, our RNA-seq analysis of virus genes also indicates that viral transcript counts for all eight genes (HA, M, NA, NP, NS, PA, PB1, and PB2) were lower in H5N8 than in H5N1 at 3 and 7 dpi. Interestingly, at 1 dpi, distinct patterns in viral transcript counts were observed in comparison to those observed at 3 and 7 dpi. The transcript counts for HA, NS1, and NA were significantly more elevated in H5N8-infected mouse lungs compared with H5N1-infected mouse lungs at 1 dpi. This raises the important question of whether H5N8's ability to replicate with better efficiency at earlier post infection time points might have led to increased numbers of DEGs being involved in host anti-viral innate immune defence functions to suppress viral replication efficiently, thereby resulting in the consequent modest pathogenicity.
RNA-seq analysis of influenza virus infection has been performed in mouse, crow, and chicken samples 10,21,22 . In mice, pandemic H1N1/2009 virus infection activated strong virus-sensing signals, such as TLRs, NLRs, RIGs, as well as the NF-κ B and JAK-STAT pathways, which play a significant role in inducing innate immunity 21 . Additionally, there are several transcriptome studies, suggesting the importance of understanding dynamics and magnitude of the innate immune responses in shaping the H5N1 pathogenesis [23][24][25][26] . These studies indicate more robust and early induction of inflammatory response observed Pathway analysis with IPA software allowed identification of pathways that were differentially expressed between H5N1and H5N8-infected mice. Genes associated with the interferon signalling canonical pathway that showed differential expression are highlighted in colour. Colour intensity indicates the degree of upregulation (red) or downregulation (green) relative to the mock-infected mice. Solid lines represent direct interactions and dashed lines indirect interactions.
in all animal models infected with H5N1 and up-regulated DEGs from H5N1-infected mice consist of genes involved in inflammation and apoptosis, which may have contributed to underlying mechanisms of severe pathogenicity. Similar to this, our data also indicate that avian influenza infection increased expression of interferon and innate immune virus-sensing pathways.
It was also interesting to note that IFN-γ was robustly expressed in the top 10 upregulated genes of H5N8-infected lungs. Previous studies demonstrate that IFN-γ treatment significantly reduced the number of T and NKT cells in the lungs at the inflammatory phase following infection, thus contributing to protection against viral pathogenesis 27 . In addition, according to Morrison et al., H5N1-infected mice also induced an increased secretion of IFN-γ 24 . Based on our data, it is possible that at early time point (1 dpi), expression of immune-related DEGs from H5N8-infected mice may have facilitated an increase in IFN-γ production, which may have led to abrupt anti-viral immunity and thus limiting viral replication. On the other hand, a minimal induction of immune-related genes from H5N1-infected mice at early time point may have caused delayed IFN-γ expression, possibly contributing to the induction of the rapid replication of virus. Additionally, interferon γ -inducible protein (CXCL10), an important chemoattractant for T lymphocytes and monocytes, was more highly upregulated in the H5N1 group at 3 and 7 dpi, and this may have led to more robust immune activation, inducing persistent activation and recruitment of immune cells (e.g. T cells, macrophages, and monocytes). In line with this, our previous in vivo results suggest that the cytokine storm-like phenotype is the hallmark of severely lethal HPAI A(H5N1) virus infection in mice 7 .
Dual RNA-seq also revealed a good correlation between viral transcript levels and host gene expression levels. Accordingly, viral transcript counts of NS1 for H5N8 were also found to be higher than those for H5N1 at 1 dpi, suggesting a strong correlation between sequenced numbers of viral transcripts and expression levels of host immune-related genes in each virus-infected lung at 1 dpi. More H5N1 genes than H5N8 genes were sequenced at 3 and 7 dpi, which is consistent with the higher viral titres in infected lungs and the higher virulence observed in vivo. The increase in the replication of H5N1virus in lungs may have resulted in the induction of immune cell recruitment with robust production of chemokines (such as CXCL10 and CXCL11); in contrast, H5N8 replication was relatively moderate and thus caused a relatively moderate induction of these chemokines at both 3 and 7 dpi. This observation may reflect a positive correlation with host responses related to virus sensing and virus-mediated inflammation, such as TREM1 signalling, acute phase response signalling, and virus receptor signalling, as shown by the IPA heat map analysis (Fig. 3B). Collectively, the activated innate immune system might have directly caused more severe pathological injury to the lung tissue, as demonstrated by increased activation of death receptor signalling in H5N1-infected mice.
The present study is the first to characterise the dual host and virus transcriptome of H5N8-infected animals and compare it with that of recent H5N1 infection. Differential expression of specific factors observed between avian influenza virus strains could explain the variation in disease pathogenicity. These findings provide a framework for future studies examining the molecular mechanisms underlying the pathogenicity of H5N8. Taken together, our data provide a comprehensive transcriptome analysis of H5N8-infected mice, and several biological pathways related to anti-viral innate immunity and inflammation responses are involved in the highest regulated DEGs. Further global immunity studies on H5N8-vaccinated animals will provide valuable insight into understanding the immunogenicity of a novel H5N8 vaccine.

Methods
Ethics statement. All animal experiment protocols performed in this study followed general animal care guidelines mandated under the Guidelines for Animal Use and Care of the Korea Center for Disease Control (KCDC). They were approved by the Laboratory Animal Research Center (approval No. CBNUA-074-0904-01), which is a member of the Institutional Animal Care and Use Committee of Chungbuk National University, and were performed in a level 3 biosafety laboratory (BSL3) in Chungbuk National University (permit No. CBNU-BSL-14-001). The methods were carried out in accordance with the approved guidelines.

Virus infection in mice.
Six-week-old female BALB/c mice (Samtako, Seoul, Korea) were used for all infection experiments. Ten mice per group were inoculated intranasally with 10 2 -10 7 EID 50 /50 mL of MDk/W452(H5N8) and EM/W149(H5N1). Body weights and survival patterns were recorded daily for 14 days. For virological and pathological examinations, groups of mice were inoculated intranasally with 10 5 EID 50 /50 mL of virus, and six mice per group were euthanised at 1, 3, 5, 7 and 9 dpi to examine the growth kinetics of the virus in mouse lungs.
RNA quality check, library construction, and sequencing. RNA quality was assessed by analysis of rRNA band integrity on an Agilent RNA 6000 Nano Kit (Agilent Technologies, Santa Clara, CA). Prior to cDNA library construction, poly (A) mRNA was enriched using 2 μ g of total RNA and magnetic beads with oligo(dT) primer. The purified mRNAs were then disrupted into short fragments, and double-stranded cDNAs were immediately synthesised. The cDNAs were subjected to end-repair, poly (A) addition, and connection with sequencing adapters using the TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA). Suitable fragments, automatically purified using a BluePippin 2% agarose gel cassette (Sage Science, Beverly, MA), were selected as templates for polymerase chain reaction (PCR) amplification. The final library sizes and qualities were evaluated electrophoretically using an Agilent High Sensitivity DNA Kit (Agilent Technologies, Santa Clara, CA), and the fragment was found to be between 350 and 450 bp. Subsequently, the library was sequenced using an Illumina HiSeq 2500 sequencer (Illumina, San Diego, CA). Three biological replicates from RNA-seq run were performed per time point. RNA-seq data from this study have been submitted to the NCBI Sequence Read Archive (SRA) under the accession number of PRJNA290088(SUB1025618) (http://www.ncbi.nlm.nih.gov/sra/). Transcriptome data analysis. Low-quality reads were filtered according to the following criteria: (i) reads containing more than 10% of skipped bases (marked as 'N's); (ii) reads containing more than 40% of bases with quality scores less than 20; and (iii) reads with average quality scores of each read of less than 20. The whole filtering process was performed using in-house scripts. Filtered reads were mapped to the mouse reference genome (Ensembl release 77) 28 using the aligner STAR v.2.3.0e 29 . Gene expression level was measured with Cufflinks v2.1.1 30 using the gene annotation database of Ensembl release 77. Non-coding gene regions were removed with the mask option. To improve the accuracy of measurement, multi-read-correction and frag bias-correct options were applied. All other options were set to default values. For differential expression analysis, gene-level count data were generated using the HTSeq-count v0.5.4p3 31 tool with the options "-m intersection-nonempty" and "-r option considering paired-end sequence. " Based on the calculated read count data, DEGs were identified using the R package called TCC 32 . The TCC package applies robust normalisation strategies to compare tag count data. Normalisation factors were calculated using the iterative DEGES/edgeR method. The q-value was calculated based on the p-value using the p adjust function of the R package with default parameter settings. DEGs were identified based on a false discovery rate (FDR) q-value threshold of less than 0.05. Genes were considered differentially expressed with displaying a change of more than two-fold with a q-value of less than or equal to 0.05 for at least one time point.
Gene ontology (GO) enrichment analysis. The GO database classifies genes according to the three categories of biological process, cellular component, and molecular function, and predicts the function of the selected genes. To characterise the identified genes from DEG analysis, a GO-based trend test was performed using Fisher's exact test; p-values < 0.001 were considered statistically significant.