Integrated analysis of microRNA and mRNA expression profiles in rats with selenium deficiency and identification of associated miRNA-mRNA network

Selenium deficiency is closely related with various type of cardiovascular disease. However, the miRNA-mRNA regulatory network in Selenium deficiency related cardiac change remains to be understand. In the present study, a reliable Selenium deficiency rat model was established and confirmed by pathological and biochemical examination. The mRNA and miRNA expression profiles were conducted by microarray technology. Gene Ontology (GO) Analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Analysis was performed to investigate the function of targeted genes, and the relationship between miRNA and mRNA was studied by network analysis. A total of 4931 mRNAs and 119 miRNAs was differentially expressed between any two groups (control group, low-selenium group and selenium supplementation group). GO and KEGG pathway analysis of selected miRNAs target genes found that selenium deficiency was related to several different biological processes. Furthermore, a miRNA-mRNA regulatory network was conducted to illustrate the interaction of miRNAs and these targeted genes. In conclusion, our present study provides a new insight that potential molecular mechanism of Selenium deficiency was a multiply miRNAs and mRNA caused biological change.

Selenium concentration detection. The selenium concentration was assessed by the flameless atomic absorption spectrophotometry method using a Z-5000 spectrophotometer (Hitachi, Ltd., Tokyo, Japan) with a cathode lamp of Se (resonance line, 196.0 nm; Photron, Victoria, Australia). Standard selenium solution was used to calibrate the results.
Histology. The hearts were removed and rinsed in phosphate-buffered saline, fixed in 4% paraformaldehyde for 24 h, embedded in paraffin and cross-sectioned into 10-µm slices. Tissue sections were stained with hematoxylin/eosin and observed by light microscopy.
Microarray analysis of miRNA and mRNA. The gene microarray analysis of the 3 groups was completed, total RNA from each sample was extracted by using Trizol reagents (Invitrogen) and the miRNeasy mini kit (Qiagen). RNA quality and quantity was measured by using nanodrop spectrophotometer (ND-1000, Nanodrop Technologies, Wilmington, DE, USA) and RNA integrity was determined by gel electrophoresis. the miRCURY ™ Hy3 ™ /Hy5 ™ Power labeling kit (Exiqon, Vedbaek, Denmark) was used according to the manufacturer's guideline for miRNA labeling. After purification with the miRNeasy mini kit (Qiagen), the Hy3TM-labeled samples were hybridized on the miRCURYTM LNA Array (v. 16.0) (Exiqon) according to the array instructions. total RNA from each sample was linearly amplified and labeled with Cy3-UTP. The labeled cRNAs were purified by RNeasy Mini Kit (Qiagen). The labeled cRNAs were hybridized onto the Whole Genome Oligo Array (4 × 44 K, Agilent Technologies). Microarray was then washed by Wash buffer kit (Exiqon) and scanned by the Axon GenePix 4000B microarray scanner (Axon Instruments, Foster City, CA). the GenePix Pro 6.0 software (Axon Instruments) was utilized to extract the data.
Differentially expressed miRNAs and mRNAs definition. Expressed data was normalized by using the Median normalization method. In this study, fold change (FC) method was utilized to identify differential mRNA and miRNA between these three groups, with following criteria: (1) Fold Change >2 or Fold Change <0.5; and (2) FDR <0.05.
Series Test of Cluster. Series Test of Cluster (STC) algorithm was utilized to analyze the gene expression dynamics and observe the changes of gene expression under different situations 13,14 . The raw expression values were converted into log2ratio. Using a strategy for clustering the different expression tendencies of differentially expressed mRNAs and miRNAs in each group, we established 8 unique expression profiles, respectively. The expression model profiles were related to the actual or the expected number of genes assigned to each model profile. Significant profiles have higher probability than expected by Fisher's exact test and multiple comparison test.
Prediction of miRNA-targeted Gene and construction of miRNA-mRNA network. The miRNA target prediction tools, TargetScan and miRnada, were utilized to further explore the targeted mRNAs, which were regulated by differentially expressed miRNAs 15 . Then We selected mRNAs from the two expected profiles, namely profile NO. 1 and profile NO. 6. These genes closely correlate with the results of regulatory gene expression, including activation and inhibition in the context of heart damage induced by Selenium deficiency. Furthermore, combining the differentially expressed miRNAs and mRNAs as well as the predicted targets for these miRNAs, a core miRNA-mRNA regulatory network was constructed by using the Cytocape software. In the miRNA-mRNA network, the circle represents mRNA and the shape of arrow represents miRNA, and their relationship was represented by one edge. The center of the network was represented by degree. Degree means the contribution one miRNA to the genes around or the contribution one gene to the miRNAs around. The key miRNA and mRNA in the network always have the biggest degrees.
SCienTifiC RepoRts | (2018) 8:6601 | DOI:10.1038/s41598-018-24826-w GO Analysis and KEGG pathway Analysis. GO analysis was applied to analyze the main function of the specific genes with significant differences in the representative profiles of miRNA target genes. Fisher's exact test was applied to determine the significant GO categories and P values were corrected by FDR. Only GOs that had a P-value of < 0.001 and an FDR of < 0.05 were chosen. Enrichment provides a measure of the significance of the function: as the enrichment increases, the corresponding function is more specific, which helps us to find those GOs with more concrete function description in the experiment.
Pathway analysis was used to find out the significant pathway of the differential genes. Pathway annotations of Microarray genes were download from KEGG (http://www.genome.jp/kegg/). Fisher exact test was also used to find the significant enrichment pathway. The resulting P values were adjusted using the BH FDR algorithm 16 . Pathway categories with p-value < 0.05 were reported.
Procedures of constructing the miRNA-mRNA network associated with selenium deficiency were shown in Fig. 1.

Statistical analysis.
All the results were expressed as mean ± standard deviation. Statistical analysis was performed with one-way analysis of variance for multiple comparisons. All the statistical analyses were performed using the software SPSS 18.0 for Windows (PASW Statistics, SPSS Inc., Chicago, IL). P < 0.05 was considered to indicate a statistically significant difference. Fig. 2.

Selenium concentration in blood samples. The Se concentration in blood was demonstrated in
As shown, compared with the control group, the Se concentration decreased significantly in the LS and SS groups (P < 0.05). After Se supplementation by intraperitoneal injection, Se concentration in blood in the SS group was increased in contrast with LS group (P < 0.05).  Histological characterization of the heart. HE staining revealed that histological changes of the rat hearts caused by selenium deficiency (Fig. 3). In the control group, well-connected intercalated discs and normally arranged filaments were seen clearly. Swelling and disorganized muscle fibers were shown in the LS group. Compared with the control group, the structure of myocardial tissue was mildly swelling in the SS group, but their arrangements were ordered relatively.

Microarray analysis and Significant Differential Gene Analysis.
To identify the possible gene expressional change in selenium deficiency rat models, mRNA and miRNA microarray was performed by using the Affymetrix probe dataset. A total of 4931 mRNAs and 119 miRNAs was differentially expressed with more than 2.0-fold change between any two groups, as shown in heat map analysis (Fig. 4). Of 4931 differentially expressed mRNAs, the expression of 2514 genes were upregulated and 2417 genes were downregulated. Of 119 differentially expressed miRNAs, the expression of 65 miRNAs were upregulated and 54 miRNAs were downregulated.

Series Test of Cluster of miRNA and mRNA expression profiles.
To further characterize mRNA and miRNA expression profiles, we establish 8 expression profiles for the 119 differentially expressed miRNAs and 4931 mRNAs by using STC method in STEM software. Each profile indicates an expression model including a set of miRNAs or mRNAs with analogous expression pattern, respectively. As shown in Figs 5, 2 expression patterns of mRNA (Profile NO. 2 and NO. 5) and 2 expression patterns of miRNA (Profile NO. 1 and NO. 6) showed statistically significance (p < 0.05), respectively. Among these patterns, we chosen the 2 significant miRNA expression patterns, which were accorded with the pattern that we expected. As we known, miRNA leading to downregulation of target gene expression at transcriptional or translational levels 17 , A negative relationship between miRNAs and their targeted genes has been reported in the literature. Then relevant mRNA expression profiles (Profile NO. 1 and NO. 6) were also chosen to further analysis. Profile NO. 1 contained 30 miRNAs and 308 mRNAs whose expression rapidly declined from LS group to control group but stable at SS group. This group of miRNA might be involved in the heart pathological changes of selenium deficiency. while Profile NO. 6 contained 10 miRNAs and 397 mRNAs whose expression remained temporally stable from SS group to control group but decreased at SS group. To the 30 miRNAs of profile NO. 1, we found 36738 target genes, to the 10 miRNAs of profile NO. 6, we found 9612 target genes, respectively.

GO Analysis and KEGG pathway Analysis of putative target genes.
With the aim of exploring the relationship between differentially expressed genes and Selenium deficiency induced heart damage, GO and KEGG pathway analysis were applied respectively. GO analysis of representative profiles of miRNA target genes were shown in Fig. 6. on the biological process level, GO terms were mainly related to metabolism and development associated process, such as lipid metabolic process and heart development. On the molecular function level, GO terms were mainly involved in protein binding and enzyme activity. On the cellular component level, GO terms were mainly included cell junction, cell projection and cytoplasm. KEGG Pathway analysis evaluating the enriched pathways for the representative profiles of miRNA targeted genes showed that some of them were involved in signal transduction pathways (Fig. 7), such as focal adhesion kinase, Rap1 signaling pathway, PPAR signaling pathway, Ras signaling pathway and Wnt signaling pathway.
Analysis of the miRNA-mRNA network. Depend on these miRNAs (in the Profile NO. 1 and Profile NO. 6) and their targeted genes, the miRNA-mRNA network was created to summarize the interaction of miR-NAs and these targeted genes. The highest degrees always mean the key miRNA in the network. As shown in Fig. 8A, 30 downregulated miRNAs which were belonged to Profile NO. 1 were in the miRNA-mRNA network. While these miRNAs which had higher degree were in the center of the network, such as rno-miR-125a-5p is the one that regulates the most genes, as many as 41 predicted target genes. These decreased miRNAs may lead to the increased expression of their target genes, for instance, Tspan5, Lgr4 and Mmd2, to regulate related signal transduction. The network of 10 upregulated miRNAs of Profile NO. 6 was illustrated in Fig. 8B, including rno-mir-207, rno-mir-152, rno-mir-133a, may repress the expression of negative modulators, for example, Fads6, Samd14 and Fbxo46, respectively. To summarize all the results, we have displayed that these crucial miRNAs may function as important regulators by lipid metabolic process, heart development, protein binding, enzyme activity, cell junction, cell projection and signal transduction, thus playing pivotal roles in heart damage caused by selenium deficiency.

Discussion
As one of the essential trace elements of mammalian, selenium regulate various biological processes such as redox balance, immunity, aging and disease 18,19 . Selenium-deficiency can cause a variety of cardiovascular diseases, most of them are induced by oxidative damage and inflammation 20 . At present, the molecular mechanism of Selenium-deficiency remains to be understood. Due to being incorporated as the unique amino-acid selenocysteine into selenoproteins, selenium is essential for the synthesis, metabolism and function of hormones 1,21 . Therefore, a comprehensive understanding of the specific role of selenium in cardiovascular disease is meaning for the treat strategy. The integrated analysis of miRNA-mRNA regulatory network was utilized to help us improve the cognition of gene regulations on selenium deficiency related cardiac disease.
Of all the differentially expressed genes, a substantial proportion were related with lipid metabolic process, heart develop, cell adhesion, blood vessel development, cell differentiation, signal transduction and positive regulation of cardiac muscle hypertrophy. These differentially expressed mRNAs might take part in the selenium deficiency mediated biological process. We also acquired the associated signal pathways which were modulated by the differentially expressed mRNAs through KEGG analysis. Among the 44 pathways with p value < 0.05, 4 pathways were related to metabolism of Lipid, including glycerophospholipid metabolism, glycerolipid metabolism, sphingolipid metabolism and fatty acid metabolism. Since selenoproteins adjust apolipoprotein E levels and regulate the gene expression of cholesterol synthesis, metabolism, and transport, they are critical to lipoprotein metabolism 22 . Selenium deficiency caused cellular responses often lead to the induction of inflammation and oxidative overload involving multiple signaling pathways 10,23 . Several signaling pathways were mediated by the differentially expressed mRNAs, such as focal adhesion, Rap1 signaling pathway, PPAR signaling pathway, Ras signaling pathway, Oxytocin signaling pathway, MAPK signaling pathway and Wnt signaling pathway. As mentioned in the literature, Selenium deficiency can support oxidation of two oxidation sensors and activate Wnt pathway 24 . Combining the function and pathway analysis for differentially expressed genes in selenium deficiency rat model, we may have concluded that cardiac cell might develop different cellular responses and gain corresponding pathological changes due to deficiency of selenium.
Pertaining to the different miRNAs, we successfully identified 30 downregulated miRNAs and 10 upregulated miRNAs in our rat models. Of these miRNAs, such as rno-miR-16, rno-miR-199a-5p, rno-miR-195, rno-miR-30e, which had been reported to be involved in selenium deficiency induced heart failure 11 . Although several research teams have conducted the clinical research to clarify the relationship between selenium deficiency and cardiovascular disease [25][26][27][28] , few studies have paid attention to the important roles of miRNAs and miRNA-mRNA interaction of molecular mechanism.    The integrated analysis of the miRNA-mRNA regulatory network was a reliable strategy to enhance the understanding of gene regulation. Specific miRNA-mRNA interactions and high-throughput analyses of the downstream effects of miRNAs on mRNAs have been the focus of many research studies 29,30 , which provided valuable insights on regulation of gene expression. As to Network-based prediction were extensively applicated in various cancer to diagnosis, predict tumor clinical phenotypes and test personalized drug targets [31][32][33][34] . We can also develop an algorithm to predict these diseases which are associated with selenium deficiency at early stage, by using machine learning method and selenium deficiency related miRNA-mRNA network. Actually, selenium deficiency caused cardiovascular disease is a complex process, several potential mechanisms such as oxidative stress, inflammation have been validated 35,36 . However, the molecular mechanism of miRNAs regulates gene (mRNA) expression of selenium deficiency has not been studied. It is necessary to identify the target with opposite expression patterns of miRNA and select the pivotal miRNA as well as their target mRNAs in selenium deficiency related cell change. With the help of bioinformatics analysis, we filtered several miRNAs which were in the center of the miRNA-mRNA regulatory network and corresponding target genes, such as rno-miR-125a-5p. In addition, more work need to be conducted to reveal the concrete regulatory mechanism of selenium deficiency related cardiovascular disease.
In conclusion, this report reveals that a number of miRNAs and mRNAs that correlated with selenium deficiency caused cardiovascular disease, and also illustrated a core functional miRNA-mRNA network which was meaningful for reveal the regulation from miRNAs to mRNAs in the selenium deficient rat model.