Differential lncRNA/mRNA expression profiling and ceRNA network analyses in amniotic fluid from foetuses with ventricular septal defects

Patient recruitment and ethics statement

In this study, pregnant women were recruited from the Second Affiliated Hospital of Fujian Medical University. All pregnant women had indications for prenatal diagnosis and were assessed by prenatal ultrasound and echocardiography. For safety and reasons of ethical constraints, amniotic fluid sampling was only performed when pregnant women had indications for prenatal diagnosis and without additional burden for the mother or fetus. Ultimately, 22 singleton pregnant women aged 19–34 years between 20 and 22+6 weeks of gestation were included in the research. Of these, 11 pregnant women who carried a VSD foetus were in the VSD group, and the other 11 pregnant women who carried a structurally normal foetus but had a high risk Down’s screening were in the control group. Both groups of pregnant women had no pregnancy complications and poor life histories such as hypertension, diabetes, smoking, or alcoholism. All cases were tested for chromosome anomalies by karyotype. Foetuses with abnormal karyotypes and extracardiac major associated anomalies were excluded. There were no differences in maternal age, the number of pregnancies, gestational age, or body mass index (BMI) between the VSD group and the control group (P > 0.05). These clinical characteristics are shown in Tables 1 and S1. This study was performed in accordance with the principles of the World Medical Association Declaration of Helsinki and licenced by the Ethics Committee of the Second Affiliated Hospital of Fujian Medical University (NO. 2021-073).

RNA extraction

AF was extracted by amniocentesis from pregnant women and specimens with gross maternal blood or meconium contamination were discarded. AF samples were centrifuged at 1,600×g for 10 min at 4 °C to remove amniocytes, and then supernatant samples were stored at −80 °C until RNA extraction. Total RNA was extracted from AF supernatant using TRIzol LS Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. Then, RNA quantity and quality were measured by a NanoDrop ND-1000.

RNA labelling and array hybridization

In this study, we selected Arraystar Human LncRNA Microarray V5.0, which is designed for the global profiling of human lncRNAs and protein-coding transcripts and can detect approximately 39,317 lncRNAs and 21,174 coding transcripts. The specific operations were as follows, each RNA sample was amplified and transcribed into fluorescent cRNA along the entire length of the transcripts without 3′ bias utilizing a random priming method (Arraystar Flash RNA Labelling Kit, Arraystar). The labelled cRNAs were purified using an RNeasy Mini Kit (Qiagen). The concentration and specific activity of the labelled cRNAs (pmol Cy3/μg cRNA) were measured by a NanoDrop ND-1000. Total of 1 μg of each labelled cRNA was fragmented by adding 5 μl 10 × Blocking Agent and 1 μl of 25 × Fragmentation Buffer, and then heating the mixture at 60 °C for 30 min. Finally 25 μl 2 × GE Hybridization buffer was added to dilute the labelled cRNA. Total of 50 μl of hybridization solution was dispensed into the gasket slide and assembled onto the lncRNA expression microarray slide. The slides were incubated for 17 h at 65 °C in an Agilent Hybridization Oven. The hybridized arrays were washed, fixed, and scanned using the Agilent DNA Microarray Scanner (part number G2505C). The microarray work was performed by Kang Chen Bio-tech, Shanghai, China.

Data analysis

Agilent Feature Extraction software (version 11.0.1.1) was used to analyse acquired array images. Then, the GeneSpring GX v12.1 software package (Agilent Technologies) was used to perform quantile normalization and subsequent data processing (Zahurak et al., 2007). Differentially expressed lncRNAs (DE-lncRNAs) and mRNAs (DE-mRNAs) between the two groups were identified through Fold Change filtering and P value with the R-based limma software package (Diboun et al., 2006). Finally, heatmap and volcano plots were generated for significantly differentially expressed genes using the “ggplot2 package” of R (Wickham, 2016). The threshold used to screen upregulated or downregulated lncRNAs and mRNAs was a | fold-change | ≥ 2.0 (P < 0.05). The microarray data discussed in this paper have been deposited in the NCBI Gene Expression Omnibus and are accessible with the GEO Series accession number GSE204935.

Bioinformatics analysis

To identify the potential role of differentially expressed genes in VSD, Gene Ontology (GO) annotation and the Kyoto Encyclopedia of Genes and Genomes (KEGG) were used to illuminate the functions of the DE-mRNAs (The Gene Ontology Consortium, 2017; Kanehisa & Goto, 2000). GO analysis provides the key functional classifications for DE-mRNAs, which include biological process (BP), cellular component (CC), and molecular function (MF). KEGG reveals the relationship networks and disease information for DE-mRNAs. Fisher’s exact test was used to classify the GO and KEGG categories. The threshold of significance was P value < 0.05 and FDR < 0.1.

Construction of the CNC network and ceRNA network

According to the above analysis, we selected core DE-mRNAs closely related to VSD for further study. Based on the expression level of all DE-lncRNAs and core DE-mRNAs, the Pearson correlation coefficient (PCC) was calculated by the “coXpress package” of R (Watson, 2006). The DE-lncRNAs and DE-mRNAs with PCC > 0.95 and P value < 0.05 were selected to construct the CNC network.

To further explore the potential role of lncRNAs, a ceRNA (lncRNA‒miRNA‒mRNA) network based on core DE-mRNAs was constructed. First, the miRcode database (v11) was used to predict miRNA targets on lncRNAs (Jeggari, Marks & Larsson, 2012), and the TargetScan database (v7.2) was used to predict miRNA targets on mRNAs (Agarwal et al., 2015). The thresholds for the algorithms were conservation scores >50 for miRcode, and target scores ≥70 for TargetScan. Then, the lncRNA‒miRNA‒mRNA interaction network was constructed using the overlapping miRNAs that harboured both lncRNA and mRNA binding targets. The Cytoscape (v2.8.2) program was used to construct and visualize the CNC network and ceRNA network (Shannon et al., 2003).

Quantitative real-time reverse transcription PCR (qRT‒PCR)

Total RNA was extracted from AF supernatant using TRIzol LS Reagent (Invitrogen, Carlsbad, CA, USA) and then cDNA was synthesized from 1 ug of total RNA using PrimeScript™ RT reagent Kit with gDNA Eraser (Takara, Dalian, China). cDNA was amplified and quantified using TB Green® Premix Ex Taq™ II (Takara, Dalian, China). The cycling parameters were as follows: the first step, 95 °C for 30 s; the second step, 95 °C for 5 s and 60 °C for 34 s, for 40 cycles; and the third step, for melting 95 °C for 15 s, 60 °C for 60 s and 95 °C for 15 s. Each sample was measured in triplicate and the sequences of the primers used are shown in Table S2. The obtained data of lncRNAs and mRNAs expression levels were calculated and normalized by the 2^−ΔΔCt method as follows: ΔCt_{(test/control)} = average value Ct_{(target gene)} − average value Ct_(GAPDH); ΔΔCt = ΔCt_{(test/control)} − average value ΔCt_(control).

Statistical analysis

Heatmaps and volcano plots based on the expression levels of genes were generated and visualized using R statistical software. Fisher’s exact test was employed to evaluate the significance of GO or KEGG enrichment. The expression levels of lncRNAs and mRNAs were calculated via an unpaired t-test between the VSD and control groups. All quantitative data are expressed as the mean ± SD and were analysed using SPSS 21.0 (SPSS Inc., Chicago, IL, USA). P < 0.05 was considered to indicate a statistically significant difference (Mushquash & O’Connor, 2006).

Differentially expressed lncRNAs and mRNAs

To evaluate the differential expression of AF-derived lncRNAs and mRNAs between the VSD and control groups, AF supernatant from five VSD foetuses and five control foetuses was obtained by a random method and subjected to microarray analysis. Box plots were used to compare the intensities and quality of microarray data (Figs. 1A and 1B). The result showed that the distributions of intensities in different samples were similar after normalization. A volcano plot was used to assess the variation in lncRNA and mRNA expression between the two groups (Figs. 1C and 1D). Finally, hierarchical clustering was used to show distinguishable lncRNA and mRNA expression profiles (Figs. 1E and 1F). With the cut-off criteria of a | fold change | ≥ 2 and P value < 0.05, 710 lncRNAs were significantly differentially expressed between the two groups, of which 240 were upregulated and 470 were downregulated in the VSD group. Similarly, 397 mRNAs were significantly differentially expressed between the two groups, of which 106 were upregulated and 291 were downregulated in the VSD group. The expression level of the up- and down-regulated DE-mRNAs and DE-lncRNAs was showed in Tables S3 and S4.

Characteristics of differentially expressed lncRNAs and mRNAs

To further analyze the DE-lncRNAs and DE-mRNAs, we summarized some general characteristics of these genes, such as type distribution, genomic location distribution, and length. As shown in Figs. 2 A and 2B, in both up- and down-regulated lncRNAs, intergenic lncRNAs are the largest types, and intron sense overlapping lncRNAs are the least type. The results of positioning showed that up- and down-regulated lncRNAs were mainly distributed on chromosomes 1 and 2, respectively, while both the up- and down-regulated mRNAs were mainly distributed on chromosomes 1 (Figs. 2C, 2D and S1). The length distribution showed that most DE-lncRNAs were distributed within 0.5–1 kb; however, the DE-mRNAs were mainly distributed at 2–4 kb in length. These results prompt that chromosome 1 and 2 may contain more genes related to the VSD process, and compared with other types of lncRNA, intergenic lncRNAs are more critical in the formation of cardiac defects.

GO and pathway analysis

To explore the mechanisms and pathways during the progression of VSD, the biological meaning and systematic features of the significant DE-mRNAs were first annotated by GO ontology. A total of 713 GO terms were significantly enriched. The top ten terms of the GO ontology results are shown in Fig. 3A. BP analysis indicated that the DE-mRNAs were enriched in heart development, cell adhesion, and regulation of cell proliferation. CC analysis indicated that the DE-mRNAs were enriched in cell junction, cell-substrate junction, and extracellular matrix. MF analysis indicated that the DE-mRNAs were enriched in kinase regulator activity and protein-containing complex binding.

KEGG pathway analysis was used to investigate the wire diagrams of molecular interactions, reactions, and relations related to the DE-mRNAs. The top ten enriched KEGG pathways are shown in Fig. 3B including the Sonic Hedgehog signaling pathway (Shh), the Phosphatidylinositol signaling system, apoptosis, and the regulation of the actin cytoskeleton. In addition, KEGG disease enrichment analysis showed that DE-mRNAs were involved in congenital malformations, cardiovascular diseases, and nervous system diseases (Fig. 3C). Of these, the most enriched disease was congenital malformations with 24 genes annotated with this term. These results show that DE-mRNAs are rich in biological activities and signal pathways related to heart development, including cell adhesion, cell proliferation and apoptosis, tissue development, and regulation of Shh signaling pathway.

Construction of the CNC network and ceRNA network

According to GO and KEGG analysis, four core mRNAs, including IDS, ID1, NR2F2, and GPC3, were selected to construct a CNC network because they are closely connected to cardiac biology. Among the coexpression analysis, four core mRNAs and 141 lncRNAs comprise the CNC network, which includes 149 pairs of coexpressing lncRNAs and mRNAs. As shown in Fig. 4A, a core mRNA could have a correlation with many lncRNAs, and most of these relationship pairs presented a positive correlation. To our knowledge, ceRNAs play important roles in the cis-regulation of lncRNAs. Therefore, we constructed a ceRNA network to reveal the regulatory mechanism between lncRNAs and mRNAs in the CNC network. In the ceRNA network, the miRcode and TargetScan databases were used to predict miRNA targets on lncRNAs and mRNAs, respectively. A total of 194 overlapping miRNAs were selected to construct the ceRNA network. To better show the specific ceRNA network, we constructed the network of up-regulated and down-regulated genes, respectively. Eventually, the down-regulated ceRNA network includes 11 lncRNAs, 129 miRNAs, and two mRNAs, while, the up-regulated ceRNA network includes four lncRNAs, 79 miRNAs, and two mRNAs. (Figs. 4B and 4C).

Validation of RNA expression by qPCR

According to the GO/KEGG analysis and CNC/ceRNA network, we examined the expression of four core mRNAs (IDS, ID1, NR2F2, and GPC3) and five lncRNAs (AL022323.5, LINC00598, GATA3-AS1, PWRN1, and LINC01551) to confirm the reliability of our microarray results. To increase the accuracy of the results, qPCR was performed in five VSD and five healthy samples used for microarray, and in six additional VSD and six additional healthy samples. The qPCR results showed that seven RNAs (IDS, NR2F2, GPC3, LINC00598, GATA3-AS1, PWRN1, and LINC01551) were consistent with the results of microarray analysis, while two RNAs (ID1, and AL022323.5) were not significantly different between the two groups (Table 2 and Fig. 5).