Ribosomal RNA‑depleted RNA sequencing reveals the pathogenesis of refractory Mycoplasma pneumoniae pneumonia in children.

Pneumonia caused by Mycoplasma pneumoniae (M. pneumoniae) is a major cause of community‑acquired pneumonia in children. In some cases, M. pneumoniae pneumonia (MPP) can develop into refractory MPP (RMPP), which shows no clinical or radiological response to macrolides, and can progress to severe and complicated pneumonia. However, the pathogenesis of RMPP remains poorly understood. The present study aimed to identify target genes that could be used as biomarkers for the clinical diagnosis of early‑stage RMPP through high‑throughput sequencing technology. The differences in long non‑coding (lnc)RNAs, mRNAs and circular (circ)RNAs were examined between whole‑blood samples from two patients with non‑refractory MPP (NRMPP), two patients with RMPP and three healthy children using ribosomal (r)RNA‑depleted RNA‑sequencing techniques and an integrated mRNA/circRNA analysis. A total of 17 lncRNAs (four upregulated and 13 downregulated), 18 mRNAs (six upregulated and 12 downregulated) and 24 circRNAs (12 upregulated and 12 downregulated) were the most significantly differentially expressed (P<0.05) between the NRMPP and RMPP groups. Upon functional analysis, the significantly differentially expressed genes encoded by the targeting mRNAs (prostaglandin‑endoperoxide synthase 2, IL‑8 and fos‑like antigen 1) were screened and identified to be enriched in the 'IL‑17 signaling pathway'. Furthermore, the key circRNAs in the NRMPP and RMPP comparative groups were primarily enriched in 'herpes simplex virus 1 infection', 'viral carcinogenesis' and 'RNA transport'. In the present study, a comprehensive analysis of the differences between the NRMPP and RMPP cases was performed based on rRNA‑depleted RNA‑sequencing techniques, and the selected genes and circRNAs may be closely associated with the complex pathogenesis of RMPP.

Introduction

Mycoplasma pneumoniae (M. pneumoniae) is one of the main pathogens that cause respiratory tract infections in humans (1,2). Outbreaks of M. pneumoniae pneumonia (MPP) occur in 3- to 7-year cycles worldwide, and epidemics in Korea occurred in 2015 and 2016 (3). M. pneumoniae causes respiratory tract infections in all age groups, accounting for up to 40% of community-acquired respiratory tract infections in children aged >5 years (1). Although MPP is usually a benign and self-limiting process, M. pneumoniae infection can develop into severe, life-threatening diseases, including refractory MPP (RMPP), acute respiratory distress syndrome and necrotizing pneumonitis (2–4). Pneumonia and extrapulmonary complications caused by MPP pose a serious threat to children's health (4). In previous years, MPP has been reported in 10–40% of cases of community-acquired pneumonia (CAP) and increases during M. pneumoniae epidemics (5,6). MPP was previously considered to be a self-limiting process, but RMPP has been reported (7), which shows no clinical or radiological response to macrolides, and may progress to severe pneumonia and cause extrapulmonary complications (6,7). The number of RMPP cases also tends to increase every year; the annual incidence of MPP in respiratory disease was <10% of cases of CAP in 2009, 20.5% in September 2010 and reached a record high of 62.5% in September 2011 (8). The usual treatment strategy used for MPP is not considered as suitable for RMPP (7,8); thus, an improved understanding of the pathogenesis underlying RMPP is urgently required in order to design a more effective treatment strategy.

Studies on RMPP have indicated that its pathogenesis is primarily associated with the macrolide-resistant mechanism of M. pneumoniae strains (9,10). Matsuoka et al (11) found that mutations in domains II and V of 23S RNA in M. pneumoniae strains result in a decrease in the affinity of antibiotics to bacterial ribosomes, which eventually leads to resistance to macrolides. In addition, abnormal immune responses caused by MPP (12) and mixed bacterial infections may cause progression of non-refractory MPP (NRMPP) to RMPP (13). However, due to its complexity, the key pathogenesis of RMPP remains unknown.

Long non-coding RNAs (lncRNAs) are >200 nucleotides in length and do not have the capacity to encode proteins (14). The mechanism of action underlying lncRNAs is not completely understood due to their poorly conserved nucleotide sequence. However, previous studies have reported that lncRNAs regulate gene expression at multiple levels through complex mechanisms (14,15). Circular RNAs (circRNAs) are a novel type of lncRNAs that form a covalently closed continuous loop and thus have no 5′-3′ polarity and no polyA tail (16). Accumulating evidence indicates that circRNAs regulate gene expression at both the transcriptional and post-transcriptional levels by serving as microRNA (miRNA) sponges (16,17). circRNAs serve a role in the development of certain diseases, including cancer, diabetes and Alzheimer's disease (16,17–19). Thus, circRNAs may serve as potential biomarkers or therapeutic targets.

The present study performed lncRNA and circRNA profiling via ribosomal (r)RNA-depleted RNA sequencing of whole-blood samples from two patients with NRMPP, two patients with RMPP and three healthy children (HC) to provide a comprehensive analysis of the differences between the HC and patients with MPP, including the main differences in RNA expression levels between patients with NRMPP and those with RMPP. The present study aimed to identify target genes and circRNAs that may serve as biomarkers for the clinical diagnosis of early-stage disease, as well as providing a theoretical basis for research into the pathogenesis of MPP.

Materials and methods

Whole-blood sample preparation

The 18 blood samples were collected from December 2018 to December 2019 in the Guangzhou Women and Children's Medical Center (Guangdong, China) from six patients with NRMPP, six patients with RMPP and six healthy children (HC; Table I). Among them, a cohort of seven children (two NRMPP, two RMPP and three HC) was used for high-throughput sequencing, and a cohort of 11 children (four NRMPP, four RMPP and three HC) was used for validation. A blood sample (2.0 ml) was collected from each patient. All patients were tested with the M. pneumoniae-IgM antibody in the serum and M. pneumoniae DNA by PCR in throat swabs on admission, and positive cases were defined as patients with MPP. All patients with MPP were treated with appropriate antibiotics (for example, Azithromycin). Cases with worsening cough, infiltrates on a chest radiograph and a fever that prolonged for >7 days were recorded as patients with RMPP. The remaining cases were defined as patients with NRMPP. All specimens were collected at an early stage of MPP on admission (≤7 days of onset was defined as an early stage of MPP) (20). The exclusion criteria were as follows: The presence of severe concomitant diseases (chronic pulmonary disease, cardiovascular disease, neoplasia, kidney or liver disease, immune function deficiency and immunodepression); and the presence of mixed infections with other microorganisms. The HC were used as healthy controls.

Table I.

Clinical characteristics of children diagnosed with MPP.

A, NRMPP
			Laboratory findings[b]								Radiology findings[c]
Case	Age	Sex	Fever duration[a], days	WBC, 1×109/l (5–12[d])	Hb, g/l (105–145[d])	Plt, 1×109/l (140–440[d])	HsCRP, mg/l (<5[d])	ALT, U/l (0–40[d])	LDH, U/l (159–322[d])	ALB, g/l (40–55[d])	Consolidation	Hydrothorax
1	1 year, 11 months	Male	4	5.7	115	399	46.6	27	295	39.5
2	2 years, 5 months	Female	7	6.4	120	264	1.7	16	390	40.1
3	11 years	Female	5	10.9	105	143	60.4	13	400	25.3	Yes
4	4 years, 1 month	Male	7	4.7	122	261	30.4	12	316	39.5	Yes
5	3 years, 2 months	Female	6	4.9	110	278	16.0	8	211	36.7
6	7 years, 2 months	Male	5	9.2	116	276	24.0	11	239	39.4	Yes
B, RMPP
			Laboratory findings[b]								Radiology findings[c]
Case	Age	Sex	Fever duration[a], days	WBC, 1×109/l (5–12[d])	Hb, g/l (105–145[d])	Plt, 1×109/l (140–440[d])	HsCRP, mg/l (<5[d])	ALT, U/l (0–40[d])	LDH, U/l (159–322[d])	ALB, g/l (40–55[d])	Consolidation	Hydrothorax
7	1 year, 9 months	Female	9	12.2	100	196	154.6	119	739	24.8	Yes
8	6 years	Female	13	5.8	102	233	155.8	23	1,470	25.3	Yes	Yes
9	3 years, 10 months	Female	7	10.1	106	495	31.9	48	332	28.1	Yes	Yes
10	9 years, 1 month	Female	12	9.7	113	254	35.0	719	1,442	32.8	Yes	Yes
11	5 years, 3 months	Male	10	7.4	119	319	25.3	25	651	28.2	Yes	Yes
12	4 years, 2 months	Male	9	5.9	101	322	73.1	45	941	26.4	Yes

Fever duration from onset to normothermia

Data extracted from the first test of the children on admission

Judged by chest radiograph or CT scan in whole course of the patient

Normal range for each laboratory finding. MPP, M. pneumoniae pneumonia; NRMPP, non-refractory MPP; RMPP, refractory MPP; WBC, white blood cell; Hb, hemoglobin; Plt, platelet; HsCRP, high-sensitivity C-reactive protein; ALT, alanine aminotransferase; LDH, lactate dehydrogenase; ALB, albumin.

The present study was approved by the Ethics Committee at Guangzhou Women and Children's Medical Center. The parents or legal guardians of the patients signed written informed consent forms and agreed to its content.

rRNA-depleted RNA sequencing

For rRNA-depleted RNA sequencing, two of the six patients with NRMPP, two of the six patients with RMPP and three of the six HC were screened out as typical cases and sent to Annoroad Gene Technology Co, Ltd. The total RNA of each sample was isolated using TRIzol® reagent (cat. no. MFCD00213058; Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. The purity, concentration and integrality of the RNA were determined using a NanoPhotometer® spectrophotometer (IMPLEN), Qubit®3.0 Fluorometer (Thermo Fisher Scientific, Inc.) and Agilent 2100 RNA Nano 6000 Assay kit (cat. no. 41105500; Agilent Technologies, Inc.), respectively. Subsequently, 3 µg RNA from each sample was loaded and the rRNA was removed using Ribo-Zero™ Gold kits (cat. no. MRZG126; Epicentre; Illumina, Inc.). NEBNext® Ultra™ Directional RNA Library Prep kit for Illumina® (cat. no. E7420S; New England BioLabs, Inc.) was used to generate the sequencing libraries, which were then sequenced through the Illumina HiSeq platform (Illumina, Inc.). The sequencing type was eukaryotic common transcriptome. The sequencing direction was P5 to P7, then P7 to P5. After removing the low-quality and polluted reads, clean reads were obtained and mapped to the reference genome sequence using Hierarchical Indexing For Spliced Alignment Of Transcripts 2 (version 2.05) (21). The detected reads were mapped to the known mRNA and lncRNA. HTSeq (version 0.6.0) was used to represent the expression level of each gene (22). The loading concentration of the final library was 2 nM. DESeq2 Rpackage (version 1.18.0) was used to perform differential expression analysis between the comparative groups (23), and the genes with P<0.05 were considered to be differentially expressed. The unmapped reads were identified as circRNA candidates using Find_circ (version 1.2) and CircRNA Identifier 2 software (version 2.0.1) (23,24). Functional analysis of the differentially expressed genes (DEGs) was performed using Gene Ontology (GO)seq (version 1.0) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology-Based Annotation System (version 2.0) (24–28).

Reverse transcription-quantitative PCR (RT-qPCR)

To verify the reliability of the sequencing results, RT-qPCR was performed to determine the expression levels of screened genes (Table SI). EditSeq software (version 7.10, DNASTAR, Inc.) was used to design the specific primers and the β-actin gene was selected as a standardization control. Briefly, total RNA from blood samples was reverse transcribed into cDNA using the QuantiTect Reverse Transcription kit (Qiagen GmbH) according to the manufacturer's instructions. Subsequently, qPCR was performed using a DNA Engine Chromo 4 real-time system (Bio-Rad Laboratories, Inc.) with a TaqMan™ Copy Number Assay kit (Thermo Fisher Scientific, Inc.). The sequences of the primers are listed in Table II. β-actin (forward, 5′-GAGGTATCCTGACCCTGAAGTA-3′ and reverse, 5′-CACACGCAGCTCATTGTAGA-3′) was used as an internal reference. Thermocycling conditions were as follows: 95°C for 10 min, followed 40 cycles of 95°C for 15 sec and 60°C for 60 sec). The expression levels were calculated using the 2−ΔΔCq method (29).

Table II.

Primers used for reverse transcription-quantitative PCR.

circRNA	Primers (5′→3′)
hsa_circ_ 0022808	F: ACTGAAGAGGATGCAGGAGC
	R: GAGGAATGTTCCCGGTCTCC
hsa_circ_ 0006793	F: AGTCCCCTGCTATCACTGGT
	R: CCAGCTTCGGTCACTGAACA
hsa_circ_ 0014390	F: TCTGTTGAAGATTTGAAGAACCCA
	R: CTGAGGGCTAGAGGACTGGT
hsa_circ_ 0014305	F: TTCTCCCTGGCGGAGGAATA
	R: GGATGGCTGGTTTGAAGCAC
hsa_circ_ 00216400	F: CTTCAGTACCAGAGCCCCAC
	R: ACCACCTCAACCGTTTCAACT

F, forward; R, reverse; circ, circular RNA.

Statistical analysis

Data were analyzed using GraphPad software (Prism 8.0; GraphPad Software, Inc.) and visualized using the ggplot2 package of R software (version 3.6.1) (30). All data are presented as the mean ± standard deviation of three independent repeats. An unpaired two-tailed Student's t-test was used to determine the significant differences in lncRNAs, mRNAs and circRNAs between the three groups. For clinical data, the non-parametric Mann-Whitney U test was used for two-group analysis of continuous variables, and the Kruskal-Wallis test followed by Dunn's post hoc test was used for three-group analysis of continuous variables. Fisher's exact test was used for the analysis of categorical variables. P-values are two-sided and were adjusted using the Bonferroni method. P<0.05 was considered to indicate a statistically significant difference.

Results

Characteristics of patients with MPP

Patients with MPP displayed typical MPP clinical symptoms and were diagnosed with pneumonia. The M. pneumoniae IgM antibody in the serum and M. pneumoniae DNA detected via PCR from throat swabs showed positive results. The clinical characteristics of the cases are listed in Tables SI and SII. Abnormal findings on chest radiographs were observed in all patients with MPP. The chest scans primarily revealed diffuse infiltration of both lungs in NRMPP cases. However, RMPP cases exhibited unequivocal focal or segmental consolidation (Fig. 1A) with pleural effusion (Fig. 1B).

Figure 1.

Imaging features of lung injury associated with MPP in children. (A) High-resolution computerized tomography scan of the chest revealing areas of airspace consolidation in the right middle lobe in a 7-year-old child with non-refractory MPP. (B) Areas of airspace consolidation and pleural effusion in the left lobe in a 4-year-old child with refractory MPP. MPP, M. pneumoniae pneumonia.

Overview of rRNA-depleted RNA sequencing analysis profiles

Upon rRNA-depleted RNA sequencing, a total of 670,528,498 raw reads (102,839,708, 98,869,356 and 93,692,076 for the HC group; 96,341,690 and 92,905,200 for the NRMPP group; and 96,437,832 and 89,442,636 for the RMPP group) were generated. Through filtering the raw data, 614,648,054 clean reads were obtained and mapped to the reference genome (Fig. 2A). The error rate of the filtered data was qualified (Fig. 2B).

Figure 2.

Quality control of rRNA-depleted RNA sequencing. (A) Proportion of various reads before filtration. (B) Q30 ratio (Error rate <0.1%) of sequencing data in each sample. rRNA, ribosomal RNA; NS, sequence-dependent ambiguous bases; MPP, M. pneumoniae pneumonia; NRMPP, non-refractory MPP; RMPP, refractory MPP; HC, healthy children.

A total of 611 lncRNAs (416 upregulated and 195 downregulated) and 692 mRNAs (598 upregulated and 94 downregulated) were significantly differentially expressed (P<0.05) between the NRMPP and HC groups (Fig. 3). A total of 937 lncRNAs (433 upregulated and 504 downregulated) and 1,027 mRNAs (593 upregulated and 434 downregulated) were significantly differentially expressed (P<0.05) between the RMPP and HC groups. A total of 17 lncRNAs (4 upregulated and 13 downregulated) and 18 mRNAs (6 upregulated and 12 downregulated) were significantly differentially expressed (P<0.05) between the RMPP and NRMPP groups (Table III). The significantly differentially expressed lncRNAs between the RMPP and NRMPP groups included ENSG00000249790, ENSG00000261026, MSTRG.215206, MSTRG.233743, MSTRG.238033, MSTRG.238419, MSTRG.268000, MSTRG.275241 (Table III).

Figure 3.

Differential analysis of the DEGs and lncRNAs. (A) Volcano map of DEGs. (B) Cluster map of differentially expressed lncRNAs. (C) Cluster map of differentially expressed mRNAs. Blue represents downregulated, black represents no change and yellow represents upregulated. DEGs, differentially expressed genes; lncRNA, long non-coding RNA; MPP, M. pneumoniae pneumonia; NRMPP, non-refractory MPP; RMPP, refractory MPP; HC, healthy children.

Table III.

Significantly differentially expressed long non-coding RNAs between RMPP and NRMPP groups.

	Relative gene expression levels
Gene	RMPP	NRMPP	P-value	Up/downregulated	Position
ENSG00000249790	635.1540951	3447.840859	3.09×10−3	Down	chr12:8788257-8795789:+
ENSG00000261026	472.7296308	4.781416405	2.65×10−2	Up	chr8:22679013-22684009:−
MSTRG.215206	2391.068595	70.25707389	7.87×10−3	Up	chr15:48211935-48296341:−
MSTRG.233743	63.78915178	0	1.70×10−2	Up	chr16:15028952-15030258:−
MSTRG.238033	0	27.68733409	2.65×10−2	Down	chr16:34582662-34587722:+
MSTRG.238419	0.404276272	87.94800239	2.34×10−6	Down	chr16:46388508-46406972:+
MSTRG.268000	5.600197692	101.7282093	3.32×10−3	Down	chr17:74562930-74563656:−
MSTRG.275241	0	25.47228864	3.27×10−2	Down	chr18:22348684-22349256:−
MSTRG.299858	0.404276272	34.38883192	2.65×10−2	Down	chr19:37282731-37293921:−
MSTRG.32896	0.404276272	54.56033482	7.16×10−4	Down	chr1:143501630-143505359:+
MSTRG.384208	0	83.06200226	1.08×10−6	Down	chr21:8219668-8220285:−
MSTRG.384231	0	78.99033548	1.45×10−6	Down	chr21:8446933-8448138:−
MSTRG.386206	0	28.50166744	2.23×10−2	Down	chr21:22451534-22451867:+
MSTRG.488193	0	25.24433402	4.13×10−2	Down	chr4:146285918-146287333:+
MSTRG.631597	321.397898	16.70145272	2.65×10−2	Up	chr8:32733246-32739937:+
MSTRG.678280	0.490427811	39.44287588	9.98×10−3	Down	chr9:101736669-101848060:+

MPP, M. pneumoniae pneumonia; NRMPP, non-refractory MPP; RMPP, refractory MPP.

Bioinformatics analysis of sequencing profiles

The differentially expressed mRNAs between the RMPP and NRMPP groups (Table IV) were ENSG00000073756, ENSG00000111788, ENSG00000122877, ENSG00000123838, ENSG00000130656, ENSG00000165949; these were identified by functional analysis using both the GO and KEGG databases. In the GO analysis, the significantly differentially expressed mRNAs were primarily enriched in ‘complement activation, classical pathway’, ‘leukocyte migration’ and ‘chemotaxis’ (Fig. 4A). In the KEGG pathway analysis, the significantly differentially expressed mRNAs were primarily enriched in the ‘IL-17 signaling pathway’ (Fig. 4B). RT-qPCR was used to verify candidate genes that may be involved in pathogenesis of RMPP, such as prostaglandin-endoperoxide synthase 2 (PTGS2), chemokine (C-X-C motif) ligand 8 (CXCL8) and Fos-like antigen 1 (FOSL1; Fig. 4C) and the primers were designed by EditSeq software.

Table IV.

Significantly differentially expressed mRNAs between RMPP and NRMPP groups.

	Relative gene expression levels
Gene	RMPP	NRMPP	P-value	Up/downregulated	Gene name	Description
ENSG00000073756	9780.85878	553.0410947	1.67×10−2	Up	PTGS2	Prostaglandin-endoperoxide synthase 2
ENSG00000111788	9.669441835	103.3909011	3.10×10−2	Down	DDX11L8	Putative ATP-dependent RNA helicase DDX11-like protein 8
ENSG00000122877	536.8492789	21.50520559	1.37×10−3	Up	EGR2	Early growth response protein 2
ENSG00000123838	240.3692975	4.57224895	2.65×10−2	Up	C4BPA	Complement component 4 binding protein α
ENSG00000130656	0	43.9105411	9.88×10−4	Down	HBZ	Hemoglobin subunit ζ
ENSG00000165949	23.52076078	359.2170269	1.08×10−6	Down	IFI27	Interferon α-inducible protein 27 mitochondrial
ENSG00000169429	5129.9448	284.2008731	9.35×10−3	Up	CXCL8	Interleukin 8
ENSG00000175592	137.6919176	5.491166033	3.22×10−2	Up	FOSL1	Fos-like antigen 1
ENSG00000185736	0.490427811	190.8606578	1.14×10−9	Down	ADARB2	Double stranded RNA-specific editase B
ENSG00000188056	0.490427811	75.33804428	1.06×10−5	Down	TREML4	Trem-like transcript 4 protein
ENSG00000211653	82.7834293	501.0131059	5.00×10−3	Down	IGLV1-40	Immunoglobulin lambda variable 1–40
ENSG00000211897	570.423636	4031.523518	8.08×10−5	Down	IGHG3	Immunoglobulin heavy constant γ3
ENSG00000236770	0.404276272	43.86941747	4.53×10−3	Down	CD300C	CMRF35-like molecule 6
ENSG00000237973	1159.431015	6.701497833	4.07×10−5	Up	COX1	Cytochrome c oxidase subunit 1
ENSG00000270550	3.664967871	81.06781283	1.81×10−3	Down	IGHV3-30	Immunoglobulin heavy variable 3–30
ENSG00000282639	0	59.68597813	1.06×10−5	Down	IGHV3-64D	Immunoglobulin heavy variable 3–64D
ENSG00000282651	0	108.602341	6.73×10−9	Down	IGHV5-10-1	Immunoglobulin heavy variable 5-10-1
ENSG00000284690	43.90687745	553.5974251	3.09×10−3	Down	CD300H	Protein CD300H

MPP, M. pneumoniae pneumonia; NRMPP, non-refractory MPP; RMPP, refractory MPP.

Figure 4.

Functional analysis of the DEGs between the RMPP and NRMPP groups. (A) GO analysis of the DEGs. (B) DEGs enriched in the ‘IL-17 signaling pathway’. (C) Validation of DEGs via reverse transcription-quantitative PCR. For the presented data, four NRMPP, four RMPP and three HC were used. Comparisons between two groups were analyzed using Student's t-test. *P<0.05 and **P<0.01 vs. NRMPP group. DEGs, differentially expressed genes; GO, Gene Ontology; MPP, M. pneumoniae pneumonia; NRMPP, non-refractory MPP; RMPP, refractory MPP; HC, healthy children; MPPS, severe MPP; MPPM, mild MPP; PTGS2, prostaglandin-endoperoxide synthase 2; CXCL8, interleukin 8; FOSL1, fos-like antigen 1; HBZ, hemoglobin subunit ζ; IGHV, immunoglobulin heavy variable.

In the circRNA/miRNA co-expression analysis, a total of 1,370 circRNAs (505 upregulated and 865 downregulated) were significantly differentially expressed (P<0.05) between the HC and MPP groups (Fig. 5A and B). The functions of circRNAs were associated with the known function of the host linear transcripts and annotated by the GO and KEGG databases. In the GO analysis of the host linear transcripts, the differentially expressed terms were classified into three categories. Under biological processes, the GO terms were primarily enriched in ‘regulation of mRNA metabolic process’, ‘nucleobase-containing compound transport’ and ‘RNA localization’ (Fig. 5D). Under the category of cellular component, the GO terms were primarily enriched in ‘nuclear speck’, ‘cell-substrate junction’ and ‘focal adhesion’ (Fig. 5E). Under the category of molecular function, the GO terms were primarily enriched in ‘ubiquitin-like protein transferase activity’, ‘ribonucleoprotein complex binding’ and ‘enhancer binding’ (Fig. 5F). In the KEGG pathway analysis, the significant DEGs were primarily enriched in ‘Herpes simplex virus 1 infection’, ‘viral carcinogenesis’ and ‘RNA transport’ (Fig. 5C). The top 11 significantly differentially expressed circRNAs between the HC and MPP groups are listed in Table V.

Figure 5.

Functional analysis of the circRNAs between the HC and MPP groups. (A) Volcano plot and (B) cluster map of significantly differentially expressed circRNAs between the HC and MPP groups. Red represents upregulated, blue represents downregulated and yellow represents no change. (C) The most enriched pathways between the HC and MPP groups. The most enriched Gene Ontology terms classified into three main categories: (D) Biological processes, (E) cellular component and (F) molecular function. circRNA, circular RNA; HC, healthy children; MPP, M. pneumoniae pneumonia; MPPS, severe MPP; MPPM, mild MPP.

Table V.

Significantly differentially expressed circRNAs between HC and MPP groups.

Gene	Log2FoldChange	P-value	Best transcript	Gene symbol	Upregulated/downregulated
hsa_circ_0002171	10.42052522	0.029534	NM_198567	C5orf25	Up
hsa_circ_0001535	7.929259862	0.0002309	NM_001101801	FAM13B	Up
hsa_circ_0007261	2.669813966	0.0254796	NM_031921	ATAD3B	Up
hsa_circ_0018432	2.28731687	0.0254796	NM_194298	SLC16A9	Up
hsa_circ_0014535	−1.810846752	0.0242991	NM_002004	FDPS	Down
hsa_circ_0002316	−1.891464967	0.0254796	NM_145243	OMA1	Down
hsa_circ_0004327	−1.973037708	0.0429052	NM_007124	UTRN	Down
hsa_circ_0019868	−2.107107635	0.0026471	NM_014631	SH3PXD2A	Down
hsa_circ_0001726	−2.465790354	0.002007	NM_032164	ZNF394	Down
hsa_circ_0026176	−2.760620461	0.0254796	NM_003217	TMBIM6	Down
hsa_circ_0021682	−4.772727352	0.0242991	NM_024662	NAT10	Down

MPP, M. pneumoniae pneumonia; HC, healthy children; circRNA, circular RNA.

A total of 156 circRNAs (85 upregulated and 71 downregulated) were significantly differentially expressed (P<0.05) between the NRMPP and RMPP groups (Fig. 6A-C). A total of 24 circRNAs were identified as the most significantly differentially expressed circRNAs between the NRMPP and RMPP groups. In the GO analysis, GO terms were primarily enriched in ‘positive regulation of myeloid cell differentiation’ and ‘positive regulation of hemopoiesis’ (Fig. 6D). The screened circRNAs (Table VI) were primarily enriched in ‘colorectal cancer’, ‘hepatitis B’ and ‘apoptosis’ (Fig. 6E). The top five upregulated circRNAs were selected for further validation by performing RT-qPCR (Fig. 6F), indicating that these circRNAs may serve as potential biomarkers for RMPP.

Figure 6.

Functional analysis of the circRNAs between the NRMPP and RMPP groups. (A) Volcano plot, (B) Lateral view of the volcano plot and (C) cluster map of significantly differentially expressed circRNAs between the NRMPP and RMPP groups. (D) Most enriched Gene Ontology terms. (E) Most enriched Kyoto Encyclopedia of Genes and Genomes pathways between the NRMPP and RMPP groups. (F) Validation of the expression of the selected circRNAs via reverse transcription-quantitative PCR. For the presented data, four NRMPP, four RMPP and three HC were used. Comparisons between two groups were analyzed using Student's t-test. *P<0.05 and **P<0.01 vs. NRMPP group. circRNA, circular RNA; MPP, M. pneumoniae pneumonia; NRMPP, non-refractory MPP; RMPP, refractory MPP; HC, healthy children; MPPS, severe MPP; MPPM, mild MPP.

Table VI.

Significantly differentially expressed circRNAs between RMPP and NRMPP groups.

Gene	Log2FoldChange	P-value	Best transcript	Gene symbol	Upregulated/downregulated
hsa_circ_0022808	23.10152809	1.38×10−6	NM_002689	POLA2	Up
hsa_circ_0006793	7.812895001	4.88×10−2	NM_001242767	MTHFD1L	Up
hsa_circ_0014390	6.98745233	3.28×10−3	NM_014847	UBAP2L	Up
hsa_circ_0014305	6.373176867	9.65×10−3	NM_023015	INTS3	Up
hsa_circ_0021640	6.166346606	1.54×10−2	NM_005898	CAPRIN1	Up
hsa_circ_0017689	6.1204205	1.71×10−2	NM_014688	USP6NL	Up
hsa_circ_0025915	6.067184847	1.80×10−2	NM_025003	ADAMTS20	Up
hsa_circ_0003781	5.86144273	12.85×10−2	NM_001242865	PRMT2	Up
hsa_circ_0024849	5.834493906	2.66×10−2	NM_001143835	NFRKB	Up
hsa_circ_0015519	5.737917011	3.35×10−2	NM_004736	XPR1	Up
hsa_circ_0005901	5.506362304	4.75×10−2	TCONS_00024174	TCONS_00024174	Up
hsa_circ_0023925	5.506362304	4.75×10−2	NM_007166	PICALM	Up
hsa_circ_0025400	−4.680048526	3.23×10−2	NM_004426	PHC1	Down
hsa_circ_0010030	−5.606173248	4.84×10−2	NM_201628	KAZN	Down
hsa_circ_0011249	−5.645431582	4.67×10−2	NM_001020658	PUM1	Down
hsa_circ_0004477	−5.753705882	3.74×10−2	NM_182751	MCM10	Down
hsa_circ_0025209	−5.791153915	3.88×10−2	NM_001033714	NOP2	Down
hsa_circ_0018110	−5.819081289	3.37×10−2	NM_004521	KIF5B	Down
hsa_circ_0010131	−5.889114841	3.05×10−2	NM_004431	EPHA2	Down
hsa_circ_0010148	−6.072531801	2.51×10−2	NM_198546	SPATA21	Down
hsa_circ_0026524	−6.476462428	1.13×10−2	NM_001417	EIF4B	Down
hsa_circ_0001890	−7.728346085	6.81×10−3	NM_001006617	MAPKAP1	Down
hsa_circ_0014349	−9.677465499	4.32×10−2	NM_002870	RAB13	Down
hsa_circ_0022807	−15.16379666	4.45×10−13	NM_002689	POLA2	Down

MPP, M. pneumoniae pneumonia; NRMPP, non-refractory MPP; RMPP, refractory MPP; circRNA, circular RNA.

Discussion

Although the macrolide-resistant mechanisms of M. pneumoniae strains and excessive immunological inflammation are the most commonly proposed mechanisms underlying RMPP (3,4), the pathogenesis of RMPP remains to be elucidated and there is still a lack of accurate assessment tools and biomarkers for RMPP. At present, the common methods for estimating the severity of RMPP are based on clinical characteristics, pulmonary imaging severity and therapeutic effect, which are unable to ensure an effective identification of early-stage RMPP (5,6). Therefore, it is necessary to identify novel tools and biomarkers for the early diagnosis of RMPP. The present study was designed to identify target genes implicated in the pathogenesis of RMPP to enable early diagnosis by comparing the differences between the cases with NRMPP and those with RMPP. To the best of our knowledge, the present study was the first to assess the differences in lncRNAs and circRNAs between NRMPP and RMPP.

circRNAs serve important roles in regulating gene expression by sequestering miRNAs as a sponge at the transcriptional or post-transcriptional levels (16). Thus, circRNAs can regulate a number of processes associated with numerous diseases, such as cancer (16,17). M. pneumoniae possesses a tip-like organelle that permits a highly oriented extracellular parasitism of the respiratory epithelium to avoid clearance by mucosal cilia and phagocytosis, and its adhesion ability is positively correlated with pathogenicity (31). Upregulation of PTGS2 promotes inflammation, which may indicate that more severe inflammation was observed in the RMPP group in the present study (32). IL-8 is a chemotactic factor that attracts neutrophils, basophils and T cells, and it is also involved in neutrophil activation (33). FOSL1 encodes the regulator protein and is involved in cell proliferation, differentiation and transformation (34). The expression of PTGS2, IL-8 and FOSL1 was significantly higher in the RMPP group compared with the NRMPP group in the present study, indicating that the upregulation of these proteins may participate in the pathogenicity of RMPP. In addition, RMPP cases exhibit a high activation level of the IL-17 signaling pathway, which may cause an autoimmune response and disease aggravation (35). Immunoglobulin heavy variable (IGHV)3-30, IGHV3-64D and IGHV5-10-1 belong to the V region of the variable domain of immunoglobulin heavy chains that participate in antigen recognition (36). In the present study, the disappearance of IGHV3-64D and IGHV5-10-1 genes and low expression of the IGHV3-30 gene in the RMPP group may be an important mechanism that leads to RMPP cases due to antigen recognition problems. However, further experiments are required to confirm these hypotheses.

In the present circRNA/mRNA analysis, circRNA function was found to be associated with the known function of the host linear transcripts. Based on the circRNA/miRNA/mRNA analysis conducted in the present study, several differentially expressed mRNAs were identified to be associated with the differentially expressed circRNAs. A total of 11 circRNAs were identified as the most significantly differentially expressed circRNAs between the HC and MPP groups. Among those, hsa_circ_0019868 [SH3 and PX domain-containing protein 2A (SH3PXD2A) gene] may be associated with MPP pathogenesis. SH3PXD2A is an adapter protein involved in the invasiveness of cancer cells that mediates the neurotoxic effect of the amyloid-β peptide (37). The higher expression in MPP groups may enhance the invasiveness of M. pneumoniae. A total of 24 circRNAs were identified as the most significantly differentially expressed circRNAs between the NRMPP and RMPP groups. Among their target genes, hsa_circ_0001890 [target of rapamycin complex 2 subunit MAPKAP1 (MAPKAP1) gene], hsa_circ_0026524 (eukaryotic translation initiation factor 4B gene), hsa_circ_0021640 [caprin-1 (CAPRIN1) gene], hsa_circ_0003781 [protein arginine N-methyltransferase 2 (PRMT2) gene], hsa_circ_0010131 [ephrin type-A receptor 2 (EPHA2) gene], hsa_circ_0025209 [NOP2 nucleolar protein (NOP2) gene] and hsa_circ_0023925 [phosphatidylinositol-binding clathrin assembly protein (PICALM) gene] may be associated with the pathogenesis of RMPP. The identified genes display the following characteristics: MAPKAP1 is involved in ciliogenesis, regulates cell proliferation and survival (38), and may serve as a novel anti-infection and antifibrogenesis genomic locus in chronic schistosomiasis japonica (39); NOP2 is involved in a ribosomal large subunit assembly and regulates cell proliferation (40); PRMT2 is involved in the regulation of proliferation and promotes apoptosis (41); EPHA2 regulates migration, integrin-mediated adhesion, proliferation and differentiation of cells (42); CAPRIN1 may regulate cell proliferation and migration in multiple cell types (43); and PICALM serves an important role in several processes, such as internalization of cell receptors, synaptic transmission and removal of apoptotic cells (44). Therefore, these genes might be involved in the pathogenesis of RMPP, but further investigations are required. In addition, the selected circRNAs (hsa_circ_0022808, hsa_circ_0006793, hsa_circ_0014390, hsa_circ_0014305 and hsa_circ_00216400) may represent valuable markers for the diagnosis of patients with early-stage RMPP and NRMPP. However, there were limited samples used in the present study; therefore, future studies should use larger sample sizes.

To conclude, the present study provided a comprehensive analysis of the expression levels of different lncRNAs, mRNAs and circRNAs between NRMPP and RMPP cases using rRNA-depleted RNA-sequencing techniques. The selected genes or circRNAs may aid with identifying the complex pathogenesis of RMPP and determining the diagnostic and therapeutic value of circRNAs in RMPP.