We carried out a single-center retrospective cohort study of children ≤16 years with pneumonia attributed to Mycoplasma pneumoniae, admitted to the Affiliated Hospital of Shandong University of Traditional Chinese Medicine (Jinan, China) during three epidemic windows: May–December 2021, November–December 2023, and October–December 2024. A total of 798 throat swab specimens were collected (338 in 2021, 174 in 2023, and 286 in 2024). Sampling windows were not identical across years. In 2021, microbiological surveillance was implemented over a broader May–December interval while the program was being established during the early post–COVID-19 control phase, under China’s normalized prevention-and-control and dynamic zero-COVID strategy. In 2023 and 2024, we intentionally focused culture and sequencing on shorter November–December and October–December windows that coincided with the local post-COVID resurgence of M. pneumoniae seen in hospital and regional surveillance.

The diagnosis of M. pneumoniae pneumonia (MPP) was based on clinical symptoms, lung imaging, and laboratory testing. Pneumonia severity was classified according to the Chinese National Health Commission’s Guidelines for CAP in Children (2019), which define severe pneumonia as the presence of at least one of the following: hypoxemia (oxygen saturation <92% on room air) and/or marked respiratory distress (age-adjusted tachypnea, chest wall indrawing, nasal flaring, grunting, or intermittent apnea); extensive radiographic involvement (consolidation affecting ≥2/3 of a lung, multilobar disease, or complications such as pleural effusion, pneumothorax, atelectasis, necrotizing pneumonia, or lung abscess); or serious extrapulmonary complications including sepsis, septic shock, or acute respiratory failure requiring intensive care support (National Health Commission of the People’s Republic of China, and State Administration of Traditional Chinese Medicine, 2019). The denominator for year-specific proportions (e.g., severe disease rate) was the number of children enrolled with clinically diagnosed MPP in that calendar year, which equals the number of study participants with specimens collected in that year.

The studies involving human participants were reviewed and approved by the Institutional Review Board of the Shandong Provincial Center for Disease Control and Prevention, Jinan, China (Approval No. 2021-24). The studies were conducted in accordance with local legislation and institutional requirements. Given the retrospective design using de-identified residual clinical specimens and metadata, written informed consent was waived by the IRB in accordance with national research-ethics guidelines.

All throat swab specimens were initially screened for Mycoplasma pneumoniae using a commercial real-time quantitative PCR (qPCR) kit (Shanghai BioGerm Medical Technology Co., Ltd.), following the manufacturer’s instructions. Specimens that tested positive were subsequently genotyped for P1–1 and P1–2 using a duplex real-time PCR assay as previously described (Zhao et al., 2015). Direct multilocus sequence typing (MLST) was then attempted for all qPCR-positive specimens by PCR amplification and Sanger sequencing of internal fragments of eight housekeeping genes according to the established M. pneumoniae MLST scheme (Brown et al., 2015); specimens with insufficient template DNA, typically reflected by high qPCR cycle-threshold (Ct) values, failed amplification and were not assigned sequence types. Concurrently, all PCR-positive specimens were inoculated into Mycoplasma Broth Base (CM1166; Oxoid, Thermo Fisher Scientific) supplemented with Mycoplasma Supplement G (SR0059; Oxoid) for enrichment culture, regardless of clinical severity. Upon indication of growth, aliquots were subcultured onto Mycoplasma Agar Base (CM0401; Oxoid) supplemented with Mycoplasma Supplement G (SR0059) for isolation and monitored for the appearance of typical “fried-egg colonies”.

Genomic DNA was extracted from 2 mL of log-phase Mycoplasma broth culture for each isolate using the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany). Sequencing libraries were prepared using the Nextera XT DNA Library Prep Kit and sequenced on a NextSeq 2000 platform (Illumina, San Diego, CA, USA) to generate 150 bp paired-end reads, achieving a target depth of at least 100× coverage. Raw sequencing reads were processed with Fastp (v1.0.1) to remove adapters and low-quality bases. The resulting high-quality reads were de novo assembled into draft genomes using SPAdes (v4.2.0) with default parameters. The quality of the draft assemblies was then assessed using QUAST (v5.0.2) (Jiao et al., 2025). Assemblies were retained for downstream analyses if they had a total length between 0.80 and 0.90 Mb, an N50 of at least 50 kb, no more than 50 contigs longer than 500 bp, and a mean depth of coverage ≥50×; all 227 genomes met these criteria.

The assembled genomes were used for downstream analysis. Multilocus sequence typing (MLST) was conducted using the MLST tool (v2.23.0), with sequence types (STs) assigned according to the scheme hosted on the PubMLST database. To identify macrolide resistance mutations, quality-filtered reads were mapped to the Mycoplasma pneumoniae 23S rRNA gene reference sequence (GenBank: X68422.1) using CLC Genomics Workbench (v25.0.3; Qiagen, Aarhus, Denmark). Variants were called from the alignments, and domain V of the 23S rRNA gene was screened for macrolide-resistance–associated substitutions at positions A2063G/C/T, A2064G/C, A2067G, and C2617G/A, as previously described (Wang et al., 2024).

We adopted a six-subclade framework for Mycoplasma pneumoniae (T1-1, T1-2, T1-3, T1-3R, T2-1, T2-2) based on recent whole-genome phylogenies and nomenclature proposals (Kenri et al., 2023). Within this scheme, the EC1 and EC2 epidemic clones described in China correspond to the T1-3R and T2–2 subclades, respectively (Li et al., 2024). Reads were mapped to the FH reference genome LR214945.1, and SNPs were called at the EC1/EC2-defining sites reported by Li et al. Isolates matching the EC1 or EC2 SNP profiles were classified as EC1 or EC2. EC status was then combined with P1 type and phylogenetic position to assign P1-type 1 EC1 isolates to the T1-3R (EC1) clade, P1-type 1 non-EC1 isolates to T1-3, and P1-type 2 EC2 isolates to T2-2 (EC2).

Reads were mapped to M. pneumoniae M129 (NC_000912) using Snippy (v4.6.0) with default quality filters (minimum site depth 10×, minimum mapping quality 60 and minimum base quality 13) to call SNPs and indels and generate the core alignment; core SNPs were extracted with SNP-sites, and recombination was detected and masked with Gubbins (v3.4.3) (Croucher et al., 2015). Pairwise SNP distances were computed with snp-dists, and a maximum-likelihood tree was inferred with RAxML (v8.2.12) under GTRGAMMA using rapid bootstrapping (-f a, 100 replicates). Transmission clusters were explored via a genomic network linking isolates at ≤3 SNPs (putative recent transmission) and, for broader connectivity, ≤11 SNPs, following prior outbreak work (with thresholds evaluated in sensitivity analyses) (Jiao et al., 2025). Trees were visualized in iTOL as cladograms with branch lengths ignored (Figures 1, 2) for improved readability, while all distance-based analyses were performed on the full branch-length trees.

Draft genomes were re-annotated with Prokka (v1.14.6) using -{{-}}-genus Mycoplasmoides -{{-}}-species pneumoniae -{{-}}-gcode 4 (translation Table 1 is appropriate for Mycoplasmoides/Mycoplasma) to generate GFF and protein FASTA files. Orthologous clusters and the gene presence/absence matrix were inferred with Panaroo (-i./*.gff -o panaroo_out -{{-}}-clean-mode moderate -a core -{{-}}-aligner mafft -t 4 -{{-}}-remove-invalid-genes), using the default clustering thresholds (minimum amino-acid identity 0.98, protein family identity threshold 0.70 and length difference cutoff 0.98) and a core-genome threshold of 0.95. Cluster-level association between clades (T1-3R (EC1) vs T2-2 (EC2)) and gene presence was tested in R (v4.5.0) using Fisher’s exact test; clusters with nominal p < 0.05 were carried forward to Benjamini–Hochberg false discovery rate correction, and those with adjusted q < 0.05 were considered significant. Functional annotation of cluster member proteins was performed with eggNOG-mapper (v2.1.13) (DIAMOND mode) against the official eggNOG v5 reference, summarizing COG categories at the cluster level; clusters without any COG assignment were labeled “Unannotated”, whereas clusters assigned to COG category S were labeled “Function unknown”.

Predicted proteins were searched against the Virulence Factor Database (VFDB; http://www.mgc.ac.cn/VFs/) using BLASTp. Matches were accepted only if sequence identity ≥85% and query coverage ≥85%, retaining the top hit per query. Calls were summarized as a binary presence/absence matrix for downstream analyses. For key virulence determinants, we refer to proteins by their conventional names and indicate the corresponding M129 locus tags at first mention.

Antimicrobial susceptibility was determined by broth microdilution and interpreted per CLSI M43-A (macrolides with defined S/R categories; tetracycline and levofloxacin with susceptible-only criteria) (Clinical and Laboratory Standards Institute (CLSI), 2011).

Continuous variables were summarized as medians (range) and compared across years (2021, 2023, 2024) using the Kruskal–Wallis test. Categorical variables were presented as counts (%) and analyzed with Pearson’s chi-square or Fisher’s exact test; trends over time were assessed where applicable. Pan-genome associations used Fisher’s exact test with Benjamini–Hochberg FDR control; clusters with nominal p<0.05 and adjusted q<0.05 were considered significant. Two-sided p<0.05 was considered statistically significant. Analyses were performed in R (version 4.5.0).

From 2021 to 2024, 798 children with clinically diagnosed MPP were enrolled. The PCR confirmation rate increased from 47.3% (160/338) in 2021 to 75.9% (217/286) in 2024, with an overall rate of 59.8% (477/798) (Table 2). Among the 477 PCR-positive specimens, 227 Mycoplasma pneumoniae isolates were cultured and subjected to genomic analysis.

Clinical severity of pediatric MPP increased in the post-COVID era despite stable demographics (median age ~7.5 years and similar sex ratios). Across the three annual cohorts, peak temperature and length of stay rose stepwise (38.5→38.8→39.0 °C; 7.0→8.5→9.5 days; Table 3). Systemic inflammatory burden increased as well: CRP surged in 2023 and remained above 2021 in 2024 (20.1→47.8→34.7 mg/L), while LDH peaked in 2023 and declined but stayed slightly higher than baseline in 2024 (338→464→349 U/L). Coagulation activity rose (D-dimer 1.05→1.55→1.70 µg/mL), and PCT peaked in 2023 before declining in 2024 yet remained above 2021 levels (1.9→2.9→2.3 ng/mL). Neutrophil fractions were stable, and WBC counts showed no consistent upward trend. Concordantly, the proportion of severe pneumonia increased year-on-year from 7.4% (25/338) in 2021 to 14.9% (26/174) in 2023 and 19.9% (57/286) in 2024 (Table 3), indicating a measurable shift toward more severe disease in recent epidemics.

At the specimen level, direct P1/MLST typing of PCR-positive throat swabs revealed a marked shift in genotype composition over time. In 2021, P1-2/ST14 and P1-1/ST3 were present at comparable frequencies (48.7% vs 44.5%), with a small contribution from ST30 (5.9%), whereas by 2023 and 2024 ST3 accounted for 65.8% and 80.1% of typed specimens and ST14 declined to 30.6% and 13.9%, respectively (Table 4). Untypeable P1–1 and P1–2 profiles remained uncommon in all three years, and genotype distributions are reported for specimens with successful P1/MLST typing.

To further characterize genotype dynamics during the resurgence, we performed P1 and MLST genotyping on all 227 cultured isolates. The isolate-based analysis showed a rapid shift in predominance (Figure 3, Table 1). In 2021, the Mycoplasma pneumoniae population was mixed, with P1-1 (predominantly ST3) and P1-2 (exclusively ST14) circulating at comparable frequencies (48.8% vs 51.2%). By 2023, P1-1/ST3 had become predominant (63.6%), reaching 84.0% (63/75) in 2024. Over the same period, P1-2/ST14 declined (33.3% in 2023 to 9.3% in 2024), and ST30 was not detected after 2021. A small number of untypeable isolates were observed (P1-1: 2 in 2023 and 3 in 2024; P1-2: 2 in 2024). Together with the specimen-level findings, these patterns support clonal expansion of the P1-1/ST3 lineage during the 2023–2024 epidemic in Jinan. Using the EC1/EC2-specific SNP panel, all 141 ST3 isolates (all P1-type 1) were classified as EC1, whereas all 73 ST14 isolates belonged to EC2. Thus, the post-COVID resurgence in Jinan reflects expansion of a P1-1/ST3 EC1/T1-3R lineage against a background of declining P1-2/ST14 EC2/T2–2 lineage.

To investigate the genomic epidemiology of the recent M. pneumoniae outbreak, we built a maximum-likelihood phylogeny from core-genome single-nucleotide polymorphisms (SNPs) for all 227 isolates. The core-genome tree resolved two major clusters that corresponded to the P1-type 1/T1-3R (EC1) and P1-type 2/T2-2 (EC2) lineages, in close agreement with P1 genotypes and EC1/EC2 clone assignments (Figure 1). Isolates from 2021 were scattered across both lineages, including a small subset of P1-type 1 isolates lacking the EC1/EC2 signatures that formed a T1-3-like branch, whereas almost all isolates from 2023–2024 grouped within the T1-3R (EC1) lineage and formed a compact, star-like cluster consistent with rapid recent expansion of a single clone. Notably, although both T1–3 and T1-3R subclades have been described in global datasets, all 141 ST3 genomes in this study carried the EC1 SNP signature and fell within T1-3R; we did not observe any ST3 isolates assigned to T1-3. The only non-EC1 type 1 genomes were six isolates sampled in 2021 (five ST30 and one MLST-untyped), forming a small T1-3-like branch that was not seen in subsequent years.

Pairwise SNP-distance profiles quantitatively mirrored this pattern (Figure 4). In 2021, a bimodal distribution indicated co-circulation of two genetically distant lineages; by 2024, the distribution collapsed to a single narrow peak with most pairs differing by <30 SNPs, consistent with a highly clonal population. Taken together, these data indicate clonal replacement, with a relatively homogeneous T1-3R (EC1) lineage displacing previously co-circulating, more diverse lineages including T1-3/T1-3R and T2-2. In a broader context, a phylogeny incorporating 308 global strains placed the Jinan outbreak cluster closest to contemporaneous isolates from other Chinese cities (e.g., Beijing), supporting regional transmission rather than a novel introduction (Figure 2).

Pairwise SNP–based networks showed strong clustering (Figure 5). Using a strict threshold of ≤3 SNPs, the network broke into many small components (42 clusters) with numerous singletons (n=75), consistent with short, localized chains. At a more permissive but still conservative threshold of ≤11 SNPs, isolates merged into 16 clusters, including one dominant supercluster containing 168 isolates (74% of all samples). This pattern indicates that most isolates belonged to a highly homogeneous clone that was widely distributed in the population. In a slowly evolving, low-diversity pathogen such as M. pneumoniae, such large connected components at a ≤11-SNP threshold primarily reflect limited genomic diversity and recent clonal expansion, and do not imply direct epidemiologic linkage among all connected cases.

A presence/absence pan-GWAS identified 37 clade-associated gene clusters (FDR q<0.05), with 19 enriched in T1-3R (EC1) and 18 in T2-2 (EC2). Mapping these clusters to COG categories revealed clear functional skews (Figure 6): T1-3R showed a higher share of clusters annotated to “Replication, recombination and repair” and uniquely contained clusters in “Amino acid transport and metabolism.” Both clades contributed one cluster in “Defense mechanisms,” whereas a single “Signal transduction mechanisms” cluster was observed only in T2-2. Notably, a large fraction of significant clusters in both clades were “Unannotated” (i.e., without COG assignment), with a higher proportion in T1-3R (EC1). In contrast, “Function unknown” (COG assigned but with unknown function) was more represented in T2-2. Together, these patterns point to clade-specific differences in genome maintenance and metabolic modules while also highlighting the substantial fraction of genes that remain poorly characterized in M. pneumoniae.

Next, we surveyed virulence determinants across the 227 genomes using the Virulence Factor Database (VFDB), a curated resource widely applied for genome-scale virulence profiling. Across the 227 Mycoplasma pneumoniae genomes, HMW1 (MPN447), HMW2 (MPN310), HMW3 (MPN452), P1 (MPN141), P40/P90 (MPN142), P30 (MPN453), P65 (MPN309), and CARDS toxin (MPN372) were detected in all isolates.

All 227 isolates harbored the 23S rRNA A2063G substitution, the canonical marker of macrolide resistance. No other known macrolide-resistance substitutions in domain V (A2063C, A2063T, A2064C, A2067G, C2617A or C2617G) were observed. In line with this genotype, macrolide MICs were uniformly high—erythromycin >16 µg/mL and azithromycin 8–>16 µg/mL (Table 5)—well above the CLSI M43-A susceptible cutoff of 0.5 µg/mL. By contrast, tetracycline and levofloxacin MICs stayed within the M43-A susceptible range (≤2 and ≤1 µg/mL, respectively); no isolate exceeded these thresholds, indicating preserved in vitro activity in this cohort.