RNA sequencing reveals novel LncRNA/mRNAs co-expression network associated with puerarin-mediated inhibition of cardiac hypertrophy in mice

Study design and transverse aortic constriction (TAC)-induced myocardial hypertrophy mouse model

Male C57BL/6 mice (aged 5 weeks old) used in this study were purchased from the Guangdong Medical Laboratory Animal Center (Guangzhou, China). The mice were fed a standard laboratory diet and housed under standard conditions. Two weeks later, a total of 15 mice were randomly assigned using computer-generated random number tables to sham group, TAC group and TAC+PUE group. The sham group mice received intraperitoneal (i.p.) injection of the saline and the surgical procedures without the constriction. TAC group mice received i.p. injection of the saline and the TAC procedures to induce cardiac hypertrophy. The TAC+PUE group mice received a daily i.p. injection of puerarin (Sigma Aldrich, St. Louis, MO, USA) 3 days before TAC and continued for 21 days after TAC. The doses administered for puerarin (100 mg/kg per day) were chosen based on previous studies (Chen et al., 2014).

TAC is widely used as a disease model of chronic ventricular hypertrophy (Bosch et al., 2020). Thereby, TAC was used to establish myocardial hypertrophy in mice as described in a previous study (Meng et al., 2018). Briefly, following anesthetization with isoflurane (Henry Schein, Melville, NY, USA), the aortic arch was exposed after a midline incision in the anterior neck. The transverse aortic arch between the left common carotid artery and the brachiocephalic artery was chosen as the site of constriction. The aortic arch was constricted by tying a 6.0 nylon suture ligature against a 26-gauge needle. After rapid removal of the needle, an incomplete constriction was formed. The successful constriction of TAC was verified using trans-thoracic echocardiography. All procedures were carried out according to the Guide for the Care and Use of Laboratory Animals (National Research Council (U.S.) et al., 2011). All animal experiments were approved by the Institutional Animal Care and Use Committee of the Second Affiliated Hospital of Guangzhou Medical University, Guangzhou, China (No. B2019-064).

Organ weight

Body weight (BW) and tibia length (TL) of each mouse were measured 21 days after the TAC procedure. The mice were sacrificed using cervical dislocation under aesthesia, and hearts were arrested in diastole with injection of potassium chloride. And then, the hearts were gained and heart weight (HW) were measured. Heart weight to body weight ratio (HW/BW) and heart weight to tibia length ratios (HW/TL) were counted.

Echocardiography

After anesthetization with isoflurane, successful ligation was confirmed by color Doppler and pulsed-wave Doppler scanning. Left ventricular posterior wall dimension (LVPWd) and interventricular end-diastolic septum thickness (IVSd) were measured by two-dimensional transthoracic echocardiography. The transthoracic echocardiography was performed by an experienced technician who was blinded to the study groups using an IE33 echocardiographic system (Philips Medical Systems, Leiden, the Netherlands).

Hematoxylin-eosin (HE) staining

After fixation with 10% formalin, the heart tissues were dehydrated through a serial alcohol gradient and embedded in paraffin wax blocks. And then, heart sections were stained with HE solution (Beyotime, China). The sections were examined under a light microscope (Nikon Technology Co., Ltd., Japan).

Western blot analysis

Western blot analysis was performed as described previously (Lin et al., 2018, 2021). Briefly, the whole protein of the heart was extracted using radio-immunoprecipitation assay (RIPA) lysis buffer (Kaiji Company, Shen Zhen, China) with protease and phosphatase inhibitors (Kaiji Company, Shen Zhen, China). Protein samples were separated by sodium dodecylsulphate polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to polyvinylidene difluoride (PVDF) membranes. The membranes were then incubated at 4 °C overnight with the following primary antibodies: β-MHC (Bioworld Technology Inc., Louis Park, MN; USA; BS70815, 1:1,000) and GAPDH (Cell Signaling Technology Inc., Danvers, MA, USA; # 5174, 1:1,000). After being washed three times with PBST, the membranes were incubated with the secondary antibodies for 1 h at room temperature. Then, the signals were detected using the Imaging System (GE, Amersham Imager 600, GE, Piscataway, NJ, USA). The relative expression level of proteins was analyzed using Image-Pro Plus 6.0 (Media Cybernetics Inc., Bethesda, MD, USA).

RNA extraction

Three mice in each group were randomly selected to harvest the heart for the RNA sequence analysis. The mice were sacrificed at 21 days after TAC and hearts were harvested. Total RNA was isolated from heart tissues using Trizol reagent (TaKaRa, Tokyo, Japan). The concentration of the RNA samples was evaluated using a NanoDrop ND-1000 instrument (Thermo Fisher Scientific, Waltham, MA, USA). The integrity of the RNA was assessed by electrophoresis on an agarose gel.

RNA sequencing analysis

Library preparation and Illumina sequencing analysis were performed as previous study (Li et al., 2017). Briefly, the RNA of heart tissue was used for library construction by KAPA Stranded RNA-seq library Prep Kit (Illumina, NEB, San Diego, CA, USA). Then, the library was pair-end sequenced by Illumina NovaSeq 6000 Sequencing system (Illumina, NEB, San Diego, CA, USA) according to the manufacturer’s protocol. The whole-transcriptome sequencing experiment was completed by Kangcheng Biotechnology Co., Ltd. (Shanghai, China). The per base quality scores indicating the sequences data is of high quality. After trimming the adaptor sequences, clean sequences reads were obtained for subsequently analysis. Then the clean reads were mapped to a reference genome using Hisat2 software. The transcripts were assembled and quantitated by StringTie. The expression level of transcripts was measured by transcripts per kilobase million (TPM) values using the R package of Ballgown (Meng et al., 2021; Zhao et al., 2021b). The unknown transcripts were further screened with Coding Potential Assessing Tool (CPAT) to distinguish the protein-coding genes from the non-coding genes. The unknown non-coding transcripts were then used to screen for novel lncRNAs. The unknown non-coding transcripts with lengths more than 200 nt and more than two exons were selected as novel lncRNA candidates (Zhao et al., 2021a). Raw sequence files have been deposited at NCBI’s Gene Expression Omnibus (Accession code: GSE176244). Differentially expressed lncRNAs with statistical significance among the three groups were observed through volcano plot filtering. Hierarchical clustering was conducted to demonstrate the distinguishable lncRNAs expression pattern among the groups. The volcano plot filtering and hierarchical clustering analysis was conducted from the TPM values using the R package of ggthemes and pheatmap accordingly.

Quantitative real-time polymerase chain reaction (qRT-PCR) analysis

The lncRNA-seq results were further validated by qRT-PCR analysis as previous study (Lin et al., 2020). Briefly, total RNA isolated from heart tissues was reverse transcribed to synthesize cDNA using SuperScriptTM III Reverse Transcriptase kit (Invitrogen: 18080-044). Then, the qRT-PCR was performed by the QuantStudio5 Real-time PCR System (Applied Biosystems, Waltham, MA, USA) with the 2 × PCR Master Mix (Arraystar: AS-MR-006-5). Airn (antisense of Igf2r non-protein coding RNA) expression was determined using the 5′-CTAAGCCACTGCCGTCCATA-3′ forward primer and 5′-TAGGCGTTCACCCTTGGTCT-3′ reverse primer. Relative expression levels of Airn were normalized with GAPDH (Forward primer: 5′-CACTGAGCAAGAGAGGCCCTAT-3′, Reverse primer: 5′-GCAGCGAACTTTATTGATGGTATT-3′). The results were calculated using the 2^–ΔΔCt method. All groups had three independent samples and all of the samples were conducted in triplicate. Cq values obtained from the RT-qPCR testing were considered acceptable for technical triplicates if they fell within mean Cq ± 2 standard deviations (SD).

Construction of the lncRNA-mRNA co-expression network

The coding/non-coding gene co-expression (CNC) network profile was constructed according to validated altered lncRNA and its related mRNAs. The CNC was established through the weighted gene co-expression network analysis (WGCNA) using the R package WGCNA (v1.69) (Langfelder & Horvath, 2008). The co-expression network between mRNA and lncRNA was constructed utilizing Cytoscape version 2.8.1 software (The Cytoscape Consortium, San Diego, CA, USA) based on the Pearson correlation analysis of lncRNA and mRNA (r > 0.95 or r < −0.95).

Gene function analysis

To explore the function of the selected genes in the co-expression network, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the targets genes were implemented using Database for Annotation, Visualization and Integrated Discovery (DAVID) (Lin et al., 2016).

Statistical analysis

All results were presented as the mean ± SD. Statistical analysis was performed using GraphPad Prism 7.0 software (San Diego, CA, USA). One-way ANOVA followed by Tukey’s post hoc analysis was used to test the differences among groups. A P value less than 0.05 was considered statistically significant.

Puerarin inhibited TAC induced-cardiac hypertrophy in mics

The ratio of heart weight to body weight or tibia length are useful indexes to quantify cardiac hypertrophy (Yin et al., 1982). As shown in Fig. 1A, the TAC group had significant high heart weight/body weight (HW/BW) ratio and heart weight/tibial length (HW/TL) ratio compared with sham group. Puerarin treatment resulted in a significant reduction of HW/BW and HW/TL ratio compared with TAC group. Echocardiography demonstrated that TAC resulted in the increase of LVPWd and IVSd, which were significantly reversed by puerarin (Fig. 1B). Besides, HE staining analysis of cardiac sections showed that the increased myocyte area in TAC mice was minimized by treatment with puerarin (Fig. 1C). Beta-myosin heavy chain (β-MHC) is a major component of cardiac myosin. The protein expression of β-MHC has historically been used as a marker of cardiac hypertrophy (Ding et al., 2019). Protein levels of β-MHC were detected by western blotting in this study. As expected, the protein expression of β-MHC increased dramatically in TAC mouse ventricular myocytes. Moreover, puerarin significantly reduced the β-MHC expression compared with TAC group (Fig. 1D).

lncRNA expression profiles and validation

RNA-seq was used to assess the expression levels of lncRNAs in heart samples of the sham, TAC and TAC+PUE groups. Overall, 32,105 lncRNAs were identified by RNA-seq among the groups. It was found that 23 lncRNAs were up-regulated and 716 were down-regulated in the TAC group compared with sham group (|log2 Fold change| > 2 and FDR < 0.05) through a volcano map (Fig. 2A). In addition, nine lncRNAs were up-regulated and one were down-regulated in TAC+PUE group compared with TAC group (|log2 Fold change| > 2 and FDR < 0.05) (Fig. 2B). Then, heat map showed that four lncRNAs downregulated in the TAC group were up-regulated by puerarin treatment (Fig. 3).

To validate the RNA-seq results, the lncRNA with most significant FDR value from the differentially expressed lncRNA profile was selected for the qRT-PCR analysis. The qRT-PCR results showed that there was significant difference among the groups in lncRNA Airn (P < 0.05) (Fig. 4).

lncRNA-mRNA network analysis

To investigate the underlying regulating mechanisms of Airn in cardiac hypertrophy and the therapeutic target of puerarin, we analyzed the co-expression of Airn and protein-coding genes by WGCNA. Pearson’s correlation coefficients of co-expression levels suggested that there was an interactive relationship between Airn and 2,387 mRNAs (r > 0.95 or r < −0.95). Those co-expressed genes including some important coding genes (Srf, Akt, Igf2, etc.) which have been reported to be involved in the pathogenesis of cardiac hypertrophy. We used Airn and the top 150 mRNAs co-expression with Airn to construct a lncRNA-mRNA visualization network using Cytoscape (Fig. 5).

Next, we carried out a functional enrichment analysis of the top 150 target genes using over-representation analysis. GO analysis found that these target mRNAs were enriched in some important biological processes and molecular function such as translational initiation, cell proliferation, insulin-like growth factor binding and poly(A) RNA binding (Figs. 6A–6C). KEGG analyses suggested that those Airn-interacted mRNAs were enriched in endocytosis, signaling pathways regulating pluripotency of stem cells and Jak-STAT pathway (P < 0.05) (Fig. 6D).