All experimental protocols included in this manuscript were approved by The Institutional Animal Care and Use Committee of Xi'an Medical University. The NPRA conditional knockout mouse model (NprA^flox/flox) was generated by the Shanghai Model Organisms Centre, Inc. This model was developed using the CRISPR/Cas9 system against the C57BL/6 mouse background, as previously described (19,20). An NprA donor vector containing flox sites flanking exon 2 of the NprA gene was constructed. Two single guide RNA (sgRNA) sites targeting intron 1 and intron 2 were transcribed in vitro. The sequences of the sgRNA target sites were 5′-GAGACGAGGACCAAGACTGCAGG-3′ and 5′-GCACAGGGGCCGTCTCATGCAGG-3′. The donor vector with two sgRNAs and Cas9 mRNA was microinjected into C57BL/6 fertilized eggs. F₀-generation mice with positive homologous recombination were identified by long PCR. The primers used for correct 5′ homology arm recombination were 5′-GACTTCAAATCCCAGCTCACAACC-3′ and 5′-CAAAGATCGAATGGAGCCCTGTGT-3′. The primers used for correct 3′ homology arm recombination were 5′-TATCTTGCGGCCCTCTGACTGTAT-3′ and 5′-ACTGGCCCTCCGTGGTTAGCA-3′. The identity of the PCR products were confirmed by sequencing. NprA flox heterozygous mice were obtained by mating F₀ generation mice with wild-type C57BL/6 mice. The genotypes of F₁-generation NprA^flox/flox mice were identified using the same genotyping methods as F₀-generation mice.

To achieve tissue-specific deletion of the NPRA gene, NprA^flox/flox mice were crossed with Myh6-Cre transgenic mice (Shanghai Model Organisms Centre, Inc.) to obtain NprA^flox/flox/Myh6-Cre mice. Both the NPRA^−/− and NPRA^−/+ mice were identified by agarose gel electrophoresis and ethidium bromide staining (Fig. S1).

All mice were euthanised by inhaling oxygen with 5% isoflurane at a rate of 1 l/min. The mice were confirmed to be deeply anaesthetized if they were immobile for 1 min. A 25% volume of CO₂ gas was constantly flowed into the chamber to bring the mice to the point of clinical death, which was confirmed by cessation of the heartbeat and respiration. The hearts were subsequently removed to obtain the myocardium. Briefly, the cardiac tissue was isolated and placed in ice-cold PBS and all fat and connective tissue were removed. The samples were stored in −80°C until use.

To measure cardiac mRNA expression of NPRA in mice, the myocardium samples (n=5 in each group) were collected and washed in PBS. Total RNA was then isolated from myocardium samples using TRIzol^® (Thermo Fisher Scientific, Inc.) from myocardium samples. RNA concentration was then measured by the NanoDrop™ ND-1000 (NanoDrop™ Technologies; Thermo Fisher Scientific, Inc.). The RNA samples were reverse transcribed into cDNA using SuperScript™ III Reverse Transcriptase (Invitrogen; Thermo Fisher Scientific, Inc.). Briefly, the mixture contained RNA, Oligo (dT)₁₈ (Takara Bio, Inc.), Random 6 mers (Takara Bio, Inc.), dNTPs Mix (HyTest Ltd.), RNase free dH₂O and was placed in a 65°C water bath for 5 min and on ice for 2 min. Following a short centrifugation in a palm centrifuge for 5–10 sec to make the mixture to gather at the bottom of the centrifuge tube, the RT reaction solution including 5X First-Strand Buffer (Invitrogen; Thermo Fisher Scientific, Inc.), DTT (Beijing Solarbio Science & Technology Co., Ltd.), RNase Inhibitor (Epicentre; Illumina, Inc.) and SuperScript™ III Reverse Transcriptase (Invitrogen; Thermo Fisher Scientific, Inc.) was sequentially added and the mixture was kept at 37°C for 1 min. The solution was incubated at 50°C for 60 min, then the reaction was terminated by incubation at 70°C for 15 min. The cDNA obtained was diluted 1:10 and stored at −20°C until use. qPCR was performed using TransStart™ SYBR Green qPCR Supermix (Takara Bio, Inc.) to determine the expression levels of NPRA. β-actin was included as an internal normalization reference gene for mRNA expression. All the PCR products were subjected to electrophoresis on a 2% agarose gel and visualized using ethidium bromide staining to confirm the presence of a single band of the expected size. The amplification efficiencies of qPCR for the target gene and the endogenous control were approximately equal and were calculated through a dilution series of cDNA. The qPCR reaction system included 2 µl cDNA, 5 µl 2C Master mix (Arraystar, Inc.) and 0.5 µl of each primer (10 µM solution) and diethylpyrocarbonate-treated water was added for a total volume of 10 µl. All qPCR reactions were performed in triplicate, including template-free controls. The thermocycling conditions consisted of an initial denaturation at 94°C for 3 min, followed by 40 cycles of denaturation at 94°C for 10 sec, and annealing at 56°C 20 sec, then extension at 72°C for 30 sec, and a final extension at 72°C for 5 min The data were analysed using the comparative Cq (ΔΔCq) method (21).

To confirm the expression of differentially expressed genes (DEGs), including mRNAs and circRNAs, briefly, the sequences of candidate mRNAs and circRNAs were obtained from the CircInteractome database (http://circinteractome.nia.nih.gov/). The GAPDH and β-actin housekeeping genes were used as internal controls. The reverse transcription was performed as described above. The thermocycling conditions were as follows: i) Initial denaturation at 95°C for 10 min; and ii) 40 cycles of denaturation at 95°C for 10 sec, and annealing at 60°C for 1 min. A melting curve analysis was performed after the amplification was complete by cooling the reaction to 60°C and then heating slowly to 99°C. The amplification products were analysed by 2% agarose gel electrophoresis and ethidium bromide staining. The relative quantification of the expression of the investigated genes was performed using a standard curve.

All primers were designed using Primer 5.0 and synthesized using a ViiA 7 Real-time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). Primer sequences are presented in Tables I and II.

For the preparation of tissue lysates, hearts from mice were dissected and washed with ice-cold PBS, then homogenized with RIPA buffer (Beyotime Institute of Biotechnology) containing 150 mmol/l NaCl as previously described (22). Cell lysates were prepared for western blotting as previously described (23). Protein concentrations were determined using a BCA protein assay kit (Beyotime Institute of Biotechnology). After adjustment of protein concentration, the lysates were boiled in SDS loading buffer for 5 min then resolved by 10% SDS-PAGE and 20 µg of protein was loaded per lane. Gels were then transferred to an Immobilon-P PVDF membrane (EMD Millipore), which was blocked by 5% non-fat dried milk in TBST (Tris Buffer Saline with 0.1% Tween 20) buffer for 2 h at room temperature and incubated with primary antibodies against mouse NPRA (cat. no. LS-C264634; 1:500; LifeSpan BioSciences, Inc.) and GAPDH (cat. no. 2118, 1:1,000; Cell Signaling Technology, Inc.) in TBS + 0.1% Tween-20 (TBST) and 5% non-fat dry milk overnight at 4°C. Membranes were subsequently washed in 5% milk/TBST and incubated with HRP-conjugated secondary antibody for 2 h at room temperature. Protein bands were visualised using a chemiluminescence detection kit (EMD Millipore).

Total RNA from each sample was quantified by agarose gel electrophoresis and further verified and quantified by Nanodrop ND-1000 spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc.). Agarose electrophoresis was used to determine the integrity of total RNA samples. A total of 1–2 µg total RNA from each sample was used for RNA-seq library preparation. Briefly, rRNA was removed from the total RNA with a RiboZero Magnetic Gold kit (cat. no. MRZG12324; Epicentre; Illumina, Inc.). The rRNA-depleted samples were then prepared for sequencing using the KAPA-stranded RNA-Seq Library Prep kit for Illumina (cat. no. KK8400; Kapa Biosystems; Roche Diagnostics). The library preparation procedure included the following steps: i) RNA molecules fragmentation; ii) first strand cDNA synthesis by reverse transcription; iii) second-strand cDNA synthesis dUTP instead of dTTP (with no Uracil-DNA glycosylase treatment); iv) end repair and 3′ adenylation of the double-stranded cDNA; v) ligation of the Illumina adapter; vi) PCR amplification; and vii) purification using magnetic beads (cat. no. KK8000; Kapa Biosystems; Roche Diagnostics).

The quality of the libraries was evaluated using an Agilent 2100 Bioanalyzer, in order to check the concentration, the fragment size distribution (400–600 bp), and levels of adapter dimer contamination. The amount was quantified using the absolute qPCR method (24). For sequencing, the barcoded libraries were mixed in equimolar amounts. Sequencing was carried out using TruSeq SR Cluster kit v3-cBot-HS (150 cycles; paired end, cat. no. GD-401-3001; Illumina, Inc.) on an Illumina HiSeq 4000 platform according to the manufacturer's instructions. The sequencing data are available in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo) under accession no. GEO140678.

As previously reported (25), image analysis and base calling were performed using the Off-Line Base Caller software v1.8 (Illumina, Inc.). Sequence quality was examined using FastQC software v0.11.7 (26). The cutadapt software (v1.2.1; http://code.google.com/p/cutadapt/) was then used to trim the fragments with a 5′,3′-adaptor and to filter the reads ≤20 bp. The next step was alignment of the reads to the reference genome using Hisat 2 software (2.1.0; http://daehwankimlab.github.io/hisat2/) (27). The average reads was 35,871,134.2 and 88.9% of the clean reads were aligned. According to a previously published method (28), the transcript abundances were estimated with StringTie v1.3.3 (http://ccb.jhu.edu/software/stringtie/), and the number of identified genes and transcripts per group was calculated based on the mean of fragments per kilobase of exon model million reads mapped (FPKM) in each group, and the transcripts with FPKM ≥0.5 were retained. The ballgown R package v 3.5.0 (www.bioconductor.org/packages/release/bioc/html/ballgown.html) was used to analyse the differentially expressed genes (DEGs) and transcripts by comparing the NPRA^−/− group to the NPRA^+/+ group. The adjusted P-values of DEGs were obtained by multiple checking and correction using the Benjamini and Hochberg method (29). Filtering was based on an absolute fold change cut-off >1.5, adjusted P≤0.05 and mean FPKM ≥0.5. CPAT v1.2.4 (http://lilab.research.bcm.edu/cpat/) was then used to determine whether transcript sequences were coding or non-coding. Principal component analysis (PCA) and Pearson correlation analysis were carried out on gene expression levels. Hierarchical clustering, Gene Ontology (GO) (30), Kyoto Encyclopedia of Genes and Genomes (KEGG) (31) pathway, scatter plot and volcano plot analysis were performed with the differentially expressed genes in the R v3.5.0 (http://www.r-project.org/), Python v2.7 or shell environment for statistical computing and graphics.

Total RNA from the samples of both groups was quantified. As described previously (32,33), sample preparation and microarray hybridization were performed according to the manufacturer's protocols (Arraystar, Inc.). Briefly, total RNA from each sample was digested with RNase R (Epicentre, Inc.) to remove linear RNA. The resulting circRNA-enriched samples were then amplified and labelled with fluorescently labelled nucleotides using an Arraystar Super RNA labeling kit (Arraystar, Inc.) utilizing a random priming method according to Arraystar Super RNA labeling protocol. The labelled circRNA molecules were hybridized to the Arraystar Mouse circRNA Array v2 (8×15 K, Arraystar). Microarray slides were scanned with a G2505C Agilent Scanner (Agilent Technologies, Inc.).

The array images were analysed by Agilent Feature Extraction software (version 11.0.1.1; Agilent Technologies, Inc.). Quantile normalization and subsequent data processing were performed using with the limma package (version 3.11; http://bioconductor.org/packages/release/bioc/html/limma.html). Differentially expressed circRNA molecules (DE-circRNAs) with statistical significance in expression levels between the two groups were identified using a volcano plot. The circRNAs that were differentially expressed between the two groups were filtered by fold change filtering (FC>2). Hierarchical clustering was subsequently conducted to illustrate the distinguishable circRNA expression profiles among samples. The raw data of microarray data are available in the GEO (https://www.ncbi.nlm.nih.gov/geo/) under accession no. GSE140798.

As described previously (34,35), the circRNA-miRNA interaction was predicted using the Arraystar miRNA target prediction software (version 1.0, Arraystar Inc.) based on TargetScan44 and miRanda, and the DE-circRNAs were annotated in detail with circRNA-miRNA interaction information. TargetScan (version 7.2; http://www.targetscan.org/), miRanda (http://www.microrna.org/), and miRDB (http://www.mirdb.org/) were used to predict the microRNA-mRNA interactions that match the seed region of mouse miRNA sequences obtained from miRBase (http://www.mirbase.org/). In addition, if a miRNA molecule could bind to both the circRNA and target mRNA, the targeted circRNA was defined as a candidate competing endogenous RNA (ceRNA) for the gene, and the circRNA-miRNA-mRNA interactions thus represent candidate ceRNA pairs (circRNA connected with miRNA and miRNA connected with mRNA). In addition, circRNA and genes in candidate ceRNA pairs are required to have the same direction of change of expression level (upregulation or downregulation) in DEG analysis because of a positive correlation in the expression of ceRNA pairs (36). Then, the DE-circRNAs and differentially expressed mRNA transcripts (DE-mRNAs) that were validated by RT-qPCR were included within the ceRNA network, while all others were excluded. The circRNA-miRNA-mRNA interaction network was constructed and visually displayed using Cytoscape (version 3.8.1; http://www.cytoscape.org/) based on the data analysis results.

Statistical analysis was performed using R statistical package (37) (v3.5.0; http://www.R-project.org/), SPSS v19.0 (SPSS, Inc.) and Prism 6 (GraphPad Software). A total of five samples were included in each group. DEGs and DE-circRNAs were obtained by comparing the NPRA^−/− group and the NPRA^+/+ group. Student's t-test was used to obtain P-values. The adjusted P-values of DEGs were obtained by multiple checking and correction based on the Benjamini and Hochberg method (29). Filtering was based on absolute FC >1.5 in DEGs and FC >2 in DE-circRNAs, adjusted P-value ≤0.05 and mean FPKM ≥0.5 were employed for filtering. RT-qPCR, statistical significance was assessed with Student's t-test. Data are presented as the mean ± SEM of three experimental repeats. P<0.05 was considered to indicate a statistically significant difference.

To obtain a Cre recombinase-mediated, loxP-directed deletion of the NprA gene, a targeting construct was designed in which the second exon of NprA was flanked by two loxP sites (Fig. 1A). In NprA^flox/flox/Myh6-Cre mice the NPRA gene was specifically disrupted in the myocardium. Following Cre recombination, exon 2 of NprA was deleted, leading to gene reading frame transcoding and protein function loss. The loss of cardiac expression of NprA was determined by RT-qPCR (Fig. 1B) and western blotting (Fig. 1C). Homozygous knockout (NPRA^−/−) and heterozygous (NPRA^−/+) mice were identified by PCR genotyping and agarose gel electrophoresis staining (Fig. S1).

The transcriptional profiles of NPRA^−/− and NPRA^+/+ mice were determined using RNA-seq. The resulting scatter plot (Fig. 2A), volcano plot (Fig. 2B), and hierarchical clustering map (Fig. 2C) highlighted significant (DEGs) between NPRA^−/− and NPRA^+/+ mice (n=5 in each group). Scatter plots are commonly used to display differences in gene expression between two datasets and to assess trends in population distribution. The values corresponding to the × and y-axes in the scatter plot are the normalized, log₂-scaled FPKM of NPRA^+/+ and NPRA^−/− mice, respectively. Symbols above the upper green line and below the lower green line indicate fold change >1.5 between the between NPRA^−/− and NPRA^+/+ mice. The volcano plot displays the relationship between the fold change of each DEG and the associated P-value in each sample group. Data clusterisation successfully discriminated NPRA^−/− from NPRA^+/+ samples. All 10 samples were included for further analysis. However, only very few of these genes were considered to be differentially expressed with an adjusted P<0.05 and fold change >1.5. Indeed, seven were upregulated, and nine were downregulated.

GO enrichment analysis was used to identify significantly enriched GO terms associated with the potential functional roles of DEGs through three categories: Biological process (BP), molecular function (MF), and cellular component (CC). Enriched terms with significant differences (P<0.05) were selected and ranked by enrichment score (-log₁₀-scaled P-value). A total of 4 BP terms, 2 MF terms, and 11 CC terms were significantly enriched in NPRA^−/− mice, relative to NPRA^+/+ mice (Fig. S2A). Intriguingly, 2 out of 4 significantly enriched BP terms were related to the ‘lipid metabolic process’ (Fig. 3A). In KEGG pathway analysis, only a single significant pathway, ‘metabolic pathways’, was observed (Fig. 3B).

Among the downregulated genes, 99 BP terms, 13 MF terms, and 11 CC terms were found by GO analysis (Fig. S2B). The top 10 BP terms in the NPRA-deficient group ranked by enrichment score are shown in Fig. 3C. Notably, three of the top 10 BP terms were related to ‘metabolic process’. These three terms are regulation of cellular metabolic process (GO: 0031323), regulation of metabolic process (GO: 0019222) and negative regulation of metabolic process (GO: 0009892).

To validate the RNA-seq results and confirm the differential gene expression profile, a subset of 12 genes which involving with metabolic process terms we got in GO analysis, including Agpat12, Alg12, Cbfa2t2, Cdkn1a, Nutf2, Rhobtb1, S100a9, Tfpi, Zbtb4, Hp, Gbp10, and Gm20521, were selected out of all 16 DEGs, and their expressions were measured using RT-qPCR. Among these genes, Alg12, Cdkn1a, S100a9, Gbp10 and Gm20521 were significantly differentially expressed in NPRA^−/− mice, compared with NPRA^+/+ mice (Fig. 4). Thus, these five genes were selected for subsequent experiments.

circRNA can act as a miRNA ‘sponge’ that regulates the expression of mRNA, and such circRNA-miRNA-mRNA interaction network plays important roles in the genesis and development of cardiovascular diseases (38). Therefore, it was hypothesized that a circRNA-miRNA-mRNA network might be at play in the context of cardiac-specific NPRA deletion, thus inducing metabolic dysfunction in the myocardium.

A murine circRNA microarray was used to analyse the circRNA profiles of the myocardium from NPRA^−/− mice, compared with NPRA^+/+ mice. The normalized intensity for the NPRA^+/+ and NPRA^−/− groups was nearly the same (Fig. 5A). DE-circRNA molecules between NPRA^−/− and NPRA^+/+ mice are displayed in a hierarchical clustering map (Fig. 5D), as well as a scatter plot (Fig. 5B) and a volcano plot (Fig. 5C). Clusterisation successfully discriminated successfully discriminated NPRA^−/− from NPRA^+/+ samples. In total, 55 upregulated and 197 downregulated DE-circRNAs were identified.

To further examine the regulatory mechanism underlying DEG profiles, all global changes in the expression pattern of circRNAs were predicted according to the complementary miRNA matching sequence. Using ceRNA analysis, a circRNA-miRNA-mRNA interaction network was constructed. To confirm the ceRNA network, the predicted circRNAs were compared to the DE-circRNAs profiles from the circRNA microarray, and the circRNAs involved in both the predicted circRNA profiles and DE-circRNA profiles were selected for validation. In addition, because of a positive correlation in the expression of ceRNA pairs, candidate circRNAs were required to be oriented in the same direction direction of change of expression level in DEG analysis. Using this approach, 14 downregulated circRNAs (mmu_circRNA_005865, mmu_circRNA_19029, mmu_circRNA_19474, mmu_circRNA_19519, mmu_circRNA_22217, mmu_circRNA_25320, mmu_circRNA_26033, mmu_circRNA_29619, mmu_circRNA_30261, mmu_circRNA_32945, mmu_circRNA_36265, mmu_circRNA_36266, mmu_circRNA_42481 and mmu_circRNA_43449) were selected for validation by RT-qPCR. With RT-qPCR analysis, 6 circRNAs were significantly downregulated in the NPRA^−/− mice, compared with NPRA^+/+ mice (Fig. 6).

Thus, the DE-circRNAs and DEGs that have been validated by RT-qPCR assay were included within the ceRNA network, and all others were excluded. Thus, with the majority of possible pathways excluded, only three target mRNA molecules (Alg12, Cdkn1a, and Gbp10) were included within the ceRNA network. The entire circRNA-miRNA-mRNA interaction network was constructed and visually displayed by Cytoscape. The top five miRNA binding sites predicted by ceRNA analysis for the six validated circRNAs in the ceRNA network are presented in Table III. The top 5 predicted miRNA response elements (MREs) in the 3′ UTR for each validated circRNA are presented in Fig. S3.

A circRNA-miRNA-mRNA interaction network was constructed, which included three genes (Alg12, Cdkn1a, and Gbp10) and six circRNAs (mmu_circRNA30261, mmu_circRNA43449, mmu_circRNA36266, mmu_circRNA32945, mmu_circRNA42481, and mmu_circRNA29619). Moreover, 162 miRNAs were predicted in this network (Fig. 7; Table IV). This network may provide insight into the potential interactions between circRNA candidates and their target genes.