Regulation of all-trans retinoic acid in various animals
High concentrations of all-trans retinoic acid (atRA) inhibit adipogenic differentiation of rat bone marrow mesenchymal stem cells. The mechanism of retinoic acid nuclear receptor ¡2 (PPAR ¡2) regulates peroxisome proliferation-activated receptor ¡2 (PPAR ¡2) and CCAAT enhancer binding protein a (C/EBP a) during adipogenic differentiation of rat bone marrow mesenchymal stem cells (Zhang et al., 2013)
. Retinoic acid, as a transcriptional regulatory factor, has a wide range of physiological functions, it regulates cell proliferation and differentiation at the physiological level, which is closely related to normal growth and development; as a teratogen, retinoic acid can induce a variety of tissue and organ deformities in a variety of organisms (Feng et al., 2018).
There is a feedback regulation cycle of migratory locust rhythm genes clk
, which are all regulated by all-trans retinoic acid and promote migratory locust diapause. Upstream all-trans retinoic acid, rarp
and rhythm genes promote the expression of downstream dh
, which in turn regulates the insulin transactivation pathway (Babar et al., 2019; Hao et al., 2019).
Retinoic acid, the first morphogen described so far in vertebrates, is a vitamin A derivative which exerts striking effects on development and differentiation. The identification of three retinoic acid receptors as members of the nuclear receptor super-family provides an explanation for the molecular action of morphogens on gene expression (de Théet_al
1990). RA levels require precise regulation by controlled synthesis and catabolism and when RA concentrations deviate from normal levels, in either direction, abnormal growth and development occurs (McCaffery et al., 2003).
In rodents, the control of neural patterning and differentiation are disrupted when RA concentrations are lowered, whereas inappropriately high concentrations of RA result in abnormal development of the cerebellum and hindbrain nuclei. The latter parallels the malformations seen in the human embryo exposed to RA due to treatment of the mother with the acne drug Accutane (13-cis RA). In cases where the child survives beyond birth, a particular set of behavioral have been described. Even adult brains might be susceptible to an imbalance of RA, particularly the hippocampus (McCaffery et al., 2003).
According to the mortality, it could be seen that the mortality rate of B. voles
after treatment with retinoic acid was 81%, which was significantly higher than that of DMSO, which was only 32%. Thus, we took retinoic acid-treated and DMSO-treated Brandt voles’ hind leg muscles for transcriptome sequencing.
Illumina HiSeq mRNA sequence
In total, 494.6 million short reads were generated from the six libraries (RA and DMSO), with 481.4 million classified as clean and at least 94% high-quality (Q30) reads selected for further analysis (Table 1). After de novo
assembly was completed, we aligned unigene sequences to protein databases, including NR, Swiss-Prot, KEGG, KOG, eggNOG, GO and Pfam (e-value<1e-5) via
Diamond software (Benjamin et al., 2015).
We then retrieved proteins with the highest sequence similarity to the given transcripts along with their functional annotations. Of the 48844 unigenes, we found 21513 that (44%) that were annotated (Table 2).
Table 1: Statistics of mortality and vole body weight change under four treatments. We used the average of surviving vole post-treatment body weight minus pre-treatment body weight as body weight change.
Table 2: Summary statistics of clean reads in the muscle transcriptomes of voles.
The UniGene database obtained by splicing was used as a database and the expression abundance of each Unigene in each sample was identified by sequence similarity alignment. 86.79% ~ 87.89% of clean reads from six samples were compared to the Unigene database by using the bowtie2 software (Langmead et al., 2012)
(Table 2) and the eXpress software (Roberts et al., 2013)
to calculate the Unigene expression FPKM value. The correlation coefficient were above 0.99, indicating high quality of the experimental replicates (Fig 1).
Differentially expressed genes (DEGs) between DMSO and Retinoic acid treatments
Fig 1: Correlation between the transcriptomes of two treatment samples. (a) Correlation coefficient analysis of RNA-seq data from six samples. (b) PCA plot showing clustering of transcriptomes of DMSO-fed and retinoic acid-fed Brandt’s voles.
A total of 4047 and 2251 up- and down-regulated transcripts were observed in the retinoic acid treatment compared to the DMSO (control) treatment (FDR≤0.001 and |log2 ratio|≥1), respectively (Fig 2).
GO and KEGG pathway enrichment analysis
Fig 2: Differential gene expression in the DMSO treatment as compared to the retinoic acid treatment. (a) Number of up-regulated (upper bars) and down-regulated (lower bars) genes examined in the comparison. (b) Clustering Heatmap of the total 6298 DEGs. Red bands indicate up-regulated genes while blue bands indicate down-regulated genes.
We analyzed the KEGG enrichment of the 6298 DEGs and found that the DEGs were enriched in 225 different signaling pathways including: retinol metabolism (ko00830), the Foxo signaling pathway (ko 04068) and Apoptosis-multiple species (ko 04215), along with 31 metabolic pathways were significantly enriched (P-value<0.05) including: ECM-receptor interaction (ko 04512), protein digestion and absorption (ko 04974) and the PI3K-Akt signaling pathway (ko 04151) (Table 3). At the same time, we analyzed 1231 genes from 6289 DEGs annotated to the GO database by GO enrichment. A total of 4861 assigned terms were obtained, of which 1845 terms were significantly enriched. We found 847 GO terms with the corresponding number of DEGs greater than two among the three GO classifications (Molecular function, cellular component and biological process) (Fig 3).
Table 3: Significantly enriched KEGG pathways comparing retinoic acid treatment with DMSO treatment.
Fig 3: GO enrichment analysis top 30. screening GO entries with a number of different DEGs greater than 2 in the three classifications, according to the -log10 P value of each entry, 10 items in descending order. Caption: The X coordinate in the figure is the name of the GO entry and the Y coordinate is -log10 P value.
Combining the analyses of functional enrichment from GO, KEGG and phenotypic data, we selected 114 genes related to cell proliferation, cell division, regulation of biological growth and development and regulation of cell death related metabolic pathways, such as cyclic AMP-responsive element-binding protein 3-like protein 1, gamma-aminobutyric acid receptor-associated protein-like 1, phosphorylase b kinase gamma catalytic chain, skeletal muscle/heart isoform and SHC-transforming protein 4 isoform X1 (Fig 4).
Validation by RT-PCR
Fig 4: 114 candidate DEGs selected based on GO and KEGG enrichment results. Red in the picture indicates that the expression of the corresponding gene is up-regulated and blue indicates that the expression of the corresponding gene is down-regulated.
DEGs of c-Jun N-terminal kinase and tubulin epsilon chain, were used for verify the RNA sequencing results. Our results showed that the relative expression of c-Jun N-terminal kinase and tubulin epsilon chain was upregulated, which were similarly to what was found in the RNA sequencing results (Fig 5).
Fig 5: The relative quantitative expression trends of two specific DEGs.