Differentially Expressed Circulating Long-Noncoding RNAS in Premature Infants with Respiratory Distress Syndrome

This prospective study enrolled 15 premature infants who were admitted to Jiangyin People’s Hospital of Nantong University between April 2019 and October 2019.

Inclusion and Exclusion Criteria: Infants were eligible for enrollment in the study if they were (1) with a gestational age less than 36 weeks; (2) admitted within 4 hours after birth; (3) appropriate for gestational age. Patients were excluded for any of the following reasons: (1) severe cyanotic congenital heart diseases; (2) congenital chromosomal diseases or severe congenital malformations; (3) severe asphyxia at delivery (5 min Apgar score <5); (4) early symptoms of sepsis [10].

Severity grading: In the present study, the severity of RDS was determined clinically using a combination of PS treatment coupled with a degree of aeration of the lungs on chest X-ray [11]. The degree of aeration of the lungs on chest X-ray was graded as follows: (1) slightly reduced radiolucency with still sharp cardiac and diaphragmatic margins; (2) markedly reduced radiolucency with retained cardiac and diaphragmatic margins; (3) severely reduced radiolucency with air bronchogram and blurred cardiac and diaphragmatic margins; and (4) almost completely white lung fields with or without air bronchogram and barely visible cardiac and diaphragmatic margin [12].

Grouping: Of the 15 included infants, 5 were neonates without RDS and 10 were newborns diagnosed with RDS (presenting as cyanosis, groan, intercostal retractions, polypnea, and nasal flaring combined with changed aeration of the lungs on chest X-ray [11]). Babies who were worsening when FiO₂ >0.30 or positive end-expiratory pressure (PEEP) > 6 cm H₂O were given PS replacement [11]. Infants were given PS only once with lung X-ray grade 1or 2 were further defined as mild RDS, while infants who needed PS re-dosing with X-ray grade 3 or 4 were defined as severe RDS. Accordingly, based on the patients’ clinical features and chest X-ray results, 10 RDS infants were further divided into mild RDS group (n=5) and severe RDS group (n=5).

Data collecting: Data provided by all patients were collected in detail using a standard data collection form, including age, gender, gestational age, weight, 5 min Apgar score, maternal gestational diabetes mellitus, antenatal glucocorticoid use, etc. The collection was completed by two individuals independently and verified by a third person. Patient information has been processed anonymously before statistical analysis.

Peripheral blood samples (2ml for each person) were collected from all infants between 1 and 6 hours after birth. Among them, it should be noted that, for RDS patients, samples were drawn before PS replacement. All blood samples were frozen in the −80°C refrigerator following a specific process which includes centrifugation at 3,000 × g for 10 min at 4°C and then separation of clear upper liquid into an RNase-free tube. Total RNA was then extracted from the blood samples using the TRIzol reagent according to the manufacturer’s instructions and a previous study [13]. After quality control of RNA, the RNA library of each sample was prepared using NEB Next Ultra RNA Library Prep Kit for the Illumina platform (BioLabs Inc., USA). The RNA sequencing analysis was performed by Genminix Informatics Co., Ltd (Shanghai, China) with the GeneChip® Human Transcriptome Array 2.0 (Affymetrix Inc., US) served as a gene expression profiling tool.

The expression profile of the lncRNAs were analyzed by Deseq package (Affymetrix Inc., US). Samples were hybridized on the Human Clariom D (Thermo Fisher Scientific) gene chip. Background-adjustment, normalization, and log-transformation of signals intensity were performed with the Signal Space Transformation-Robust Multi-Array Average algorithm (RMA). Raw data were analyzed by the transcriptome analysis console (TAC) 4.0 software (Applied Biosystems, Foster City, CA, USA) awaiting further analysis [14]. The differentially expressed lncRNAs and mRNAs were screened according to the criteria of gene differential expression with |log2-fold change| (FC) more than 2 times and adjusted P<0.05. The differentially expressed lncRNAs were afterward clustered by a heatmap package while the hierarchical clustering diagram was drawn to show the results.

For lncRNA expression analysis, total RNA was transcripted to cDNA using a Reverse Transcription Kit (PrimeScript RT Master Mix, Takara Bio Inc., Otsu, Japan), then real-time quantitative PCR (qRT-PCR) validation was performed using the SYBR method (SYBR® Premix Ex Taq™, Takara Bio Inc., Otsu, Japan) according to the product instructions. An aliquot of 1 µg total RNA was added to each reaction mixture. qRT-PCR was performed on an ABI 7500 thermal cycler (Applied Biosystems; Thermo Fisher Scientific, Inc. US) with SYBR Green (Roche Diagnostics Co., Ltd. GER). The thermocycling conditions were as follows: 95°C for 5 min, followed by 40 cycles of 95°C for 20 sec and 55°C for 20 sec. At the end of each run, a melting curve analysis was performed at 72°C to monitor primer dimers and formation of non-specific products. For data analysis, the comparative Ct method (2^ΔΔCt) was used. Results were expressed as fold changes of gene expression adjusted to housekeeping gene GAPDH [15]. All primers used in the present study were shown in table 1.

Multivariate statistical analysis was used to calculate the Pearson correlation coefficient between differentially expressed lncRNAs and mRNAs. The greater the correlation coefficient, the greater possibility that there was a regulatory relationship between certain lncRNAs and mRNAs. The co-expression network was constructed with the Pearson correlation coefficient r> 0.99 and P< 0.05 as the filtering standard in this study.

GO and KEGG pathway analysis were applied to predict functions of the differentially expressed genes. The GO project offers a controlled vocabulary to label gene and gene product attributes in any organism (geneontology.org). GO results were mainly classified into three subgroups namely biological process, cellular component, and molecular function. GO analysis provides an interpretation of the relevance of genes differentially expressed between the groups. Fisher’s exact test and the χ² test were performed to calculate the P-value and false discovery rate of each GO term function. KEGG (kegg.jp/kegg/pathway.html) pathway analysis is a functional analysis tool, mapping a set of genes that may be associated with a certain lncRNA to potential pathways. The enrichment probability of a differentially expressed gene set in a term entry was represented by an enrichment score (EC), with a higher EC indicating a higher significance of the entry. The EC was calculated as the negative base 10 log of the P value. The input used in the bioinformatics analysis was the differential mRNA genes co-expressed with lncRNA that were screened in the lncRNA expression profile.

For clinical results (clinical characteristics), data were analyzed using SPSS 17.0 software. Quantitative data are expressed as mean ± standard deviation (SD). One-way variance analysis was applied to detect differences among the three groups. In terms of qualitative data, the Pearson Chi-square test was performed. Significant differences were considered as P < 0.05.

This present study was comprised of 15 premature infants in total, 5 cases without RDS for control (named as Group 1, G1), 5 cases with mild RDS (named as Group 2, G2) and 5 with severe RDS (named as Group 3, G3). The recruitment procedures are shown in Figure 1, and clinical characteristics of the infants are summarized in table 2. There were no significant differences in gestational age, birth weight, 5 min Apgar score, gender, mode of delivery, twin pregnancy, mother of gestational diabetes, and prenatal glucocorticoid (P>0.05).

Figure 1.

Affymetrix Human GeneChip was utilized to determine the expression spectrum of lncRNAs. As a result, the G1 vs. G2 comparison showed a total of 10112 differentially expressed lncRNAs, while the G2 vs. G3 comparison showed a total of 4663 differentially expressed lncRNAs. Of them, 135 lncRNAs were indicated to be differentially expressed among all three groups (G1, G2, and G3) after fold-change filtering (adjusted P value<0.05 and |log2-fold change|>2). The information of the top 10 upregulated and downregulated lncRNAs are listed in table 3. A hierarchical clustering map is presented to distinguish lncRNA expression profiles among the three groups. (Figure 2)

Figure 2.

Furthermore, differentially expressed mRNAs were compared for target prediction. Of them, the comparison between G1 and G2 showed a total of 2520 differentially expressed mRNAs, while the comparison between G2 and G3 showed a total of 530 mRNAs. The comparison of the three groups showed a total of 616 differentially expressed mRNAs. The lncRNA-mRNA co-expression network was constructed and showed a complex interaction between lncRNAs and mRNAs. Our analysis finally identified a total of 278 mRNAs closely related to 108 upregulated lncRNAs and 27 downregulated lncRNAs. These mRNAs with FC> 2 are shown in table 4.

GO and KEGG analysis were further performed to annotate the biological functions of differentially expressed mRNAs. The GO analysis indicated that the mRNAs co-expressed with 108 upregulated lncRNAs were associated with 247 GO terms. The top 25 enriched terms are shown in Figures 3A and 3B.

Figure 3.

Additionally, a KEGG pathway analysis was performed to investigate the possible roles of the lncRNA-associated mRNA genes. The most significant pathways enriched in the set of upregulated protein-coding genes included PI3 kinase/Akt (PI3K-Akt), RAS, and mitogen-activated protein kinase (MAPK) signal pathways, while the most significant KEGG pathways of the downregulated protein-coding genes were mainly related to metabolic pathways, etc. The bubble diagrams of the top KEGG pathways of mRNAs co-expressed with upregulated and downregulated lncRNAs are shown in Figure 3C and 3D.

Following the screening, four lncRNAs including ENST00000470527.1, ENST00000504497.1, ENST00000417781.5, and ENST00000440408.5 were further confirmed by qRT-PCR. Compared with G2 and G1, the expression levels of these four lncRNAs were increased in G3, which is consistent with the results of RNA sequencing. The relative expression levels are shown in Figure 4.

Figure 4.

In-depth bioinformatics analysis of lncRNAs showed all the four lncRNAs were involved in the MAPK signaling pathway by down-regulating gene GRB2 and MECOM. Moreover, ENST00000417781.5 and ENST00000440408.5 may regulate the MAPK signaling pathway by down-regulating gene IGF2, while ENST00000440408.5 and ENST00000504497.1 may target the MAPK signaling pathway by up-regulating gene EFNA1 and PLA2G4F, respectively. This study also showed that the above four lncRNAs participate in PI3KAkt and RAS signaling pathway by down-regulating GRB2, while ENST00000417781.5 and ENST00000440408.5 could regulate PI3K-Akt pathway by down-regulating IGF2 and up-regulating ITGB8 and TCL1B.

In addition, three of lncRNAs including ENST00000417781.5, ENST00000470527.1, and ENST00000504497.1 could target the RAS signaling pathway by up-regulating RASAL1, while ENST000004177, as well as ENST00000440408.5 could regulate RAS pathway by down-regulating IGF2. ENST00000440408.5 and ENST00000504497.1 may be involved in the RAS signaling pathway by up-regulating EFNA1 and PLA2G4F, respectively. LncRNAs including ENST00000417781.5, ENST00000470527.1, and ENST00000504497.1 could participate in the TGF-β signaling pathway by promoting gene expression of AMH and inhibiting TGIF2. In addition, ENST00000470527.1 and ENST00000504497.1 could be involved in the TGF-β pathway by up-regulating GDF7 and down-regulating GDF6, while ENST00000440408.5 may down-regulate FST to be involved in TGF-β pathway. The pathway regulatory network of four validated lncRNAs are shown in Figure 5.

Figure 5.