Sequence analysis, assembly, and functional annotation
The transcriptome sequencing results of brain tissue are shown in Table 1. After filtering the raw reads, 143,936,974, 142,787,930, 139,867,950 and 138,606,658 clean reads were obtained for the four brain states, respectively, and the corresponding base numbers were 22,640,739,600 bp, 22,425,822,300 bp, 22051944000 bp and 21,898,530,000bp (Table 1). The proportions of clean/raw reads of the four states were between 94.49% and 95.51%, suggesting the high quality of the RNA-Seq data available for further analyses.
Table 1. Sequencing and assembly statistics of brain samples for V. sinensis
The clean reads of this brain area and liver were assembled together, and a total of 403,707 unigenes sequences were obtained. The unigenes sequences were used as the reference genome in the subsequent analysis of the brain. The lengths of 403,707 unigenes ranged from 201 bp to 19,402 bp, with an average length of 511 bp, an N50 value of 615 bp, and an N90 value of 244 bp.
Among the 403,707 unigenes, 83,237, 60,588, 23,820 and 9,176 were respectively annotated by the NR, Swiss-Prot, KOG, and KEGG databases (Supplementary Fig. S1). For KOG annotation in particular, the term of signal transduction mechanisms was the most highly represented (Supplementary Fig. S2).
Degs Analysis
A correlation analysis was conducted based on the TPM values. From the correlation results (Supplementary Fig. S3), the sample had good repeatability and could be used for gene differential expression analysis.
To clarify the differences in the transcription levels of genes in different activity states of V. sinensis, we compared the transcriptome data of two adjacent time points and obtained four comparisons of DEGs. In addition, we also focused on the analysis of the rest vs awake (before and after feeding) and sleep vs activity comparisons, and a total of six pairwise comparisons were obtained (Table 2). In the six pairwise comparisons, a total of 672 DEGs were detected. In the rest vs activity comparison, the number of DEGs was the highest, with 136 genes, and the lowest number of DEGs was detected in the awake vs active comparison, with 60 genes. The numbers of upregulated genes in the comparison groups were higher than the numbers of downregulated genes in the two comparisons. In terms of the definitions of upregulation and downregulation of specific genes in each pairwise comparison, an upregulated gene was defined as one that had a higher expression level in the former time point than in the latter, and vice versa for a down-regulated gene. In the awake vs active comparison, this may be because bats are always active after awakening to prepare for subsequent hunting behavior followed by a highly active predation state, so the numbers of DEGs in the two states were less, as similar genes were expressed. A heatmap of the hierarchical clustering of all DEGs indicated the large differences in gene expression between different time points (Fig. 1).
Table 2. The number of DEGs in the six pairwise comparisons of brain tissue
Comparable group
|
Up
|
Down
|
Total
|
Rest vs sleep
|
63
|
59
|
122
|
Sleep vs wake
|
65
|
69
|
134
|
Wake vs activity
|
32
|
28
|
60
|
Rest vs activity
|
63
|
73
|
136
|
Rest vs wake
|
63
|
57
|
120
|
Sleep vs activity
|
47
|
53
|
100
|
Go And Kegg Analysis Of Degs
The six pairwise comparisons were analyzed by GO and KEGG enrichment, and the results suggested that the bats had significant differences in the metabolic and transport activities at different time points during the day (Supplementary Tables S1–S4). The rest (after feeding) vs awake (before feeding) comparison represented the two extreme states before and after predation. There were 63 genes that showed high expression levels in the after feeding state and low expression levels in the before feeding state (Table 2). In the cell component of GO, multiple complexes were significantly enriched by these 63 genes, including Prp19 complex, spliceosomal complex, and DNA packaging complex. In the category of molecular function, these genes were significantly enriched in terms of transferase activity, for example, histone methyltransferase activity (H4-R3 specificity), aryl sulfotransferase activity, and protein-arginine omega-N asymmetric methyltransferase activity. In the category of biological process, these genes were significantly enriched in hemoglobin synthesis-related terms such as positive regulation of hemoglobin biosynthetic process, regulation of hemoglobin biosynthetic process, and hemoglobin biosynthetic process (Fig. 2). Regarding the KEGG results, various metabolic pathways were significantly enriched by these 63 genes, including drug metabolism-cytochrome P450, metabolism of xenobiotics by cytochrome P450, and glutathione metabolism (Fig. 3). The above results indicated that hemoglobin synthesis, transferase activity, and metabolic activity of V. sinensis in the after feeding state may be higher than during the before feeding state.
There were 57 genes that showed low expression levels in the after feeding state and high expression levels in the before feeding state. In the cell component category, these genes were significantly enriched in the terms of proton transport. The three top GO terms with the most significant p values were mitochondrial proton-transporting ATP synthase complex, proton-transporting ATP synthase complex, and proton-transporting two-sector ATPase complex, catalytic domain. In the molecular function category, these genes were significantly enriched in terms that were closely related to various enzyme transport pathways such as cation transporting ATPase activity, beta mannosidase activity, and ATPase coupled ion transmembrane transporter activity. For the biological process category, transport-related terms were significantly enriched by these 57 genes, including hydrogen ion transmembrane transport, hydrogen transport, and proton transport (Fig. 2). Regarding the KEGG results, a variety of endocrine-related pathways were significantly enriched by these 57 genes, including the renin-angiotensin system, renin secretion, and cortisol synthesis and secretion (Fig. 3). The above results indicated that the intensity of transport and endocrine activities of bats in the before feeding state may be higher than in the after feeding state.
In addition, sleep and activity were the two states that best reflected the difference between the day and night activities of V. sinensis. There were 47 genes that showed high expression levels in the sleep state and low expression levels in the active state (Table 2). In the cell component category, these genes were significantly enriched in organelle-related terms such as cytoplasmic part, Golgi apparatus, and organelle membrane. In the molecular function category, these genes were significantly enriched in terms related to enzyme activities, including inositol tetrakisphosphate phosphatase activity, cytochrome c oxidase activity, and heme-copper terminal oxidase activity. In the biological process category, these 47 genes were significantly enriched in nuclear-related terms including nuclear matrix organization, nuclear matrix anchoring at nuclear membrane, and cytoskeletal anchoring at nuclear membrane (Fig. 2). For the KEGG results, the 47 genes in the sleep state were only significantly enriched in the ribosome pathway (Fig. 3).
In addition, there were 53 genes that were expressed at low levels in the sleep state and highly expressed in the activity state. In the cell component category, terms that were related to proton transport were significantly enriched by the 53 genes, including mitochondrial proton-transporting ATP synthase, catalytic core, and proton transporting ATP synthase. In the molecular function category, genes were significantly enriched in the terms related to DNA binding, including RNA polymerase II transcription factor activity (ligand-activated sequence-specific DNA binding), transcription factor activity (direct ligand regulated sequence-specific DNA binding), and chromatin DNA binding. In the biological process category, these genes were significantly enriched in multiple terms related to endocrine activities, and the three top GO terms with the most significant p values were response to corticotropin releasing hormone, cellular response to corticotropin releasing hormone stimulus, and regulation of type B pancreatic cell proliferation. In addition, these genes were also significantly enriched in terms related to the circadian clock, including entrainment of circadian clock, circadian regulation of gene expression, and circadian rhythm (Fig. 2). Regarding the KEGG results, various hormone secretion pathways and signal transduction pathways were significantly enriched by the 53 genes, including cortisol synthesis and secretion, MAPK signaling pathway, and PI3K-Akt signaling pathway. In addition, there were two circadian clock-related pathways that were significantly enriched, namely, circadian rhythm and the circadian rhythm-fly. Importantly, both circadian clock-related pathways contain the Per2 gene (Fig. 3). The above results showed that the two states of bats with completely opposite circadian rhythms do involve circadian clock-related genes and regulate the regular changes of some circadian physiological processes.
Trend Analysis Of Degs
To clarify the periodic expression and oscillation changes of genes in the brain area of V. sinensis during a 24-hour period, we examined the intersection of the six groups of DEGs, i.e., screened out genes that were expressed at all time points, to obtain a total of 557 genes. Then, we conducted trend analysis on these 557 genes, and the results are shown in Supplementary Fig. S4. Among the 20 modules obtained, the DEGs were significantly clustered in module 6 that contained 27 genes. These 27 genes had similar expression trends at the four time points during the day. The expression levels of these genes were decreased in the resting state, were the lowest in the sleep state, and gradually increased in the waking state, reaching the highest expression level in the activity state (Fig. 4a). These genes in module 6 were used in enrichment analysis, and the results are shown in Fig. 4 and Supplementary Table S5. Numerous circadian clock-related terms and pathways were significantly enriched by these 27 genes. Through the above analysis, we believe that the change trend of the genes in module 6 was consistent with the biological habits of bats. These 27 genes showed obvious circadian rhythms and were closely associated with bat activity states. Bats in an active state require high expression levels of these genes, so the physiological and biochemical processes that these genes participate in play important roles in the activity and periodic awakening of bats.
The important physiological and biochemical processes enriched by genes of module 6 are shown in Fig. 4. The GO results for 27 genes in module 6 suggested that a total of 161 terms were significantly enriched; 238 terms were significantly enriched in the biological process category; 26 terms were significantly enriched in the molecular function category, and 18 terms were significantly enriched in the cell component category (Supplementary Table S5). In the cell component category, the top three terms with the most significant p values were transcription factor complex, varicosity, and nuclear membrane. In the molecular function category, terms related to DNA molecular binding were significantly enriched, and the top three terms with the most significant p values were RNA polymerase II transcription factor activity, transcription factor activity, and steroid hormone receptor activity. In the biological process category, terms related to endocrine activities were significantly enriched by genes of module 6, and the top three GO terms with the most significant p values were response to corticotropin releasing hormone, cellular response to corticotropin releasing hormone stimulus, and regulation of type B pancreatic cell proliferation. Additionally, terms related to circadian clocks were also significantly enriched, including entrainment of circadian clock and the entrainment of circadian clock by photoperiod (Fig. 4b). The KEGG results of module 6 are shown in Fig. 4. The top three pathways with the most significant p values were cortisol synthesis and secretion, aldosterone synthesis and secretion, and Cushing syndrome. The circadian rhythm and circadian rhythm-fly pathways were also significantly enriched (Fig. 4c). The above results suggested that there were circadian clock-related genes in module 6 that have important molecular functions in the circadian rhythm.
Important Circadian Clock Genes
Combining the results of DEGs and trend analysis, we found that the genes in the sleep vs activity comparison and the genes in module 6 were enriched in pathways related to the circadian clock (Supplementary Table S3 and Table S5). We further analyzed the genes in the pathways related to the circadian clock and found that they all involved Per2. In the GO results of the sleep vs activity comparison and trend module 6, Per2 was significantly enriched in terms related to the circadian clock. Therefore, we believe that Per2 may be an important circadian clock gene in V. sinensis. We marked the Per2 gene in the figure illustrating the circadian rhythm pathway to show its important role in the pathway (Supplementary Fig. S5).