miR-21 is upregulated, promoting fibrosis and blocking G2/M in irradiated rat cardiac fibroblasts

Experiment design

We performed WGCNA and DEG analysis of gene expression in X-ray-irradiated CFs. Because the degree of influence of irradiation on cells varies with the X-ray dose (Zhao et al., 2019), we used two doses: 1 Gy and 8.5 Gy. The cells were divided into a control group (un-irradiated, 0 Gy), 1 Gy group, and 8.5 Gy group. Additionally, to investigate which genes were regulated by miR-21 after radiation, interference of miR-21 in CFs was performed. Groups were defined as the anti-miR-21 group, anti-miR-21+1-Gy group, and anti-miR-21+8.5-Gy group. These six groups of cells were sequenced.

Each group was subjected to WGCNA and DEG analysis separately. The identified DEGs were compared with the module eigengenes (MEs) from WGCNA, and consensus genes from the MEs with DEGs (—log2FC— ≥ 1, adjusted p-value ≤ 0.05) were selected. Each group of consensus genes is referred to as a cluster; forinstance, in a consensus process, the same gene between DEGs in the 8.5-Gy group compared to 0-Gy group and in the brown module was denoted as the 8.5 vs 0-Gy cluster (Supplemental Information 1). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were used to analyze the functions of these genes.

Furthermore, the role of miR-21 after radiation was investigated; we chose new identified consensus genes (new-pick genes) in the target cluster showing significant down-regulation, and the expression level was attenuated in the anti-21+X-ray cluster according to log2FC (X-ray cluster)-log2FC (anti-miR21+X-ray cluster) >1 (Supplemental Information 1).

Culture and isolation of cardiac fibroblasts (CFs) and cardiomyocytes (CMs)

Healthy new born Wistar rats (days 0–2) purchased from Laboratory animal Center of Gansu University of Chinese Medicine. Parent rats were raised in the Laboratory center in a comfortable room with a 12 h dark/12 h light cycle and free access to food and water. Every sacrificed newborn rat comes from the same parents above. Cultivation of CFs and cardiomyocytes (CM) from new-born rats were cultured as described previously with some modifications (Thum et al., 2008). Briefly, rats were sacrificed after isoflurane (2%) inhalation and cervical dislocation. Primary cardiac fibroblasts and cardiomyocytes were separated from 4-6 mice cardiac tissue each time. The collected cells were plated in MEM-F12 containing vitamin B12, NaHCO₃, and 10% fetal calf serum. The cultures contained mostly primary cardiac fibroblasts, as >95% of cells were stained with antibodies for fibroblast-specific antigen prolyl-4-hydroxylase (P4HB, ab137110, Abcam, Cambridge, UK) and Anti-Vimentin antibody (ab92547, Abcam, Cambridge, UK), and >95% of cells were negative for the cardiomyocyte-specific marker α2-actinin (ACTN2, clone EA-53, Sigma) (Figs. 1A–1I). We cared for the animals and sacrificed the rats in strict accordance with animal welfare laws and regulations and animal welfare ethics requirements. This study was approved by the Ethical Committee (the certificate number: GZY2018-115, Date: 03-05-2018) of Gansu University of Chinese Medicine in Lan Zhou, China.

RNA isolation and gene expression

Total RNA containing small RNA was extracted using the miRVanaTM RNA Isolation Kit AM1561 (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNAs were transcribed to cDNA using the SuperScript™ III Reverse Transcriptase kit (Invitrogen).Stem-loop quantitative real-time PCR (sqRT-PCR) assay was performed to validate the miRNA, Q-PCR was conducted to detect mRNA expression levels using the following conditions: 95 °C for 5 min, 35 cycles of 95 °C for 15 s, 60 °C for 30 or 60 s, and 72 °C for 20 s. The relative expression level of miRNA/mRNA was normalized to U6/glyceraldehyde 3-phosphate dehydrogenase and fold-change in the miRNA/mRNA level was calculated with the 2-^△△Ct method. Each sample was analyzed in triplicate. All primers were obtained from RiboBio (Guangzhou RiboBio Co., Guangzhou, China).

Assay for transfection, apoptosis, and cell cycle

Rat CFs were transfected with antagomir of miR-21 (100nM, RiboBio) using HiPerFect Transfection Reagent (Cat No./ID:301705, Qiagen, Hilden, Germany). The percentage of CFs undergoing apoptosis was determined by staining with annexin V (Annexin V: FITC Apoptosis Detection Kit I, 556547, BD Biosciences, Franklin Lakes, NJ, USA) followed by fluorescence-activated cell sorting analysis. The cell cycle were evaluated by propidium iodide (PI) staining.

Ionizing radiation treatment

Irradiation of cell cultures containing 1 ×10⁶ log phase cells was performed with a precision X-RAY irradiator (X-RAD 225, Precision X-ray, North Branford, CT, USA) at a dose rate of 2 Gy/min. Doses of 1 Gy and 8.5 Gy were administered at room temperature, and control cells were sham-irradiated.

Illumina sequencing

Next-generation sequencing was performed with the Illumina NovaSeq6000 instrument (San Diego, CA, USA) at Beijing CapitalBio Technology Company (Beijing, China). Normalized gene expression values were calculated using the expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM) method (Mortazavi et al., 2008).

WGCNA construction

Gene hierarchical cluster network was constructed using the WGCNA R package (Langfelder & Horvath, 2008) and visualized in a dendrogram. All samples were clustered, and outlier samples (cutHeight = 3000) were excluded from data analysis (Fig. 2A). A β-soft power threshold of 10 was selected to ensure that the network satisfied a scale-free topology (R² > 0.9) based on the linear regression model fitting index obtained from the functions ’pickSoftThreshold’ operation (Figs. 2B–2C). Coexpression modules were detected using the function ‘blockwiseModules’ with default settings and modified parameters (power, 10; minModuleSize, 50; mergeCutHeight, 0.25). The grey modules indicate nonsense areas. Based on each module eigengene, correlation analysis was performed to identify modules that were significantly associated with the measured traits (apoptosis and cell cycle).

Identification of differentially expressed genes

‘DESeq2’ package in R was used to identify DEGs (Love, Huber & Anders, 2014). To control the false discovery rate, p-values were adjusted. Values of —log2FC—≥ 1 and padj ≤ 0.05 were considered as different.

Bio-function enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) (DennisJr et al., 2003) was used to perform GO function and KEGG pathway enrichment analysis. The GO terms of biological processes (BP), cellular components (CC), and molecular functions (MF) were assessed. A p-value <0.05 was considered to indicate a significant difference for the GO terms and KEGG pathways.

UpSet veen analysis

An UpSet veen map was constructed using online omicshare soft ( http://www.omicshare.com/tools/Home/Soft/getsoft).

Identification of protein–protein interaction networks

The online Search Tool for the Retrieval of Interacting Genes (STRING database, Version11.0; http://string-db.org/) was used to build a protein-protein interaction (PPI) network (Szklarczyk et al., 2019) (threshold of > 0.4).

Statistical analysis

All experimental data statistical analyses were carried out using GraphPad Prism software (version 7.0; GraphPad, Inc., La Jolla, CA, USA), and the data are expressed as the mean ± standard error of mean. Statistical differences between two groups were calculated by two-tailed Student’s t-test, and p < 0.05 was considered significant.

Upregulation of miR-21 in irradiated CFs

The sqRT-PCR results revealed upregulation of miR-21 in irradiated CFs, which persistently increased up to 7 days after 8.5 Gy irradiation (Fig. 1K); however, miR-21 expression did not increase with increasing radiation intensity. In addition, there were more apoptotic cells after irradiation, which increased with irradiation intensity, compared to non-irradiated control cells. Importantly, many more apoptotic cells were observed after interference of miR-21 with irradiation (Figs. 1L–1Q). Cell cycle analysis showed that X-rays can stall cells in the G2/M phase; however, this effect was strongly attenuated by treatment with the miR-21 inhibitor (Figs. 1R–1W). These results indicate that miR-21 has an important function in CFs.

Data processing and construction of weighted coexpression network and identification of key modules

To explore the co-expression patterns of CFs genes after irradiation, RNA-seq FPKM data from 18 samples were evaluated with the WGCNA package (Langfelder & Horvath, 2008) (Fig. 2). Results from the 8.5-Gy-3 sample were outliers and thus were excluded (Fig. 2A). We found 24 different MEs according to their degree of connectivity (Fig. 2D). The number of genes contained in each module is shown in Supplemental Information 2.

To analyze the expression of ME in each group, module genes of 17 samples were analyzed. We drew heat maps of gene expression (Fig. 3A). In the Brown module, gene expression was down-regulated after irradiation and severely down-regulated in the 8.5-Gy group but increased in the anti-21 and anti-21+1-Gy groups. In the 8.5-Gy+anti21 group, gene expression was down-regulated compared to in the anti-21 group but up-regulated compared to in the 8.5-Gy group. These results indicate that the Brown module has a central module relationship with miR-21. The Yellow and Turquoise modules showed the opposite gene expression pattern as the Brown module; after irradiation, gene up-regulation in the Turquoise module was intensified after interference with miR-21 (Fig. 3A, Figs. 4A–4C).

Features (apoptosis and cell cycle) relationship with modules and identification of key modules

To clarify the relationship between apoptosis and cell cycle with the modules, 24 different MEs were analyzed and a heat map was generated (Figs. 2D, 3E). Three modules showed notable results (Figs. 3B–3D). Brown module showed a significant negative correlation with G2/M phase, with a correlation index of −0.96 (p = 2 ×10⁻⁹); Yellow module, associated with late apoptosis, exhibited a correlation index of 0.80 (p = 1 ×10⁻⁴); and Turquoise module which is related to cell death and S phase showed a correlation index of 0.86 (p = 9 ×10⁻⁶; p = 1 ×10⁻⁵).

Correlation between modules and identification of key modules

Among the 24 different MEs, interaction associations were analyzed, and a tree map was constructed (Fig. 3F). The results showed that modules were independent of each other.

DEGs and consensus DEGs in target modules

To identify the key genes in MEs, we performed DEG analysis and evaluated consensus genes showing differential expression within different modules. DEGs (Supplemental Information 3) in the 1-Gy, 8.5-Gy, 1-Gy+anti21, and 8.5-Gy+anti21 (vs 0-Gy) groups were used to identify genes within the three WGCNA modules (Brown, Yellow, and Turquoise) (Figs. 4A–4C). After each module and DEGs were screened and filtered, the consensus genes were used to draw an UpSet veen map (Figs. 4D–4F). Subsequently, we selected 1-Gy vs 0-Gy cluster genes and 8.5-Gy vs 0-Gy cluster genes from the three modules for GO and KEGG analysis.

GO functional analysis of X-ray irradiated CFs

In the consensus study, we performed GO enrichment analysis of the consensus genes following evaluation with the DAVID web-based search tool. The enrichment is shown in Table 1 and (Figs. 5–7).

In the Brown module, among BP, the cluster of 8.5Gy vs 0Gy was significantly associated with cell division, including chromosome segregation and mitotic nuclear division, and with DNA replication, including DNA replication and DNA replication initiation (Fig. 5A). In CC enrichment analysis, the terms nucleus, nucleoplasm, kinetochore, and cytoplasm were enriched (Fig. 5B). In MF analysis, the 8.5Gy vs 0Gy clusters were significantly associated with biomacromolecule binding, including chromatin binding, ATP binding, nucleotide binding, poly(A) RNA binding, and so on (Fig. 5C). However, in the 1Gy vs 0Gy clusters, only the nuclear envelope integral membrane protein 1 gene was up-regulated, and thus GO and KEGG analysis were not performed.

For genes in the Yellow module, among the BP, the 1Gy vs 0Gy cluster was significantly associated with the cellular response to organic substance, skeletal muscle cell differentiation, response to insulin, and negative regulation of apoptotic process (Fig. 6A). In the cluster of 8.5Gy vs 0Gy, the same results were observed as in the 1Gy vs 0Gy cluster, which were associated with cellular response to organic substance, negative regulation of apoptotic process, and mainly enriched in muscle contraction (Fig. 6E). In CC enrichment analysis, terms related to not only the nucleus but also the extracellular space were enriched in the 1Gy vs 0Gy cluster (Fig. 6B). The extracellular matrix was also enriched in the 8.5Gy vs 0Gy cluster, which is closely related to cellular fibrosis (Fig. 6F). In MF analysis, the 1Gy vs 0Gy cluster was significantly associated with transcriptional activity, RNA polymerase II core promoter proximal region sequence-specific binding, and transcription regulatory region DNA binding (Fig. 6C), whereas the 8.5Gy vs 0Gy cluster was also significantly associated with transcriptional activity, RNA polymerase II core promoter proximal region sequence-specific binding, and protein heterodimerization activity (Fig. 6G).

For genes in the Turquoise module, among the BP, the 8.5Gy vs 0Gy cluster was mainly enriched in cellular response to interferon-gamma and chemokine-mediated signalling pathway (Fig. 7A). In CC enrichment analysis, terms related to the extracellular space were enriched in 8.5Gy vs 0Gy (Fig. 7B). In MF analysis, our results showed that 8.5Gy vs 0Gy were significantly associated with chemokine activity and CCR2 chemokine receptor binding (Fig. 7C). However, in the 1Gy vs 0Gy cluster, only the BTG anti-proliferation factor 2 gene was up-regulated and chromosome 14 open reading frame 132 (D430019H16Rik) was down-regulated.

KEGG pathway analysis

KEGG pathway maps of biological functions were obtained from the DAVID web-based search tool (Huangda, Sherman & Lempicki, 2009). In WGCNA, the Yellow and Brown modules showed a greater correlation with apoptosis, cell cycle, and cellular fibrosis compared to the other modules. We performed KEGG analysis to explore differences in the pathways in each cluster, See Table 1 and (Figs. 5D, 6D, 6H, 7D).

In the Brown module, the 8.5Gy vs 0Gy cluster was mainly enriched in the cell cycle, DNA replication, spliceosome, and DNA repair (Fig. 5D). In the Yellow module, the 8.5Gy vs 0Gy cluster was significantly associated with the P53 signalling pathway, HTLV-1 infection, viral carcinogenesis, and MAPK signalling pathway (Fig. 6H), and the 1Gy vs 0Gy cluster was significantly associated with HTLV-1 infection (Fig. 6D). In the Turquoise module, the 8.5Gy vs 0Gy cluster was significantly associated with the chemokine signalling pathway (Fig. 7D).

GO functional analysis and KEGG analysis of CFs after anti-miR-21 treatment and irradiation

We observed the response of cells to radiation after interfering with miR-21. As shown in Figs. 5E–5H in the Brown module, after 8.5 Gy irradiation (compared to 0 Gy), some genes were down-regulated; whereas, after 8.5Gy irradiation and interference miR-21 (compared to 0 Gy), only a slight decreased in the genes involved in the former group was observed. We selected these genes for GO and KEGG analysis (Table 1). BP analysis showed that functional genes were mainly concentrated in microtubule-based movement, chromosome segregation, and cell division, and were involved in regulation of the G2/M transition in the mitotic cell cycle, which agrees with our cell cycle results. CC analysis showed that functional genes were mainly centralized in the midbody, chromosome, centromeric region, and cytoplasm. MF analysis demonstrated that functional genes were mainly concentrated in microtubule binding, protein kinase binding, microtubule motor activity, and ATP binding. KEGG analysis clearly revealed that these genes were associated with the cell cycle.

PPI network analysis

To further investigate the function of the cluster genes at the protein level, we constructed a PPI network of the clusters (Figs. 8A–8E). Genes and GO terms of clusters of three modules; see Table 2 and Figs. 8A–8E. We further explored which genes and biological processes were being regulated by miR-21 elevation post-irradiation. Our study suggested that compared to the control (0 Gy), the Brown module gene was mainly down-regulated in irradiated cells, and there was a higher probability that the Brown module contained the target gene of miR-21 (Table 3).

Experimental validations

For the preliminary verification of the fibrosis of CFs after irradiation, the expression levels of Grem1, Clu, Gdf15, Ccl7, and Cxcl1 were determined. The Q-PCR results showed that the mRNA expression levels of Grem1, Clu, Gdf15, Ccl7, and Cxcl1 were significantly increased in the 8.5-Gy group in contrast to the controls (0 Gy) (Fig. 8F), showing 2.300 ± 0.073-, 1.630 ± 0.159-, 2.616 ± 0.149-, 3.225 ± 0.138-, and 2.115 ± 0.085-fold higher expression, respectively. We also validated the genes predicted online (Gdf15, Rsad2) as miR-21 targets (Fig. 8G). The results showed that the expression of these two genes was increased by 2.616 ± 0.149- and 2.115 ± 0.085-fold after irradiation, and increased by 5.207 ± 0.111- and 49.853 ± 0.058-fold after interference miR-21 and following radiation (P < 0.05).