Transcriptome Analysis of the Regulation of Retinoic Acid on the Growth and Development of Brandt’s Vole

Brandt’s vole is distributed throughout Russia, Mongolia, Inner Mongolia and China. And is one of the most harmful rats to grassland ecology, particularly in the temperate steppe of China. Brandt’s voles not only overgraze the forage grass, resulting in a substantial decrease in the livestock carrying capacity of the region with high population levels, but also accelerate vegetative degradation and desertification. For example, after vole caves are abandoned, the ecological environment needs 3-5 years to recover, which seriously affects the development of animal husbandry (Wang et al., 2010). Integrated control of Brandt’s voles is a complex animal husbandry ecological project. At present, the primary control measures that are practiced include artificial capture and stocking of natural enemies, such as snakes and cats and the use of chemicals such as anticoagulant rodenticide (Hou et al., 1991). In order to improve the prevention and control measures, we aimed to uncover different strategies that might help assist in the management of Brandt’s voles from the perspective of limiting biological growth and development. Toward this goal, we explored the mechanisms that regulate the growth and development of Brandt’s voles.

Retinoids are important biological molecules that act in cell proliferation, epithelial cell growth and maturation, apoptosis, immunological functions and are vital in establishing organ development during embryogenesis (Niederreither et al., 2008). Lactase (lactase-phlorizin hydrolase, or LPH) is a member of the β-galactosidase family of enzymes and is involved in the hydrolysis of the disaccharide lactose into constituent galactose and glucose monomers (Day et al., 2000). Quercetin has antioxidant and anti-inflammatory effects which might help reduce inflammation, kill cancer cells, control blood sugar and help prevent heart disease (Ali et al., 2018; Kim et al., 2018; Kim et al., 2019). These three chemical molecules all play a regulatory role in apoptosis in the growth and development of mammals.

The purpose of this study was to find the most significant treatment for controlling Brandt’s vole populations by comparing the mortality and body weight changes of Brandt’s voles under various chemical treatments. Furthermore, we screened differentially expressed genes related to the growth, development and death of mice under this treatment to explore a new mechanism for regulating the growth and development of Brandt’s vole.

Sample preparation
 
To study the function of retinoic acid in female B. vole development, we selected female mice for test (Nowickyj et al., 2008). 270 B. voles, 108 of them were female voles, were captured in the Naranbulag Soum, Abagage, Xilingol League, Inner Mongolia, China, from where Brandt’s vole populations began expanding, from June 12, 2018 to July 08, 2018.

Female B.Voles divided into 35 cages, each cage for 3 to 4 voles as a group. B. Voles are stored at 23±1°C temperature room. The size of a small iron tank containing B. Voles is (27 x 28 x 22 cm) with sawdust for comfortable environment. Every day prepare food (carrot) and water twice (morning and evening).

LPH inhibitor, retinoic acid, quercetin and DMSO (Standard treatment concentration 10 mg/ml) were given them through oral with the help of Animal oral gavage feeding needles and 1 ml syringe was used. Depending on bodyweight, each vole treated with a standard of 0.5 ml/30 g. The doses were given them in the interval of 24 hours, 72 hours and 120 hours. Each treatment has five repeats, each repeat with 5 voles. After each treatment, the hind legs muscles of surviving Brandt’s voles were randomly sampled and mixed. All the samples were frozen at -80°C for transcriptomic analysis.
 
RNA extraction, library construction and sequencing
 
Total RNA was isolated from all four treatments using the RNA plant mini kit with on-column DNase digestion (Qiagen, Hilden, Germany) following the manufacturer’s instructions. RNA degradation and contamination were detected on 1% agarose gels. RNA concentration was measured using a Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, Carlsbad, CA, USA). Additionally, RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

A total amount of 3 μg RNA per sample was used as input material for cDNA preparations. Finally, three samples with RNA integrity number (RIN) values>8 were used for construction of cDNA libraries. Sequencing libraries were generated using the NEBNext Ultra™ RNA Library Prep Kit for Illumina (NEB, USA) following manufacturer’s instructions.

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and 125 bp paired-end reads were generated.
 
Sequence reads, mapping, assembly and annotation
 
Raw data (raw reads) of fastq format were firstly processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapters, reads containing ploy-N and low-quality reads from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality. Transcriptome assembly was accomplished based on left.fq and right.fq using Trinity with min-kmer-cov set to 2 by default and all other parameters set default (Hao et al., 2017) .

Gene function was annotated based on following seven databases: Nr (NCBI non-redundant Protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (Protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (A manually annotated and reviewed protein sequence database), KO (KEGG Ortholog database) and GO (Gene Ontology), using Diamond software (Benjamin et al., 2015).
 
DEGs analysis
 
Prior to DEG analysis, for each sequenced library the read counts were adjusted using the edge R program package through one scaling normalized factor. DEG analysis of two samples was performed using the DEGseq R package (Anders et al., 2010). P value was adjusted using Padj value (Storey et al., 2003). Padj value£0.05 and |log 2 (fold change) ³1 was set as the threshold for significantly differential expression.
 
KEGG enrichment analysis of DEGs
 
KEGG is a database resource for understanding high-level functions and utilities of the biological systems (Kanehisa et al., 2008) at the cell, the organism, or the ecosystem level. We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways (Mao et al., 2005).

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, cry and per, 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).
 
Mortality statistics
 
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_al1990). 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).





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
 
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
 
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).





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
 
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).

DEGs of Brandt’s voles treated with retinoic acid were enriched in multiple KEGG signaling pathways including: Ko04512:ECM-receptor interaction, ko04140:Regulation of autophagy, ko00051:Fructose and mannose metabolism, ko04974:Protein digestion and absorption, ko04977:Vitamin digestion and absorption, ko04151:PI3K-Akt signaling pathway, ko04911:Insulin secretion, ko04713:Circadian entrainment, ko04215:Apoptosis-multiple species, ko04068:FoxO signaling pathway, ko04210:Apoptosis, ko03320:PPAR signaling pathway, ko04910:Insulin signaling pathway, related to birth, growth and rearing. Enrichment of the ko04620: Toll-like receptor signaling pathway revealed that retinoic acid can be involved in the regulation of growth and development of Brandt’s vole by regulating these networks.

None