Single-bacterial cell insights into mechanisms of ceftriaxone resistance in Neisseria subflava

Domain: biology
Published: 2026-03-27
Editorial outcome: minor_revision

Abstract
The contribution of airway pathobionts to chronic respiratory disease is increasingly recognized, yet the evolutionary processes that shift commensals to pathogens remain poorly understood. Here we investigate how antibiotic pressure drives adaptation in Neisseria subflava, a common airway commensal associated with bronchiectasis. Using serial passage under ceftriaxone exposure, we observe a >300-fold increase in resistance, accompanied by enhanced biofilm formation and genetic reprogramming. Whole-genome sequencing reveals recurrent mutations in the adhesin gene ataA, while single-cell transcriptomics identifies six functionally distinct clusters indicating adaptive programs in growth, metal homeostasis, oxidative stress, and cell-wall remodeling. Notably, biofilm integrity is maintained through compensatory upregulation of comP and bamE, which promotes phagocytic evasion and resistance in experimentally evolved strains and clinical isolates. Iron availability further stabilizes biofilm and modulates antibiotic tolerance, underscoring metal homeostasis as a contributory adaptive axis. Together, these findings reveal a multifaceted strategy by which N. subflava exploits antibiotic selection to transition towards pathogenicity. By integrating experimental evolution with single-cell resolution, we establish a framework for understanding the commensal-to-pathobiont transition, with broad implications for the airway microbiome and antimicrobial resistance in chronic respiratory disease.

Introduction
The dynamic interplay between host and microbiome has reshaped our understanding of airway pathology, highlighting the microbiome as a critical determinant of respiratory health and disease1–3. Advances in high-throughput molecular profiling and computation have enabled holistic analyses of airway microbial communities2,4, uncovering novel mechanisms of pathogenesis that extend beyond canonical pathogens5–11. Within this framework, the concept of pathobionts, microorganisms that are typically commensal but acquire pathogenic potential under specific conditions, has emerged as a central theme in chronic airway disease9,11–14.

Bronchiectasis exemplifies this complexity. It is a chronic inflammatory disorder of the airways, characterized by recurrent infection, het erogeneous clinical presentation, and notable geographic variation. Large-scale microbiome studies implicate both bacteriome and mycobiome composition, interkingdom interactions, host immune responses, and antimicrobial resistance determinants in its pathophysiology15–19. Notably members of the Neisseriaceae, particularly Neisseria subflava, have recently been associated with bronchiectasis in Asian populations13,20–22. Long regarded as a commensal, N. subflava can disrupt epithelial barriers, impair trans-epithelial electrical resistance, and provoke inflammation, thereby driving airway remodelling and immunopathology distinct from established pathogens such as Pseudomonas aeruginosa18. Related observations in Chronic Obstructive Pulmonary Disease (COPD) further highlighting the complexity of pathobiont biology and the influence of the environment on their pathogenesis23.

Understanding the commensal-to-pathobiont transition is therefore a critical challenge in microbiology10,18. While Neisseria meningitidis illustrates the devastating potential of such shifts in the upper airway, the mechanisms underlying more transitions in organisms such as N. subflava remain unresolved24–26. Evolutionary pressures within the host, particularly antibiotic exposure, are hypothesized to accelerate adaptation and virulence in commensals, paradoxically, fostering dysbiosis and resistance while treating acute infections10,27,28. Yet, how airway commensals adapt under antibiotic pressure at the genetic and functional levels remains poorly defined29.

Here, we address this knowledge gap by combining experimental evolution with single-bacterial cell transcriptomics to investigate the adaptive responses of N. subflava under ceftriaxone selection. We identify convergent genomic and transcriptional changes involving biofilm maintenance, iron-responsive programs, and cell-wall remodeling that together support increased antibiotic tolerance features consistent with pathogenic potential. By delineating the molecular strategies underlying the commensal-to-pathobiont transition, this work provides new insights into the evolution of the airway microbiome and its implications for chronic respiratory disease.

Results
Experimental evolution of antibiotic resistance and biofilm formation in Neisseria subflava
Neisseria subflava, a commensal of the oral cavity and upper airway, has limited documented antibiotic resistance30,31. To establish a baseline antibiotic susceptibility, we performed dose-response assays using seven antibiotic classes on the clinical isolate SG-KL01 (SG001) obtained from a bronchiectasis patient with a high N. subflava burden18 (Supplementary Fig. 1a-k). Due to the absence of established clinical breakpoints and specific PK/PD parameters for N. subflava, minimum inhibitory concentration (MIC) values served as the primary reference. Initial profiling showed relatively high MIC values (> 0.9 μg/ml) for gentamicin, tetracycline, kanamycin, and ofloxacin, whereas extended-spectrum cephalosporins (ESCs), such as ceftriaxone (Cef) and cefotaxime (CTX), exhibited consistently low MICs, confirming their clinical utility against Neisseriaceae infections (Supplementary Table 1).

Ceftriaxone is a third-generation cephalosporin widely used in respiratory infections and a key agent in the treatment of Neisseriaceae-associated infections32,33. Despite its broad-spectrum activity, serial ceftriaxone exposure of three independent SG001 lineages (Ca, Cb, Cc) over 26 generations (Fig. 1a) resulted in a 300-fold increase in ceftriaxone MICs (> 5 μg/ml versus 0.02 μg/ml in the parental strain) (Fig. 1b). These evolved lineages also displayed cross-resistance to other β-lactams (ampicillin and penicillin G) while retaining sensitivity to antibiotics outside this class (Fig. 1c, Supplementary Fig. 2a-h, Supplementary Table 2), indicating β-lactam-specific resistance development. Concomitantly, all three C26 lineages (C26-a, C26-b, and C26-c) exhibited pronounced colony wrinkling and increased biofilm formation relative to both wild-type and 9th-generation strains (C9-a, C9-b, and C9-c) as confirmed by Congo Red and crystal violet assays (Fig. 1d). Resistance was evident as early as generation 9 (C9-a, C9-b, and C9-c) with MICs of > 0.08 μg/ml in these strains (Supplementary Fig. 2i), suggesting that biofilm development and resistance co-evolve under antibiotic resistance pressure34–36.

Systems-level analysis of ceftriaxone-resistant Neisseria subflava strains
To investigate the molecular mechanisms underlying resistance, we performed long-read genome sequencing of evolved strains (Fig. 1e). All assemblies were supported by high coverage and were polished with Illumina short reads, ensuring high confidence in mutation calls. We identified predominant mutations in the ataA gene (QFG95_10090), accounting for 63–77% of the total mutations (Fig. 1f, Supplementary Fig. 2j, Supplementary Data 1). The ataA gene, involved in bacterial adhesion and biofilm formation37, harbored diverse mutational events, including nonsynonymous and synonymous single-nucleotide variations (SNVs), as well as frameshift insertions and deletions. Mutations in other genes, including mltA (QFG95_07970) and QFG95_05805, indicated parallel genomic adaptations under antibiotic pressure.

To confirm independent evolutionary trajectories, we analyzed a SNP-based dendrogram, which supported related but distinct lineages among the replicates (Supplementary Fig. 3a). Because SNPs, especially synonymous ones, can recur under similar selective pressures, we further profiled structural variants (SVs) in C9 lineages using high-resolution long-read sequencing. Large insertions, deletions, duplications, and complex rearrangements were hierarchically clustered by genomic position and variant types, revealing lineage-specific SV patterns unique to each replicate population (C9-a, C9-b, and C9-c) (Supplementary Data 2, Supplementary Fig. 3b-c). Such SV patterns are highly unlikely to arise from parallel evolutionary processes, unlike SNPs, thereby supporting independent evolutionary trajectories and excluding the possibility of cross-contamination among replicates38–40. Consequently, the observed convergence in SNP patterns across the lineages thus likely reflected shared selective pressures rather than technical artifacts.

Phenotypically, resistant strains demonstrated enhanced, rather than impaired, biofilm formation, suggesting the presence of compensatory regulatory pathways beyond ataA mutations. Bulk transcriptomic profiling supported this, revealing upregulation of genes involved in biofilm development, pilus synthesis (comP), iron metabolism (pubB), pyrimidine biosynthesis (pyrDA), carbohydrate metabolism (frk2), and ataA itself in the evolved strain compared to wild-type (Fig. 1g, Supplementary Data 3). Functional gene ontology analyses indicated the coordinated induction of genes involved in cell adhesion and two-component regulatory systems (Supplementary Fig. 2k, Supplementary Data 4), which are essential for biofilm formation and environmental stress responses. Importantly, bulk RNA-seq of the additional evolved lineages (C26-b and C26-c) revealed transcriptional profiles highly consistent with those of C26-a (Supplementary Fig. 3d), reinforcing a common transcriptomic rewiring across independently evolved replicates. Notably, in addition to comP and bamE, several genes that are directly involved in biofilm matrix construction, such as QFG95_07045 (murein-L,D-transpeptidase), murD, milC, murE, and pgpA, were upregulated across resistant lineages. The metal-ion transporter mntB (Mn/Fe uptake) was likewise increased, consistent with enhanced metal homeostasis in the resistant populations. Additionally, enrichment analysis further confirmed convergence on biofilm-formation programs in C26-b and C26-c (Supplementary Fig. 3e). Although the composition of upregulated genes differed across lineages, numerous loci involved in biofilm formation, iron metabolism, and cell-envelope biogenesis were recurrently elevated, albeit with lineage-specific magnitudes. Together with our genomic analyses, these results indicate that sustained antibiotic exposure drives complex adaptive responses in N. subflava, where genomic mutations and transcriptomic rewiring converge to promote antimicrobial resistance. The convergence of biofilm reinforcement and envelope remodeling highlights the plasticity of this commensal species and its capacity to transition toward antibiotic-resistant pathobiont behavior.

Bacterial heterogeneity drives the adaptive evolution of antibiotic resistance in Neisseria subflava
To explore biofilm-related gene expression and functional diversity at single-cell resolution, we performed single-cell RNA sequencing (scRNA-seq) on N. subflava cultured under planktonic (shaking) conditions (Fig. 2a). This approach enabled comparison of the transcriptional states of wild-type and ceftriaxone-resistant (C26) strains and highlighted population heterogeneity during antibiotic adaptation. Recent evidence indicates that biofilm and stress-response genes can be activated under planktonic conditions in response to environmental stress or sub-inhibitory antibiotic concentrations, producing a “primed” transcriptional state that precedes biofilm formation41–43. Consistent with this, our data revealed pronounced activation of gene modules associated with antibiotic resistance, with upregulation of biofilm-associated genes in the evolved strain, even under planktonic growth conditions.

We successfully captured single-cell transcriptomes from 9,661 N. subflava cells (3,542 wild-type and 6,119 C26) with comparable sequencing depths across the two groups. Wild-type bacteria averaged 96 ± 126 reads (SD) per cell and expressed 87 ± 43 genes (SD) per cell, while resistant bacteria had similar average reads per cell (194 ± 113) and slightly more genes per cell (95 ± 43) (Supplementary Fig. 4a).

UMAP visualization resolved six transcriptionally distinct populations with differential representation across conditions (Fig. 2b-c, Supplementary Data 5). WT bacteria were enriched in clusters 1 (11.5%), 2 (48.1%), and 3 (34.6%), whereas resistant bacteria primarily occupied clusters 1 (34.7%), 4 (34.4%), 5 (16.3%), and 6 (7.9%). Cluster 1, the most abundant shared population, was characterized by outer membrane porin (porB) and ribosomal genes (rplC, rpsC). Clusters 2 and 3, predominantly found in WT bacteria, were marked by genes involved in RNA processing and protein synthesis (tuf, rhlE), reflecting active growth and metabolic activity. Cluster 3 uniquely expressed lolA, adhesins, and chaperones (dnaK, groL), which play critical roles in outer membrane biogenesis and stress response.

In contrast, resistant-associated clusters exhibited divergent gene expression linked to ceftriaxone resistance. Clusters 4 and 5 shared transcriptional similarities and both upregulated comP (pilin formation), bamE (outer membrane protein assembly), bfrA and bfrB (bacterioferritin linked to oxidative stress mitigation)44, and nirK (a nitrite reductase for detoxification)45. Cluster 5 displayed enhanced metabolic capacity (glyA, acn B, purl, leuS), suggesting an adaptive rewiring of central metabolism to sustain energy demands under antibiotic stress. Lastly, cluster 6 was characterized by sulfur metabolism genes (cysP, cysI, cysJ, cysK), which potentially support alternative redox balance and cell wall modification, enabling evasion of β-lactam targeting.

In concordance with the bulk transcriptomic analysis, scRNA-seq profiling of C26-resistant N. subflava strains revealed key molecular adaptations that confer ceftriaxone resistance, with marked overexpression of comP, bamE, nirK, ataA, as well as carboxymuconolactone decarboxylase family proteins, transferrin-binding proteins and murein-L,D-transpeptidase (Fig. 2d-e, Supplementary Data 6). Although ataA carries a nonsense mutation predicted to produce a truncated, non-functional protein, it nonetheless showed transcriptional upregulation in C26 compared with the WT strain. This suggests a potential compensatory or autoregulatory role in biofilm dysregulation under antibiotic selection. Collectively, these data identify comP and bamE as key genes enriched in resistant-associated clusters, strongly implicating them in biofilm-associated ceftriaxone adaptation.

Lineage-specific functional adaptations for sustaining resistance
Gene Ontology (GO) enrichment analysis of the cluster-specific marker genes demonstrated that bacterial clusters of the resistant strain possess enhanced oxidoreductase activity, which may neutralize the oxidative stress induced by ceftriaxone while simultaneously activating stress response regulons that promote survival. Interestingly, the peptidoglycan biosynthesis pathways were enriched in clusters 3 and 5, but exhibited divergent genetic programs, with cluster 5 exhibiting important antibiotic adaptations including metal ion binding and enhanced DNA repair capacity, both relevant to countering ceftriaxone-induced damage to cell wall precursors and DNA (Fig. 3a, Supplementary Data 7). KEGG analysis indicated heightened macronutrient metabolism in clusters 4 and 5, supporting the energetic demands of efflux activity and cell-wall repair under cephalosporin pressure (Fig. 3b, Supplementary Data 7). This metabolic rewiring mirrors resistance strategies observed in pathogenic Neisseria species, where increased metabolic flux through the TCA cycle and oxidative phosphorylation supports β-lactam resistance46. The specialized functions of clusters 4 and 5, particularly their overexpression of denitrification genes (nirK, ncd2, norB), competence and adherence factors (comP, bamE), SOS and stress response (sodB, cspA, rnr, recQ), and metal ion acquisition systems (bfrA, bfrB, lutA1, pqqL), collectively reinforce envelope integrity, limit oxidative damage, and maintain metal homeostasis under ceftriaxone exposure (Fig. 3c). By contrast, cluster 6 adopted a metabolically conservative program characterized by sulfate assimilation (cysP, cysI, cysJ, cysK), redox activity (nuoF, nuoD, nqrB) and selected efflux transporters (tet(M), acrB, potA, tet35) (Fig. 3b-c). As sulfate assimilation precedes cysteine biosynthesis, this transcriptional profile indicates elevated cysteine production to bolster redox balance and antioxidant defenses during β-lactam-induced envelope stress.

Although bfrA and bfrB genes primarily function as iron-storage proteins, they may secondarily contribute to enhancing biofilm stability and stress tolerance47–49. Inductively Coupled Plasma Mass Spectrometry (ICP-MS) confirmed increased intracellular Fe levels in resistant versus wild-type bacteria, consistent with elevated expression of iron-related genes, while Mg levels remained unchanged (Supplementary Fig. 4b). We also compared the growth kinetics of wild-type and C26 strains. The C26 strain also showed accelerated growth rates and higher stationary-phase density, consistent with enhanced metabolism and proliferative potential (Supplementary Fig. 4c). Additionally, we established a bfrSA-gfp genomic reporter system to directly assess iron-dependent gene regulation. Upon iron supplementation (50 μM FeCl₃), we observed robust induction of GFP fluorescence driven by the bfrA promoter, confirming enhanced iron-responsive gene activity in the evolved C26 strain (Supplementary Fig. 4d). These results provide direct support for the adaptive role of metal homeostasis identified by single-cell analyses, reinforcing the crucial role of iron regulation in cephalosporin resistance.

Among these modules, comP and bamE showed peak expression in the resistance-associated cluster 5, highlighting their central roles within the dominant adaptive lineage and linking biofilm features with antibiotic tolerance. To assess whether planktonic culture conditions could obscure biofilm-specific expression, RT-qPCR was performed under both planktonic (shaking) and sessile (biofilm-promoting) conditions. The evolved C26 strain maintained similarly elevated expression of both genes under planktonic and sessile conditions, whereas wild-type cells showed no significant difference (Supplementary Fig. 4e). Thus, the resistant lineage maintains a constitutively “primed” biofilm transcriptome across growth modes. At the population level, comP and bamE are the principal genetic drivers linking biofilm enhancement to ceftriaxone resistance. ComP encodes a DNA uptake minor pilin that mediates natural competence and facilitates incorporation of extracellular DNA (eDNA) into the matrix, a structural component that reinforces Neisseria biofilms and stabilizes surface communities under antibiotic stress50–52. BamE, a lipoprotein of the β-barrel assembly machinery (BAM), ensures efficient outer-membrane protein assembly, preserving the integrity and architecture of the biofilm surface layer and reducing β-lactam permeability53–55. Together, comP-driven eDNA enrichment and bamE-mediated envelope stabilization establish a biofilm phenotype that shields cells from cephalosporin penetration and maintains a microenvironment conducive to horizontal gene exchange.

Single-cell RNA-seq revealed functional heterogeneity underpinning ceftriaxone resistance in N. subflava, which confers a survival advantage by enabling specialized roles, enhancing resilience, and supporting long-term adaptation29,56. This resistance is centred on a biofilm-oriented programme driven by comP-mediated extracellular DNA enrichment and bamE-dependent outer-membrane stabilisation. This programme converges with the functional enrichment profiles of clusters 4-6, where stress-response activation, cell-wall remodelling and metal-ion homeostasis form an integrated network that maintains biofilm integrity under sustained β-lactam pressure. These findings support a model in which resistance evolution is intrinsically linked to biofilm optimisation, providing immediate protection from antibiotic challenge while preserving a genetic reservoir for long-term adaptation.

Elevated expression of comP and bamE drives antibiotic resistance and immune evasion
Building on these observations, we next dissected the functional contributions of comP and bamE as the principal biofilm-linked drivers, alongside bfrA and bfrB as bacterioferritin-associated iron homeostasis modules. To assess gene expression at single‑cell resolution, we integrated chromosomally integrated fluorescent reporters (mCherry for comP and GFP for bamE) in both wild-type and C26 strains. Fluorescence microscopy revealed markedly increased reporter intensity in the evolved C26 strain, with a larger fraction of cells showing high comP and bamE expression compared with the wild type (Fig. 4a).

MIC assays confirmed that the overexpression of comP and bamE enhanced ceftriaxone resistance, with MIC values of 0.2 μg/ml and 0.1 μg/ml, respectively, representing a 5-10-fold increase compared to vector controls. Conversely, deletion of comP or bamE reduced MICs by approximately two-fold, confirming their essential roles in resistance (Fig. 4b). Both genes also promoted biofilm maturation. Their overexpression enhanced biofilm mass, whereas deletion markedly impaired biofilm formation (Fig. 4c, Supplementary Fig. 5a). Overexpression of bfrA and bfrB did not affect ceftriaxone resistance and had a minor impact on biofilm mass, suggesting their limited role. Macrophage phagocytosis assays revealed enhanced resistance to phagocytic clearance in strains overexpressing comP and bamE, consistent with their roles in biofilm formation and immune evasion (Fig. 4d).

The clinical relevance of these findings was supported by analysis of patient cohorts. In bronchiectasis patients from Shenzhen, China, N. subflava was detected in 18.4% of samples (n = 38) (Supplementary Fig. 5c), with carriers generally being older and showing lower leukocyte counts (Supplementary Fig. 5d). In an independent cohort, nine clinical isolates were examined by multilocus sequence typing (MLST), which grouped them into three sequence types: ST-a (n = 6), ST-b (n = 1), and ST-c (n = 2). ST-a isolates exhibited significant upregulation of comP and bamE, which correlated with enhanced biofilm formation and higher ceftriaxone MICs (Fig. 4e-f; Supplementary Fig. 5b). Whole-genome sequencing further revealed frequent ataA mutations across clinical isolates, particularly in ST-a isolates, suggesting that comP and bamE upregulation may compensate for impaired ataA function. This suggests comP and bamE upregulation potentially compensates for impaired ataA function, promoting adaptation under antibiotic pressure. Consistent with this, ST-a isolates formed more robust biofilms than ST-b or ST-c strains, mirroring the phenotypes observed in experimentally evolved resistant lineages (Fig. 4g).

Taken together, our results demonstrate direct involvement of comP and bamE in ceftriaxone resistance, biofilm formation, and immune evasion in N. subflava. By contrast, iron-storage genes bfrA and bfrB play ancillary roles in this emerging pathobiont.

Iron enhances biofilm formation and antibiotic resistance in ceftriaxone-evolved N. subflava
Although the precise roles of bfrA and bfrB were not resolved functionally, their upregulation in C26 antibiotic-resistant strains suggests that iron homeostasis could modulate adaptation under antibiotic stress. We therefore examined the effect of iron (Fe³⁺) supplementation on bacterial survival and biofilm formation in wild-type (WT) and ceftriaxone-resistant (C26) strains.

In colony-forming unit (CFU) assays, Fe³⁺ supplementation did not alter the viability of wild-type strain at its MIC (0.02 μg/ml) but markedly improved survival of the C26 strain at its higher MIC (6 µg/ml) (Fig. 5a). Fe³⁺ supplementation had no impact on WT viability at its MIC (0.02 μg/ml), indicating iron availability does not influence WT growth at low antibiotic pressure. In contrast, Fe³⁺ markedly enhanced viability of the C26 strain at its higher MIC (6 μg/ml), significantly reducing ceftriaxone susceptibility. In the absence of iron supplementation, the viability of the C26 strain sharply declined, underscoring the critical role of iron in sustaining resistance under antibiotic exposure (Fig. 5b).

Biofilm formation showed that Fe³⁺ supplementation increased biofilm biomass in both strains under antibiotic-free conditions, with the evolved C26 strain exhibiting greater biofilm formation compared to WT. Under ceftriaxone treatment at MIC levels, biofilm biomass in the WT strain remained unaffected by iron supplementation, whereas the C26 strain showed significantly increased biofilm biomass, indicative of reduced antibiotic susceptibility (Fig. 5c-d).

Further assessment of bacterial survival across a gradient of antibiotic pressures highlighted strain-specific advantages provided by Fe³⁺ supplementation. At ceftriaxone concentrations slightly above WT MIC (0.02–0.04 μg/ml), WT viability was minimal, irrespective of iron availability. In contrast, Fe³⁺ supplementation substantially enhanced the viability of the C26 strain even at higher ceftriaxone concentrations (6-8 μg/ml), further confirming iron’s role in resistance under significant antibiotic stress (Fig. 5e-f).

These findings highlight iron’s dual contribution to biofilm formation and bacterial viability under antibiotic challenge, specifically benefiting the evolved, resistant C26 strain. Although iron broadly supports biofilm development, only the antibiotic-adapted strain effectively utilizes iron-enhanced biofilms to sustain ceftriaxone tolerance. This highlights ceftriaxone resistance as a modifiable trait in antibiotic-adapted N. subflava, influenced by iron availability via supplementation or chelation.

Methods
All experiments were performed using standard microbiological techniques unless otherwise specified. A complete list of bacterial strains, plasmids, and reagents used in this study is provided in Supplementary Table 3.

Ethics statement
The study involving human clinical samples was approved by the Institutional Review Board of Shenzhen Third People's Hospital (Ethics approval number: 2024-153). Informed consent was obtained from all participants involved in the study.

Bacterial culture media
Neisseria subflava was initially recovered from glycerol stocks stored at -80°C and cultured on Columbia blood agar plates (Columbia agar supplemented with 5% defibrinated sheep blood) at 37 °C in a 5% CO₂ incubator for 14 h. Subsequently, a single colony was transferred to Brain Heart Infusion (BHI) liquid medium (Cat#024053, HKM, China) supplemented with 5% fetal bovine serum (FBS)for culture. Bacteria were incubated at 37 °C in a 5% CO₂ incubator for an additional 12–16 h to ensure optimal growth.

MIC determination (Minimum Inhibitory Concentration)
The overnight bacterial cultures were re-inoculated at a 1:1000 dilution into fresh medium and allowed to grow to the mid-log phase in a shaking incubator (cat#MQT-60R, Minquan, China) at 37 °C, 220 rpm. The bacterial solution was added to a transparent 96-well plate at a final concentration of OD 0.01 per well, and a gradient dilution of antibiotic concentration using BHI, 600 with four parallels for each concentration. After culturing the bacteria in various concentrations of antibiotics for 24 h, the absorbance at 600 nm was measured to determine the antibacterial effect. Fitting curves by comparing different concentrations of antibiotics and absorbance to obtain IC90 (Inhibitory Concentration 90). Use GraphPad Prism 8 (GraphPad Software) to calculate S-type dose response (slope variation) nonlinear fitting.

Biofilm formation assay
Overnight cultures were diluted to OD of 0.01 and the 96-well plates were cultured at 37 °C for 600 24 h, then cell plates were washed three times with PBS to remove unbound bacteria. Then biofilm was evaluated by 0.05% crystal violet (CV: 150 μl; Sangon Biotech, China) staining for 30 minutes84. Gently rinsed unbound dyes with sterile dd H O and blown dried at room temperature. 2 After that, 200 μl of 30% glacial acetic acid was added to each well. Relative biofilm formation was determined by measuring the absorbance at 550 nm with Tecan Infinity 200 PRO microplate reader. Biofilm inhibition experiments were conducted in three independent trials, with four replicates per condition.

Congo red plate assay
Columbia agar was autoclaved, cooled to 50-55 °C, mixed with 5% sterile defibrinated sheep blood, supplemented with 10 μg/mL Coomassie brilliant blue G-250 and 20 μg/mL Congo red (CR), and poured into plates to solidify. Fresh overnight bacterial culture (5 μl per spot, 5 replicates) was spotted onto the plates, which were incubated at room temperature in the dark for 72 h, and then colony staining characteristics were observed and recorded.

Growth curve assay
To assess growth kinetics of wild-type and evolved isolates, overnight cultures were initially adjusted to an OD₆₀₀ 1.0 in BTHI broth. The standardized suspensions were diluted into fresh BHI medium supplemented with 5% FBS to a starting OD₆₀₀ of approximately 0.02. Aliquots of 200 μl were transferred to flat-bottom 96-well microtiter plates (Corning), with at least three technical replicates per strain. Plates were incubated at 37 °C in TECAN Spark microplate reader, which recorded OD₆₀₀ every 30 minutes for 48 hours with continuous double orbital shaking between reads. Blank wells containing medium only were included for background correction. Growth curves were plotted and analyzed using GraphPad Prism 8. Parameters including lag phase duration, exponential growth rate, and maximal OD were extracted to compare relative growth 600 fitness between strains.

Preprocessing for single-cell sequencing
Wild-type and antibiotic-resistant strains were cultured at 37 °C with shaking at 220 rpm. Cultures were grown for approximately 6 hours until the optical density at 600 nm (OD ) reached 1.0. To 600 minimize technical artifacts and preserve in vivo transcriptomic states, all subsequent steps were performed on ice (4 °C), and all buffers were pre-chilled. The time from harvest to fixation was strictly limited to under 10 minutes for all samples. Wild-type and C26 strains were processed in parallel under identical conditions to ensure comparability. 1 ml of each bacterial culture was collected and centrifuged at 6,010 × g for 5 minutes. The supernatant was discarded, and the pellet was washed with 1 ml of phosphate-buffered saline (PBS). The washed bacterial pellet was resuspended in 1 ml of ice-cold 4% paraformaldehyde (PFA) in PBS. The suspension was placed in a rotary mixer at 4 °C and rotated for 9 hours to cross-link intracellular RNA, DNA, and proteins. After fixation, the samples were centrifuged at 10,000 × g for 10 minutes at 4 °C to remove the PFA. The pellet was washed twice with 1 ml of PBS-T (PBS containing 0.1% Tween-20 and 1 µl RNAse inhibitor). The bacterial pellet was then resuspended in 200 µl of PBS-T, with an additional 5 µl of RNAse inhibitor added. The tubes were sealed and stored at 4 °C until further use. Prior to single-cell sequencing, quality control of the samples was performed to ensure a cell count of ≥ 10 million, minimal debris, no visible clumping, and an aggregation rate of ≤ 20%.

Single-cell RNA sequencing
Single-cell encapsulation, droplet generation, and nucleic acid capture were conducted using the VITAcruizer DP400 instrument (Cat#E20000131, M20 Genomics), following the VITApilote High-Throughput Single-Bacterial Transcriptome Kit protocol (Cat#R20114124, M20 Genomics). Fixed bacterial cells underwent permeabilization and cell wall digestion, followed by reverse transcription with random primers to produce full-length cDNA. After adaptor ligation, cell barcodes and unique molecular identifiers (UMIs) on microbeads were combined with adaptor-ligated cDNA and encapsulated in droplets. Droplet encapsulation and extension generated labeled cDNA, which was then PCR-amplified. The final library, incorporating P5 and P7 adaptors, was sequenced on an Illumina NovaSeq 6000 platform.

Single-cell data preprocessing
Single-cell RNA sequencing data were processed using VITAseer V1.2 software (M20 Genomics), designed specifically for transcriptomic data generated by the VITA platform. Briefly, raw sequencing reads were trimmed using Cutadapt (Pv3.7). For each Read1, unique molecular identifiers (UMIs) and cell-specific barcodes were extracted, with acceptable barcodes defined as those with a Hamming distance ≤ 2 nt and a unique match. Read2 was aligned to the reference genome of the species of interest using STAR (v2.7.10a) under default parameters, retaining only uniquely mapped sequences for downstream analysis. Subsequent data analysis was performed using Seurat (v4.3.0.1) under default parameters unless specified otherwise. Cells expressing fewer than 10 genes were excluded, and genes present in fewer than one cell were filtered out. Data were log-normalized using NormalizeData function, followed by FindVariableFeatures function to identify highly variable features for dimensional reduction. The data were scaled with ScaleData function. Dimension reduction was performed by principal component analysis (PCA) to reduce dimensions to 8. Uniform Manifold Approximation and Projection (UMAP) was constructed using these principal components. Clustering was performed using the Louvain algorithm (FindClusters) at a resolution of 0.5, and differential gene expression was assessed with a log-fold change threshold of 0.25 and adjusted p-value of < 0.05 (FindMarkers function).

GO and KEGG enrichment analysis of DEGs
Based on the hypergeometric tests, GO and KEGG enrichment analysis of DEGs was performed to identify the significantly enriched terms using clusterProfiler package (v4. 16.0) in R (v4.5.1).

Genomic DNA extraction and genome sequencing
Neisseria subflava SG-KL01 (SG001) were grown to 109 CFU/ml, and total DNA was extracted from each 2 ml bacterial suspension using the TIANamp Bacterial DNA Kit (Tiangen Biotech, Beijing, China) according to the manufacturer’s instructions. Qualified genomic DNA was fragmented with G-tubes (Covaris, Woburn, MA, USA), followed by end-repair using the NEBNext FFPE DNA Repair Kit (M6630, NEB, USA) and NEBNext End Repair/dA-Tailing Module reagents (E7546, NEB, USA). The fragmented DNA was subsequently ligated to adapters to construct sequencing libraries using NEBNext Quick Ligation Module reagents (E6056, NEB, USA). After purification with the AMPure XP system (Beckman Coulter, Brea, CA, USA), library sequencing was performed on the Oxford Nanopore PromethION long-read sequencing platform (Oxford Nanopore Technologies, Oxford, UK) following standard protocols. For Illumina sequencing, qualified genomic DNA was first randomly sonicated, then subjected to end-repair, A-tailing, and adapter ligation using the NEBNext® Ultra™ DNA Library Prep Kit for Illumina (NEB, USA) according to the manufacturer’s protocol. DNA fragments of 300-400 bp (in length) were enriched by PCR. Finally, PCR products were purified using the AMPure XP system (Beckman Coulter, Brea, CA, USA). Library size distribution was analyzed with the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA), and quantification was performed by real-time PCR. Genome sequencing was conducted on the Illumina NovaSeq 6000 platform using paired-end technology (PE150).

Genome assembly and quality control
Raw ONT (.fast5) signals were base-called with Oxford Nanopore Guppy v6.4.6 (HAC), followed by quality control to remove reads with Q < 7 and trim adapter-containing reads, and the filtered data were used for assembly and alignment. Raw Illumina reads were processed with FASTP v0.20.0 to (i) remove adapters (trimming adapters and downstream sequence; reads < 50 nt after trimming were discarded), (ii) remove poly-A reads, (iii) remove reads with N > 10%, and (iv) remove low-quality reads where Q ≤ 20 bases constituted ≥ 50% of the read; the resulting clean reads were aligned to the assembly with BWA v0.7.17-r1188 and used for iterative polishing with Pilon v1.23 to obtain the final genome sequence. Per-base coverage was obtained with samtools depth (v1.16.1; ONT) and bedtools genomecov (v2.30.0; Illumina), and coverage breadth at thresholds (≥1×–≥200×) was derived from these per-base depths (Supplementary Data 8). Inter-assembly variant calling was performed by whole-genome alignment with MUMmer v3.1 and regional structural variation interpreted with SyRI v1.4; SNPs/indels were derived from MUMmer/nucmer alignments (parameters: --maxmatch -c 100 -b 200 -l 50), and SVs were detected with SyRI (parameters: --unic 100 -k).

Bulk RNA-seq library preparation, sequencing and analysis
For RNA extraction, overnight cultures of SG001 were diluted to an optical density (OD) at 600 nm of 0.01 and grown at 37 °C with shaking at 220 rpm for 12 h. Total RNA was extracted from whole cells with EASYspin Plus Bacterial RNA Rapid Extraction Kit (cat#RN43, Aidlab, China), according to the manufacturer’s protocol. RNA purity and quantity were evaluated using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). RNA integrity was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). After rRNA depletion, RNA was fragmented in fragmentation buffer and subjected to first- and second-strand cDNA synthesis, followed by end repair, A-tailing and adaptor ligation (with USER treatment where applicable). Libraries were PCR-amplified, purified with AMPure XP beads, and size-selected to 370-420 bp prior to final QC. Qualified libraries were sequenced on an Illumina NovaSeq 6000 to generate 150-bp paired-end (PE150) reads. Adapters and low-quality bases were trimmed with cutadapt; reads < 50 nt or with > 10% N were discarded, and reads mapping to rRNA were removed to obtain clean reads. Per-sample metrics (Raw reads, Adapter trimmed, rRNA removed, Clean rate, Q20/Q30, GC%) were summarised. Clean reads were aligned to the Neisseria subflava reference genome using STAR (default parameters), retaining uniquely mapped reads. Gene-level raw counts were obtained with featureCountsR; TPM values were additionally computed for display/visualisation only and were not used for statistical testing. Differential expression was performed on raw counts using edgeR (negative-binomial modelling with multiple-testing correction); unless otherwise specified, genes with p < 0.05 and |log fold change| ≥ 1 were considered differentially expressed, and FDR values were reported. For visualisation, the TPM matrix of DEGs was transformed as log (TPM + 1), row-wise Z-scored, and subjected to hierarchical clustering of genes and samples to generate heatmaps. 2

Bacterial strains and genomic sequencing
The clinical isolate Neisseria subflava strains (SG001–009) were isolated from patients with bronchiectasis or other chronic lung diseases in the Singapore-Kuala Lumpur region (SG-KL, Malaysia)18. Genomic DNA was extracted from overnight cultures using the TIANamp Bacteria DNA Kit (cat#GDP302-02, Tiangen, China) according to the manufacturer’s protocol. High-throughput sequencing was conducted on the Illumina NovaSeq 6000 platform, generating paired-end reads with a read length of 150 bp. Raw sequencing reads were assessed for quality using FastQC (version 0.11.9) and subsequently trimmed using Trimmomatic (version 0.39) to remove low-quality bases and adapter sequences. De novo genome assembly was performed using SPAdes (v3.15.4) under default settings, and genome annotations were generated using Prokka (v1.14.6).

ICP-MS analysis of intracellular metal content
Bacterial cultures of wild-type and evolved (C26) Neisseria subflava strains were grown to mid-log phase in BHI broth at 37 °C with shaking. Approximately 1 × 10⁹ colony-forming units (CFU) were harvested by centrifugation at 5,000 × g for 10 min at 4 °C. Washed pellets were resuspended in 200 μl of ultrapure water and subjected to acid digestion by adding 500 μl of 65% trace-metal-grade nitric acid (HNO₃) in metal-free Teflon tubes. Digestion was carried out at 95 °C for 1 hour in a closed heating block or using a microwave digestion system where available, following validated protocols for biological samples. After cooling to room temperature, the digests were diluted to 5 ml with ultrapure water and filtered through a 0.22 μm PTFE filter to remove particulates. Metal ion concentrations (Fe, Mg) were quantified by inductively coupled plasma mass spectrometry (ICP-MS, Agilent 7900) under standard operating conditions with helium collision mode. Calibration was performed using certified multi-element standards. Three biological replicates were analyzed for each strain.

RAW264.7 cell phagocytosis resistance assay
RAW264.7 murine macrophage-like cells were maintained in DMEM supplemented with 10% FBS at 37 °C in a humidified incubator with 5% CO₂. Cells were seeded at a density of 5 × 10⁵ cells per well in 12-well plates and incubated overnight to allow adhesion. Neisseria subflava strains were grown overnight in BHI medium at 37 °C with shaking. Cultures were diluted to an OD600 of 0.6 in fresh BHI medium and grown to mid-log phase. Bacterial cells were washed twice with PBS, resuspended in antibiotic-free DMEM, and added to the RAW264.7 cells at a multiplicity of infection (MOI) of 100:1. Plates were centrifuged at 300 × g for 5 minutes to synchronize infection and then incubated for 1 hour at 37 °C with 5% CO₂. Following infection, cells were treated with 100 µg/mL gentamicin for 30 min at 37 °C to eliminate extracellular bacteria after PBS washing, then washed again and lysed with 0.1% Triton X-100 in PBS. The lysates were serially diluted. Dilutions were plated on BHI agar and incubated at 37 °C for 24 h to determine colony-forming units (CFU). CFU counts were normalized to the initial RAW264.7 cell density, and results were expressed as mean ± SD from three independent experiments.

Sequence typing (ST) analysis
Multilocus sequence typing (MLST) was conducted using the MLST scheme for Neisseria species. Seven housekeeping genes (abcZ, adk, aroE, fumC, gdh, pdhC, and pgm) were amplified and aligned against the PubMLST database. The resulting ST profiles were clustered using goeBURST to determine phylogenetic relationships. SG001–006 were assigned to the ST-a group, SG007 was classified as the ST-b group, and SG008 and SG009 were grouped into the ST-c group.

RT-qPCR
Total RNA was extracted using EASYspin Plus Bacterial RNA Rapid Extraction Kit, (cat#RN43, Aidlab, China). RNA concentrations were measured using NanoDrop 2000c (Thermo Fisher, USA) and normalized among the experimental sample set. cDNA was reverse transcribed using the PrimeScriptTM RT reagent Kit with gDNA Eraser (TaKaRa Bio, Japan) according to the manufacturer’s instructions. The SYBR Green Pro Taq HS qPCR Premix (Cat#AG11701, AG, China) and reactions were run on ABI Step One Plus Real-Time PCR equipment in accordance with the manufacturer’s instructions (Applied Biosystems, Foster City, CA, USA). The levels of target genes were normalized to the levels of gene actin (Gyrb) mRNA and relative expression was determined according to the ΔΔCt method. Three candidate genes (comP, bamE, bfrA, and bfrB) were randomly selected from the transcriptome data to verify the accuracy of the transcriptome data; the primers are listed in Supplementary Table 4.

Gene Knockout in Neisseria subflava
Approximately 700 bp sequences flanking the upstream and downstream regions of the target gene were PCR-amplified from N. subflava genomic DNA, using primers listed in Supplementary Table 7. The gentamicin resistance cassette, containing its own promoter and terminator, was PCR-amplified from a plasmid template. These three fragments (upstream arm, gentamicin cassette, downstream arm) were fused into a single linear knockout cassette by overlap extension PCR and verified by agarose gel electrophoresis. For transformation, N. subflava was streaked on blood agar plates and incubated overnight at 37 °C. Fresh colonies were scraped and resuspended in ice-cold 0.3 M sucrose, washed three times, and finally resuspended in 0.3 M sucrose to adjust the suspension to an OD600 of approximately 1.5. Fifty microliters of electrocompetent cells were mixed with ~1 µg of purified knockout cassette and transferred into a pre-chilled 0.2 cm electroporation cuvette. Electroporation was performed using a Bio-Rad Gene Pulser at 2.5 kV, and cells were immediately recovered in 1 ml of pre-warmed BHI broth with 15% FBS at 37 °C with gentle shaking (110 rpm) for 2-6 h. Transformants were selected on blood agar supplemented with 15 µg/ml gentamicin and incubated at 37 °C for 36-48 h. Gentamicin-resistant colonies were screened by PCR using primers flanking the homologous arms, and correct gene replacement was further confirmed by Sanger sequencing. Colonies were screened by colony PCR with primers outside the homologous regions, and positive mutants were confirmed by Sanger sequencing.

Fluorescent Tagging via Homologous Recombination
For in situ fluorescent tagging (comP::mCherry, bamE::gfp, and bfrA::gfp), the double-crossover homologous recombination strategy described above was employed, with the fluorescent marker gene (mCherry or gfp) inserted in-frame at the target locus along with a gentamicin resistance cassette for selection. Fusion constructs consisted of relatively long homologous arms flanking the target site, along with external primers to enhance transformation efficiency. Electroporation, recovery, and selection were performed as described for gene knockout mutants. Correctly tagged clones were verified by fluorescence microscopy, PCR, and Sanger sequencing. Strains were cultured overnight on blood agar, then inoculated into BHI broth and grown at 37 °C with shaking (220 rpm) to mid-log phase (OD600 ≈ 0.5-0.8). Bacterial cells were harvested by centrifugation, washed once with PBS, and resuspended in PBS. A 5 μl aliquot was spotted onto a glass slide for microscopy. In bacterial live-cell localization experiments, the ECLIPSE Ti2-I microscope was used to capture both mCherry and GFP channel fluorescence, as well as bright-field images of bacteria. NIS-Elements Viewer 5.27 was employed for subsequent image processing. For each condition, three independent fields of view were analyzed. Primers for fluorescent reporter fusions (comP::mCherry, bamE::gfp, bfrA::gfp) are provided in Supplementary Table 6.

Fluorescence Imaging of bfrA::gfp
To assess iron-dependent expression of the bfrA::gfp reporter, single colonies were inoculated from blood agar into BHI broth and cultured at 37 °C with shaking (220 rpm) to mid-log phase (OD600 ≈ 0.5-0.8). Cultures were divided and supplemented with FeCl₃ (0 or 50 μM final concentration), followed by continued incubation for 6 hours. Bacterial cells were harvested by centrifugation, washed with PBS, and resuspended in PBS. A 5 μl aliquot was placed on a glass slide for microscopy. The ECLIPSE Ti2-I microscope was used to capture fluorescent images of bacteria, and NIS-Elements Viewer 5.27 was employed for subsequent image processing. Three independent fields of view were analyzed per group.

Cloning procedures and overexpression of genes
DNA cloning and plasmid construction were performed using standard molecular biology techniques. PCR primers were designed to incorporate BamHI and HindIII restriction sites at their 5’ ends, enabling the subsequent digestion and ligation of the target gene into the pMCS5 vector, which harbors a tetracycline resistance gene. Seamless cloning was performed according to the manufacturer’s guidelines. The recombinant plasmid was introduced into Escherichia coli DH5α and Neisseria subflava via electroporation. Positive clones were selected on 5% sheep blood agar plates supplemented with 15 µg/ml tetracycline (TC), and all plasmids were verified by Sanger sequencing. Neisseria subflava strains harboring the recombinant plasmid were cultured in medium containing 15 µg/ml TC until reaching logarithmic growth, at which point gene expression was induced with 1 mM IPTG. The overexpressed gene is in Supplementary Table 5.

Bacterial isolation from bronchoalveolar lavage fluid (BALF)
Bronchoalveolar lavage fluid (BALF) samples were collected from 38 patients with bronchiectasis undergoing fiberoptic bronchoscopy at Shenzhen Third People's Hospital, based on clinical indications. During the procedure, sterile 0.9% saline was instilled into the right middle lobe or lingular segment of the bronchus and recovered using low-pressure suction (≤ 100 mmHg). Samples were immediately placed on ice to preserve microbiome composition, then centrifuged at 3,000 × g, 4 °C for 10 min, and the pellet was resuspended in sterile PBS. To maximize bacterial recovery, aliquots were plated onto chocolate agar and sheep blood agar (37 °C, 5% CO₂, 48 h) for common respiratory pathogens, MacConkey agar (37 °C, aerobic incubation, 48 h) for Gram-negative bacteria, and CDC anaerobic agar (37 °C, anaerobic conditions, 5-7 days) for anaerobic species. After single-colony purification, individual bacterial isolates were selected for species identification using 16S rRNA sequencing and confirmed by Sanger sequencing.

Statistical analyses
All data analysis was expressed as the mean ± SD. Each experiment was repeated a minimum of three times, with statistical analysis based on biological replicates. Detailed information regarding replicates and statistical methods for each experiment can be found in the figure legends and the Source data file. Differences were considered significant when the P value was < 0.05. GraphPad Prism 8 software (GraphPad Software Inc., San Diego, CA, USA) was used to analyze the experimental data.

Data availability
These datasets are publicly accessible and provide detailed information to support data validation, reproducibility, and further analysis by the research community. The transcriptomic and single-cell RNA sequencing datasets generated in this study have been deposited in the NCBI Sequence Read Archive (SRA). The single-cell RNA-seq and bulk RNA-seq datasets for wild-type and C26-a strains are available under BioProject accession PRJNA1209901, and the whole-genome sequencing data are under PRJNA1208745. To support lineage-level validation, additional bulk RNA-seq datasets for independently evolved C26-b and C26-c strains have been deposited under BioProject accession PRJNA1292455. All data are publicly available as of the date of publication. Source Data are provided with this paper.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript “Single-bacterial cell insights into mechanisms of ceftriaxone resistance in Neisseria subflava” presents a\nmultifaceted study that aims to understand mechanisms of evolved resistance to antibiotics by Neisseria subflava. This\nbacterium has recently been implicated with airway disease and presents a unique way to track genetic changes that allow\norganisms typically existing as commensals to become pathogenic. In this study the authors use a lab-based evolution\napproach to evolve Neisseria in the presence of a relevant antibiotic in order to identify genetic changes that are\naccumulated. Triplicate cultures are used to try and uncover mutations that are repeated and are therefore deterministic. The\nstudy then checks whether evolved strains have resistance to multiple antibiotics and attempts to characterize the\nmechanisms of resistance. This characterization is first done using physiological measurements that demonstrate that the\nevolved cultures have greater biofilm producing capacity despite all cultures having mutation that inactivate a key biofilm\nregulator. Then the authors use traditional bulk RNAseq analysis to identify gene expression changes differentiating one of\nthe lineages from the parental strain. After this analysis the authors use a bacterial single cell sequencing approach to try\nand understand the behavior of the culture. This single cell RNAseq characterization is the section of the manuscript that is\nused to make many of the less substantiated and more bold claims but unfortunately this section has several major issues\nwith study design, analysis and interpretation. Importantly, the author's central claim is that the core mechanism of resistance\nis through biofilm formation, but the single cell analysis is all done on planktonic cultures grown with much shaking. In\naddition, the authors point to varied states and trajectories that are supposed to be different between the parent strain and\nevolved culture (which apparently was only done for one of the lineages) - but the changes appear to be minimal and the\ntrajectory analysis, which is notoriously a non-reliable single cell analysis method even in higher organisms that have\ndefined developmental processes, is not validated by any dynamic measurement and makes several assumptions that are\nnot founded. Furthermore, data that is critical to assess the quality and interpretation of the results is missing and primary\ndata is not present in the manuscript, supplementary figures or supplementary tables, which make it impossible to\nunderstand how conclusions are drawn. After the single cell analysis the authors use a genetic approach to try and validate\na few of the findings, focusing on biofilm-associated genes comP and bamE and iron-acquisition related genes bfrA and bfrB\nfrom mutagenesis and the RNAseq data, but this approach too has some weaknesses. Overall the manuscript presented\nhere used many techniques and boldly attempted to provide a comprehensive look at evolution of a pathobiont, but in its\ncurrent form it is missing too much pertinent information and contains several critical shortcomings that need to be\naddressed. There are many technical concerns which kept my revision to the most important major points, but there are also\nnumerous small mistakes in figures and text that would need to be addressed.\nMajor points:\n1) Section on lab-based evolution and phenotypic (biofilm) characterization:\na. Three replicates were used for the evolution, in keeping with standard practices in lab-based evolution to compare\nevolved trajectories in several parallel cultures but in most studies more than 3 are used (6 or 10 are current standards). Is\nthere a reason 3 was chosen? Also, the mutational signatures between replicates looks eerily similar in supplementary\nfigure 2J, are there any measures taken to make sure that no cross contamination occurred? Often times looking at repeated\nsynonymous mutations can identify whether this is an issue\nb. The authors refer to their evolution approach throughout the manuscript as “Directed Evolution” – a term that is used for\nmolecular evolution strategies used very differently than in the current manuscript to enhance activity of enzymes by\nfocusing on distinct parts of a functional gene for specific improvement. What they do to here is “lab-based evolution”. This\nmay sound like a semantic distinction but it is a very important difference as these are two highly distinct fields.\nc. There is critical information missing about the mutational analysis. The methods section does not actually discuss how\nlong reads were processed, what read coverage was and how similar mutations were called – these are all needed to\nunderstand what was done and the quality of the data and require a methods section – especially dealing with long read\nsequencing which (*using some platforms at least*) can have lower fidelity\nd. In initial characterization of biofilms all the data shown is from one replicate, it is important to see how different the\nreplicates are from each other.\n2) Section on traditional RNAseq transcriptomics:\na. The traditional RNAseq it appears triplicates were used (though this is not mentioned in the methods but there are 3\ncolumns in the figure for wildtype and the evolved strain). One question is which of the evolved strains was used and why all\nevolved strains were not tested?\nb. Gene enrichments are provided but no statistics are given and the list of genes used to define each category in the\ngeneset enrichment is not given. These are needed to determine what each gene set actually represents.\n3) Section on Single Cell RNAseq\na. A major flaw in logic from the first sections to single cell and bulk transcriptomics is that the authors hypothesize that\nbiofilm is important for antibiotic survival in the evolved strains, but both bulk and single cell transcriptomics were done with\nplanktonic cultures grown for only a few hours with high speed shaking, conditions that are not particularly suitable for\nbiofilms.\nb. The authors need to provide any statistics about the scRNAseq data quality (for example reads per cell, genes captured\nper cell, etc).\nc. The authors need to provide supplemental information with the differentially expressed genes in each cluster. For now\nthere is a heatmap but this does not provide pertinent information like P-values, log change, and genes beyond the selected\nfew genes per cluster.\nd. The authors make a claim that the difference in abundance in the evolved strain for the “metal binding cluster” are different\n– but in reality a change from 16.6% cells to 19.5% cells in a couple of replicates *(not clear how many if any replicates are\nused) does not seem statistically significant, and even if it is statistically significant biologically this would not likely be\nsignificant.\ne. Overall the UMAP of the wildtype and evolved strains look strikingly similar (as close to identical as I could imagine single\ncell data looking) – this makes it hard to claim that these are truly distinct.\nf. The authors ascribe biological function to clusters (“proliferating”, “iron scavenging” etc) – these cannot be made just\nbased on transcriptomic data, the authors need to validate these with physiological measurements. This is especially true for\nany definitions that rely on growth rates etc.\ng. Technique wise there is concern about the transcriptomic data because cells were processed for at least 10 minutes\nbefore fixation – bacterial mRNA half life is on the order of 1-2 minutes so many changes can be due to adaptation to PBS,\ncentrifugation, room temperature etc.\nh. The scRNAseq only compared wildtype to one of the 3 lineages in evolution, and not all 3, yet authors want to draw\nconclusions about general mechanisms of evolution. At the very least they need to change the language to explain that this\nmay not be a general result\ni. The pseudotime analysis relies entirely on a pipeline (Monocle2) that was developed for mammalian systems and often\nprovides erroneous temporal projections. Pseudotime analysis needs to be validated with timelapse or other kinetic\nmeasurements of gene expression that show that the proposed trajectories are actually meaningful. In bacterial transcription the\nactivation of genes is often pulsatile and there is no way to know if pulses are coordinated over time without dynamic\nmeasurements.\nj. The pseudotime comparison of the Wildtype and evolved lineage leave many questions about underlying data quality and\nthe analysis commands used, without seeing the actual pseudotime trajectory maps themselves it is hard for a reviewer to\njudge if the data is actually reasonable. In the limited data presented a troubling aspect is that there are duplicates of porB.\nUnfortunately none of the monocle steps or QC outputs (from Monocle or the other single cell analyses) are included in the\nsupplementary and the methods do not provide needed information to understand the data quality and exactly what was\ndone.\n4) Sections on Genetic Follow Up\na. The authors use plasmid-based fluorescent reporters, ideally chromosomal integration would be used to show\nheterogeneity and not plasmids.\nb. The differences seen in qPCR (figure 4e) appear to be statistically significant but it is not at all clear if such a small fold-\nwise effect (less than 2-fold) would actually have biological meaning.\nc. To show involvement of comP and bamE they were overexpressed with a strong promoter, it would be nice to see the\neffect of deletion or a weaker promoter.", "reviewer_2": "This study investigates ceftriaxone resistance in a Neisseria commensal, N. subflava, which can act as a pathobiont. The\nauthors evolve resistance by serial passage and antibiotic selection, and then characterize derived mutations and\ntranscriptional responses of the evolved lines in response to ceftriaxone. The authors find resiant strains almost always\nharbored ataA gene mutations. Resistant strains also upregulated genes involved in biofilm formation, comP and bamE. Iron\nsupplementation increased resistance further. The authors find heterogeneity among the transcriptional response of distinct\ncolonies. As not a lot of work has been done on antibiotic resistance in commensal Neisseria, I see this as an important\nstudy for the field. Commensals can be used as a model for the pathogenic Neisseria; and findings in this study may shed\nlight on undescribed resistance mechanisms in pathogens.\nSome minor comments:\n-Please add additional text in the discussion describing why there may be different lineages emerging with distinct\ntranscriptional responses within homogeneous cell cultures? Is this a function of where cells are in the biofilm? Or antigenic\nvariation? Why do the authors think this difference exists?\n-Do the authors find any of their evolved mutations in natural Neisseria populations? In vitro evolved cell lines may not be fit\nand survive in nature. Unless these mutations (or similar) are found in nature, the impact of these findings may be limited.\n-Throughout the text subflava is sometimes called “subflora”. Please correct."}, "author_response": "We sincerely thank the reviewers for their constructive comments. Below, we provide a point-\nby-point response to each of the concerns raised. All changes in the revised manuscript are\nhighlighted in blue font for clarity.\n----------\nReviewer #1 (Remarks to the Author):\n----------\nC1: The manuscript “Single-bacterial cell insights into mechanisms of ceftriaxone resistance\nin Neisseria subflava” presents a multifaceted study that aims to understand mechanisms of\nevolved resistance to antibiotics by Neisseria subflava. This bacterium has recently been\nimplicated with airway disease and presents a unique way to track genetic changes that allow\norganisms typically existing as commensals to become pathogenic. In this study the authors\nuse a lab-based evolution approach to evolve Neisseria in the presence of a relevant antibiotic\nin order to identify genetic changes that are accumulated. Triplicate cultures are used to try and\nuncover mutations that are repeated and are therefore deterministic. The study then checks\nwhether evolved strains have resistance to multiple antibiotics and attempts to characterize the\nmechanisms of resistance. This characterization is first done using physiological measurements\nthat demonstrate that the evolved cultures have greater biofilm producing capacity despite all\ncultures having mutation that inactivate a key biofilm regulator. Then the authors use traditional\nbulk RNAseq analysis to identify gene expression changes differentiating one of the lineages\nfrom the parental strain. After this analysis the authors use a bacterial single cell sequencing\napproach to try and understand the behavior of the culture. This single cell RNAseq\ncharacterization is the section of the manuscript that is used to make many of the less\nsubstantiated and more bold claims but unfortunately this section has several major issues with\nstudy design, analysis and interpretation. Importantly, the author's central claim is that the core\nmechanism of resistance is through biofilm formation, but the single cell analysis is all done\non planktonic cultures grown with much shaking. In addition, the authors point to varied states\nand trajectories that are supposed to be different between the parent strain and evolved culture\n(which apparently was only done for one of the lineages) - but the changes appear to be minimal\nand the trajectory analysis, which is notoriously a non-reliable single cell analysis method even\nin higher organisms that have defined developmental processes, is not validated by any\ndynamic measurement and makes several assumptions that are not founded. Furthermore, data\nthat is critical to assess the quality and interpretation of the results is missing and primary data\nis not present in the manuscript, supplementary figures or supplementary tables, which make\nit impossible to understand how conclusions are drawn. After the single cell analysis the authors\nuse a genetic approach to try and validate a few of the findings, focusing on biofilm-associated\ngenes comP and bamE and iron-acquisition related genes bfrA and bfrB from mutagenesis and\nthe RNAseq data, but this approach too has some weaknesses. Overall the manuscript presented\nhere used many techniques and boldly attempted to provide a comprehensive look at evolution\nof a pathobiont, but in its current form it is missing too much pertinent information and contains\nseveral critical shortcomings that need to be addressed. There are many technical concerns\nwhich kept my revision to the most important major points, but there are also numerous small\nmistakes in figures and text that would need to be addressed.\nR1: We sincerely thank Reviewer #1 for the comprehensive and insightful critique of our\nmanuscript. In response, we have made extensive revisions to address the major concerns raised,\nas summarized below: 1) We have expanded transcriptomic validation to include C26-b and\nC26-c strains, which confirmed the reproducibility of biofilm-related gene expression patterns\noriginally observed in C26-a. 2) Single-cell RNA-seq analyses were further refined, enabling\nthe identification of six transcriptionally distinct clusters. This was complemented by RT-\nqPCR validation of comP and bamE expression dynamics under both planktonic and sessile\nconditions. 3) We revised the pseudotime analysis and incorporated additional analyses using\nSlingshot and Palantir algorithms, revealing distinct evolutionary trajectories between wild-\ntype (WT) and resistant strains. 4) To dissect functional mechanisms, we fluorescently tagged\nand genetically manipulated comP and bamE to test their roles in resistance-associated\nphenotypes.\nWe are grateful for the reviewer’s insightful critique, which has guided substantial\nrevisions and additional validation experiments. We believe that the revised manuscript now\nmore rigorously and transparently presents the evolution and transcriptional adaptation of\nN. subflava under antibiotic selection.\nSection on lab-based evolution and phenotypic (biofilm) characterization:\nC2: Three replicates were used for the evolution, in keeping with standard practices in lab-\nbased evolution to compare evolved trajectories in several parallel cultures but in most studies\nmore than 3 are used (6 or 10 are current standards). Is there a reason 3 was chosen? Also, the\nmutational signatures between replicates looks eerily similar in supplementary figure 2J, are\nthere any measures taken to make sure that no cross contamination occurred? Often times\nlooking at repeated synonymous mutations can identify whether this is an issue\nR2: We sincerely thank the reviewer for raising these important points on replicate number,\nreproducibility, and cross-contamination safeguards. We agree and have addressed each point\ndirectly, as follows:\n1) In our study, we chose to derive three independent replicate lineages, following a\nstandard minimum benchmark used in experimental evolution to compare parallel evolutionary\ntrajectories. This design allowed us to perform deep multi-omic characterization, including\nwhole-genome sequencing, bulk RNA-seq, and single-cell RNA-seq (WT vs. C26-a), as well\nas detailed phenotypic profiling of biofilm and antibiotic resistance. Given the breadth and\nresolution of these downstream analyses, we opted for a focused number of replicates to ensure\nrigorous multi-omic comparisons across molecular and phenotypic dimensions. Despite using\nonly three replicates, we observed remarkably consistent evolutionary outcomes by generation\n26, including convergent phenotypes (e.g., increased biofilm formation and ceftriaxone\nresistance) and shared genomic alterations, most notably in the ataA gene, across all three\nlineages. Nonetheless, to further assess reproducibility at the transcriptomic level, we\nperformed additional bulk RNA-seq analyses on lineages C26-b and C26-c, in parallel with\nC26-a. As presented in Supplementary Fig. 3c-d, all three lineages showed highly similar\nexpression profiles, confirming that the observed transcriptional responses are reproducible\nand not specific to a single lineage. This confirms the consistency of the selection conditions\nand emergent adaptation in our model isolate SG001 by demonstrating a consistent adaptive\nresponse in 3 biological replicates (with 3 technical replicates) across genetic, phenotypic and\ntranscriptomic responses in response to ceftriaxone selection. We also note that we observed\nthe biofilm and resistance phenotypes in the intermediately evolved strain C9 (similarly) across\nthree independent replicates (C9-a, C9-b, C9-c) (Fig. 1d) affording additional reassurance. We\nhave now explicitly highlighted the number of replicates in our revised manuscript.\n2) We share the reviewer’s concern regarding cross-contamination, given the observed\nsimilarity in mutational signatures (Supplementary Fig. 3a-c). To ensure independence and\nintegrity of the replicates, we adopted stringent laboratory procedures: (i) Independent\nhandling of each lineage with separate consumables and workspace. (ii) Regular archiving of\nbacterial stocks every three generations, stored separately and clearly labeled. (iii) Validation\nthrough colony PCR before sequencing to verify lineage-specific mutations. We appreciate\nthe opportunity to clarify these important points.\n3) Following the reviewer’s suggestion, we examined synonymous SNP patterns across\nlineages and complemented this with SV-level analysis using long-read–informed assemblies.\nSNP-based clustering using Jaccard distances (Supplementary Fig. 3a), a method widely used\nin bacterial genomics, shows related but independent trajectories. Recent studies highlight that,\nunder strong selection, both synonymous and nonsynonymous mutations can arise in parallel\nacross independent populations, reflecting convergent adaptation rather than contamination1,2.\nStructural variants (SVs), which include large insertions, deletions, duplications, and other\ncomplex rearrangements, typically arise independently and are rarely subject to recurrent\nparallel mutation events observed with SNPs, making them highly reliable indicators for\nlineage differentiation3,4. Specifically, SV profiles are lineage-specific (C9-a/b/c) and separate\ncleanly by replicate in hierarchical analyses (Supplementary Data 2; Supplementary Fig.\n3b). Hierarchical clustering of SV data demonstrated unambiguously separate lineages, with\neach replicate displaying unique combinations of large-scale genomic alterations. Venn\ndiagrams in Supplementary Fig. 3c illustrate the overlap of mutations across the three lineages.\nMost of the shared variants were SNPs, with ataA (QFG95_10090) alone accounting for 63-\n77% of the total mutational load, highlighting ataA as a recurrent evolutionary hotspot under\nceftriaxone pressure, in contrast to structural variants that were more lineage-specific. This\nstrong lineage-specific resolution provided by SVs effectively rules out cross-contamination\nduring the experimental propagation process. Therefore, while we acknowledge the apparent\nsimilarity in SNP accumulation patterns, we argue that these similarities likely reflect parallel\nevolutionary pressures and mutation biases common in bacterial evolution, rather than cross-\ncontamination.\nIn sum, we (i) justify and make explicit the use of three independent lineages, (ii) add\nadditional bulk RNA-seq in C26-b/c to demonstrate reproducibility, and (iii) provide\northogonal genetic evidence (lineage-specific SVs) that rules out cross-contamination while\nacknowledging expected convergent SNPs under common selection. The manuscript and\nsupplement have been revised to surface these details with clear figure/panel and line\nreferences (lines 139-151). Thank you again for bringing this important issue to our attention\nand allowing us to enhance the clarity and robustness of our evolutionary analysis.\nC3: The authors refer to their evolution approach throughout the manuscript as “Directed\nEvolution” – a term that is used for molecular evolution strategies used very differently than in\nthe current manuscript to enhance activity of enzymes by focusing on distinct parts of a\nfunctional gene for specific improvement. What they do to here is “lab-based evolution”. This\nmay sound like a semantic distinction but it is a very important difference as these are two\nhighly distinct fields\nR3: We thank the reviewer for highlighting the use of the term “directed evolution”. To avoid\nconfusion with protein-engineering contexts, we have replaced it with “experimental\nevolution” throughout the manuscript. This term more accurately reflects our approach of\nserial passaging under antibiotic selection without targeted mutagenesis. Accordingly, we have\nrevised the manuscript to use the term “experimental evolution” consistently, which reflects\nthis semantic distinction and accurately conveys our experimental approach.\nC4: There is critical information missing about the mutational analysis. The methods section\ndoes not actually discuss how long reads were processed, what read coverage was and how\nsimilar mutations were called – these are all needed to understand what was done and the\nquality of the data and require a methods section – especially dealing with long read sequencing\nwhich (*using some platforms at least*) can have lower fidelity.\nR4: We appreciate the reviewer’s request for greater methodological transparency regarding\nlong- and short-read sequencing. We have revised the Methods section (lines 641-657) to\nexplicitly describe: (i) ONT base-calling with Guppy v6.4.6 (HAC) and subsequent QC (Q <\n7 filtering and adapter trimming); (ii) Illumina QC with FASTP v0.20.0 and explicit thresholds;\n(iii) read-coverage computation and breadth metrics using samtools depth (v1.16.1; ONT) and\nbedtools genomecov (v2.30.0; Illumina), with summaries in Supplementary Data 9 (all\nIllumina datasets covered ≥ 200× coverage across ≥ 95% of the genome, and ONT median\ndepths ranged 400–600× across samples); and (iv) inter-assembly variant calling by whole-\ngenome alignment with MUMmer v3.1/nucmer (parameters: --maxmatch -c 100 -b 200 -l 50)\nfor SNPs/indels and SyRI v1.4 (parameters: --unic 100 -k) for structural variants.\nC5: In initial characterization of biofilms all the data shown is from one replicate, it is important\nto see how different the replicates are from each other.\nR5: We thank the reviewer for emphasizing the importance of demonstrating replicate\nconsistency in biofilm experiments. Biofilm quantification has been expanded to include all\nC9-a/b/c and C26-a/b/c lineages, each performed with independent biological and technical\ntriplicates. Biofilm formation was quantified using the crystal violet staining assay in 96-well\nplates after static incubation, complemented by visual evaluation on Congo Red agar plates for\ncolony morphology and extracellular matrix production. Images from the Congo Red staining\nassays have been added to Fig. 1d, and this is described in the revised results section (lines\n120-123) and Methods (lines 538-544).\nSection on traditional RNAseq transcriptomics:\nC6: The traditional RNAseq it appears triplicates were used (though this is not mentioned in\nthe methods but there are 3 columns in the figure for wildtype and the evolved strain). One\nquestion is which of the evolved strains was used and why all evolved strains were not tested?\nR6: We appreciate the suggestion to clarify which lineages were profiled by bulk RNA-seq.\nOur initial bulk RNA-seq experiments presented in the original manuscript were performed\nusing three biological replicates each of the wild-type (WT) strain and the evolved strain\nlineage C26-a. Our rationale for initially focusing on C26-a was its representative phenotype,\nwhich included marked biofilm enhancement, ceftriaxone resistance, and clear genomic\nsignatures.\nIn response to the reviewer’s concern about reproducibility across independently evolved\nlineages, we performed additional bulk RNA-seq on C26-b and C26-c. As presented in\nSupplementary Fig. S3. d-e of the revised manuscript, these additional RNA-seq analyses\nconfirm highly consistent transcriptional profiles across all three evolved lineages (C26-a, C26-\nb, and C26-c). Specifically, the heatmaps of differentially expressed genes and relevant gene-\nsets demonstrate reproducible upregulation of pathways involved in biofilm formation and\nassociated metabolic and regulatory pathways. These results strongly reinforce our original\nfindings and provide robust support that the observed transcriptomic adaptations are\nreproducible across independently evolved replicate lineages. We have updated the revised\nmanuscript to clearly indicate these new data (lines 161-179).\nC7: Gene enrichments are provided but no statistics are given and the list of genes used to\ndefine each category in the geneset enrichment is not given. These are needed to determine\nwhat each gene set actually represents.\nR7: We thank the reviewer for requesting more transparent enrichment statistics. We have now\nexplicitly provided detailed statistical data (including p-values, adjusted p-values, and\nlog (Fold change)) for the gene sets highlighted in our manuscript, particularly those associated\n2\nwith the “cell adhesion” and “two-component regulatory system” pathways. This detailed\ninformation is presented in the revised Supplementary Data 3 and Supplementary Data 4.\nWe sincerely thank the reviewer for prompting us to clarify this crucial aspect of our data\npresentation, thus significantly enhancing the clarity, reproducibility, and rigor of our analyses.\nSection on Single Cell RNAseq:\nC8: A major flaw in logic from the first sections to single cell and bulk transcriptomics is that\nthe authors hypothesize that biofilm is important for antibiotic survival in the evolved strains,\nbut both bulk and single cell transcriptomics were done with planktonic cultures grown for\nonly a few hours with high speed shaking, conditions that are not particularly suitable for\nbiofilms.\nR8: We sincerely thank the reviewer for this important point. While biofilm-associated survival\nis central to our hypothesis, technical limitations preclude reliable single-cell profiling of\nbiofilm-derived bacteria, as extracellular polymeric substances impair cell capture and RNA\nrecovery. To ensure data quality, we analyzed bacteria uniformly under planktonic conditions:\n1) Technical limitations of single-cell transcriptomics: The primary limitation is technical;\nsingle-cell RNA sequencing platforms currently have limited compatibility with biofilm-\nderived samples. Specifically, the abundant extracellular polymeric substances (EPS) produced\nin mature biofilms severely hinder efficient single-cell isolation, encapsulation, and processing,\nleading to poor cell capture efficiency and compromised RNA quality. To ensure high-quality,\nreproducible, and biologically meaningful single-cell transcriptomic data, we therefore chose\nto process cells uniformly in planktonic conditions.\n2) Ensuring comparability and reproducibility: All transcriptomic analyses (bulk RNA-\nseq and single-cell RNA-seq of wild-type and C26 strains) were conducted under rigorously\ncontrolled, identical planktonic conditions. This approach allowed us to accurately identify\nintrinsic transcriptional signatures associated specifically with the evolved strain, capturing\nearly molecular responses indicative of a “biofilm-adapted” transcriptional state.\n3) Rationale for biofilm-related gene expression in planktonic conditions: Accumulating\nevidence from recent studies indicates that biofilm-associated genes can indeed be\ntranscriptionally primed or induced even in planktonic cultures, particularly under\nenvironmental stress or antibiotic pressure. For instance, planktonic aggregates of\nStaphylococcus epidermidis exposed to environmental stress showed characteristic biofilm\ngene signatures5. Likewise, sub-inhibitory antibiotic stress in planktonic cells triggered\nbiofilm-like transcriptional responses, including enhanced expression of EPS and stress-\nresponse regulators6. Similar phenomena were observed in planktonic organisms, where stress\ninduced biofilm-associated gene expression patterns7, consistent with selective pressures\nencountered during our antibiotic-driven evolution experiments.\n4) Experimental validation and supporting evidence: To explicitly address the reviewer’s\ncomment, we performed further RT-qPCR validation of two representative biofilm-associated\ngenes (comP and bamE) under both planktonic (shaking) and sessile (biofilm-promoting)\nconditions, comparing expression levels in evolved (C26) and wild-type strains\n(Supplementary Fig. 4g). These measurements showed that in C26, comP and bamE were\ncomparably upregulated under both planktonic and sessile conditions, whereas wild-type cells\nshowed no mode-dependent change. This indicates an intrinsic transcriptional predisposition\ntoward biofilm formation in the evolved lineage. We have now explicitly included and clarified\nthese important experimental validations together with explanatory details in our revised\nmanuscript (lines 184-191 and 279-298).\nC9: The authors need to provide any statistics about the scRNAseq data quality (for example\nreads per cell, genes captured per cell, etc).\nR9: We thank the reviewer for noting the need for quantitative QC metrics. We now provide\nfull scRNA-seq quality summaries, including total cells, median reads per cell, and detected\ngenes per cell, in Supplementary Fig. 4a. These demonstrate comparable data quality between\nWT and C26 samples. in the revised main text (lines 193-197).\nC10: The authors need to provide supplemental information with the differentially expressed\ngenes in each cluster. For now there is a heatmap but this does not provide pertinent information\nlike P-values, log change, and genes beyond the selected few genes per cluster.\nR10: We thank the reviewer for this important suggestion. We have substantially expanded our\nsupplementary materials to ensure transparency and allow deeper evaluation of our single-cell\ndata. We now provide complete cluster-specific DEG tables listing gene identifiers, log₂ fold-\nchange, p-value, adjusted p-value, and detection rate. All data are accessible in Supplementary\nData 7, and filtering criteria are detailed in the Methods section (line 599-610).\nC11: The authors make a claim that the difference in abundance in the evolved strain for the\n“metal binding cluster” are different – but in reality a change from 16.6% cells to 19.5% cells\nin a couple of replicates *(not clear how many if any replicates are used) does not seem\nstatistically significant, and even if it is statistically significant biologically this would not\nlikely be significant.\nR11: We appreciate the reviewer’s concern that cluster frequency changes appeared modest.\nWe have also now reanalyzed the data to assess whether the observed differences were\ninfluenced by the Canonical Correlation Analysis (CCA) integration step, which can reduce\nsignal strength in prokaryotic datasets due to their lower transcript counts compared to\neukaryotes. Specifically, we compared the log FC of genes expressed in at least 5% of cells\n2\nfrom both WT and C26-resistant strains, with and without CCA integration. While the\nproportion of statistically significant genes remained similar between the two approaches\n(Figure R1), we noted that the expression levels of several essential genes, including ribosomal\ngenes (rpoB, rplE, rpsA, rpmL), shikimate pathway genes (aroA), and DNA transcription-\nrelated genes (gyrB), were attenuated in the CCA-integrated dataset (Figure R2). These finding\nsuggests that CCA, while useful in reducing noise, may also obscure true biological differences\nin bacterial datasets. We appreciate the reviewer’s suggestion, which has helped us to clarify\nthis limitation and provide a more balanced interpretation of our analysis.\nFigure R1: Log (FC) of C26 vs WT strains with and without CCA integration step.\nIn view of this issue, we re-analyzed the single-cell RNAseq data without CCA integration.\nThe new analysis showed that WT and C26 resistant strains formed distinct clustering patterns,\ncomprising six populations with unique transcriptomic signatures. Cluster 5, which exhibited\nmetal-binding properties, was more abundant in the resistant strains (16.3%) compared to WT\nstrains (0.5%).\nTo further evaluate the functional significance of our findings in relation to the metal-\ndependent cluster, we performed ICP-MS analysis to directly measure intracellular metal ion\nconcentrations. Notably, the evolved C26 strain contained significantly elevated intracellular\niron levels relative to the wild type (Supplementary Fig. 4d). Combined with antibiotic\nsensitivity testing, this confirms an iron-dependent resistance phenotype, consistent with the\ntranscriptomic signature of a metal-binding cluster. The specific details are described in lines\n266-269 of the revised manuscript. Metal quantification protocols are described in detail in the\nMethods section (lines 696-709). Thus, although the transcriptional shifts appeared modest,\ncomplementary physiological assays demonstrate that iron homeostasis is indeed a key\nadaptive axis in the resistant strain.\nC12: Overall the UMAP of the wildtype and evolved strains look strikingly similar (as close to\nidentical as I could imagine single cell data looking) – this makes it hard to claim that these are\ntruly distinct.\nR12: We reanalyzed our scRNA-seq data without the CCA integration step, which we now\nrecognize may have obscured subtle but biologically meaningful differences (Fig. 2b-c). The\nupdated UMAP visualization resolved six transcriptomically distinct populations, differences\nin distribution that are now clearly visible upon omission of the CCA analysis: WT bacteria\nwere predominantly found in clusters 1 (11.5%), 2 (48.1%), and 3 (34.6%), while resistant\nbacteria primarily occupied clusters 1 (34.7%), 4 (34.4%), 5 (16.3%), and 6 (7.9%). Notably,\ncluster 1, the major shared population, was characterized by outer membrane porin (porB) and\nribosomal genes (rplC, rpsC), suggesting conserved basal functions. However, clusters 2 and\n3 were enriched in WT bacteria, whereas clusters 4-6 were enriched in resistant bacteria. These\nfindings confirm that the strains are transcriptomically distinct, albeit with partial overlap in\ncluster 1.\nC13: The authors ascribe biological function to clusters (“proliferating”, “iron scavenging” etc)\n– these cannot be made just based on transcriptomic data, the authors need to validate these\nwith physiological measurements. This is especially true for any definitions that rely on growth\nrates etc.\nR13: We thank the reviewer for this important point. We agree that functional assignments\nbased solely on transcriptional profiles are predictive and require validation. Accordingly, we\nnow provide population-level validations aligned to the transcript-defined trends and clearly\nindicate their locations: (i) growth kinetics show accelerated growth and higher stationary-\nphase density in C26, consistent with a proliferative program (Supplementary Fig. 4e; lines\n269-272; Methods, lines 545-556); (ii) a bfrA-gfp genomic reporter demonstrated marked\nupregulation of bfrA expression upon iron supplementation (Supplementary Fig. 4f; lines\n272-278; Methods, lines 785-794), supporting the assignment of a “metal-binding” cluster;\nand (iii) ICP-MS confirms elevated intracellular iron in C26, providing an orthogonal readout.\nNevertheless, we recognize that direct, cluster-resolved functional validation remains\ntechnically challenging in bacteria due to the absence of unique surface markers and the\ncomplexity of bacterial gene regulatory networks. As a result, our functional assignments\nremain inferentially supported by population-level validation rather than definitive single-\ncluster assays. We have discussed these methodological limitations and the rationale for our\napproach in the revised manuscript (lines 474-495). We thank the reviewer for comments that\nstrengthened the rigor and transparency of the study.\nC14: Technique wise there is concern about the transcriptomic data because cells were\nprocessed for at least 10 minutes before fixation – bacterial mRNA half life is on the order of\n1-2 minutes so many changes can be due to adaptation to PBS, centrifugation, room\ntemperature etc.\nR14: We thank the reviewer for this important point. We agree that short bacterial mRNA half-\nlives can make handling steps consequential. We took several measures to minimize such\neffects. We processed all samples on ice (4 °C) with pre-chilled, isotonic buffers, and kept\nharvest-to-fixation <10 minutes; WT and C26 were handled in parallel under identical\nconditions to preserve comparative integrity. We now describe these safeguards in Methods\n(“Preprocessing for single-cell sequencing”) (Methods lines 557–576). While brief handling\nmay still influence absolute levels, the relative contrasts we report are robust and supported by\nboth bulk and single-cell data.\nC15: The scRNAseq only compared wildtype to one of the 3 lineages in evolution, and not all\n3, yet authors want to draw conclusions about general mechanisms of evolution. At the very\nleast they need to change the language to explain that this may not be a general result.\nR15: We appreciate this comment and fully acknowledge this limitation. Our single-cell RNA-\nseq analysis was performed only on the C26-a lineage, selected as a representative evolved\npopulation despite the extensive transcriptional diversity observed across all lineages in bulk\nRNA-seq. At the same time, bulk RNA-seq across C26-a/b/c shows highly concordant\npathway-level changes (biofilm, iron, envelope programs), indicating that the major adaptive\nthemes are reproducible across independently evolved lineages.\nAccordingly, all single-cell findings, such as the nature of cellular heterogeneity, cell-state\ntransitions, and adaptive gene expression programs, should be interpreted within the specific\ncontext of C26-a, and may not necessarily reflect the full spectrum of evolutionary responses\npresent in other evolved lineages. We have revised the manuscript throughout to clarify this\nscope, explicitly stating that mechanistic insights at the single-cell level pertain to C26-a alone.\nNevertheless, the consistent bulk RNA-seq signatures across all lineages provide strong\nsupport for shared major adaptive themes, and our approach illustrates the value of integrating\nboth bulk and single-cell perspectives. We agree that further single-cell analyses of additional\nevolved lineages will be required to fully assess the diversity and reproducibility of these\nadaptive mechanisms. This limitation is now clearly stated in the revised Discussion section\n(line 474-495).\nC16: The pseudotime analysis relies entirely on a pipeline (Monocle2) that was developed for\nmammalian systems and often provides erroneous temporal projections. Pseudotime analysis\nneeds to be validated with timelapse or other kinetic measurements of gene expression that\nshow that the proposed trajectories are actually meaningful. In bacterial transcription the\nactivation of genes is often pulsatile and there is no way to know if pulses are coordinated over\ntime without dynamic measurements.\nR16: We thank the reviewer for this thoughtful comment. To address concerns, we reanalyzed\nthe single-cell data using two distinct algorithms, Slingshot and Palantir, both of which have\nbeen successfully applied to bacterial and other non-mammalian contexts8,9. Palantir’s utility\nin bacterial single-cell RNA-seq analyses has been validated in previous studies, such as\nmapping developmental stages in Staphylococcus aureus. We independently applied the\nPalantir algorithm, which uses a diffusion map-based approach distinct from Slingshot’s\nminimum spanning tree method. Both algorithms consistently revealed distinct evolutionary\ntrajectories between WT and resistant strains. The WT population showed a single linear\ntrajectory ending at cluster 3, representing a homeostatic state. In contrast, resistant strains\nexhibited a branched trajectory with two distinct endpoints: cluster 5 (metal sequestration) and\ncluster 6 (sulfur utilization) (Supplementary Fig. 4b-c). This branching pattern suggests\ncephalosporin pressure drives multiple, stable resistance strategies through divergent genetic\nadaptations (revised in lines 217-226).\nC17: The pseudotime comparison of the Wildtype and evolved lineage leave many questions\nabout underlying data quality and the analysis commands used, without seeing the actual\npseudotime trajectory maps themselves it is hard for a reviewer to judge if the data is actually\nreasonable. In the limited data presented a troubling aspect is that there are duplicates of porB.\nUnfortunately none of the monocle steps or QC outputs (from Monocle or the other single cell\nanalyses) are included in the supplementary and the methods do not provide needed\ninformation to understand the data quality and exactly what was done.\nR17: We thank the reviewer for raising important points regarding our pseudotime analysis and\nthe underlying data quality. We have now substantially enhanced the rigor and transparency of\nthis analysis, as detailed below:\n1) Trajectory visualization and cross-method validation: We have now included detailed\nvisualizations of all inferred trajectories in Supplementary Fig. 4b and Supplementary Data\n6, providing comprehensive visualizations of the inferred pseudotime trajectories for each\nbacterium as determined by the Slingshot and Palantir algorithms. As shown in Supplementary\nFig. 4a, violin plots of reads per cell and genes per cell demonstrate comparable data quality\nbetween WT and C26 populations. Low-quality cells were excluded prior to analysis, and these\nQC procedures are now explicitly described in the revised Methods section (lines 611-619).\nWe used Slingshot and Palantir algorithms (not Monocle) to infer lineage trajectories. The\nresulting pseudotime maps are presented in Supplementary Fig. 4b-c, which consistently\nreveal distinct evolutionary trajectories between WT and resistant populations. The WT cells\nfollowed a single linear trajectory, culminating in cluster 3 (homeostatic state), whereas the\nresistant strains exhibited a branched trajectory with endpoints in clusters 5 (metal\nsequestration) and 6 (sulfur utilization).\n2) PorB duplication correction and data transparency: We corrected duplicate porB entries\nand updated the gene expression matrix. All computational steps, QC outputs, and key data\ntables, including pseudotime, branch assignments, and gene expression, are now fully provided\nin the Supplementary Materials (Supplementary Data 6). We trust these comprehensive\nadditions and clarifications fully address the reviewer’s concerns regarding data quality,\nanalytical robustness, and the reproducibility of our pseudotime analysis.\nSections on Genetic Follow Up:\nC18: The authors use plasmid-based fluorescent reporters, ideally chromosomal integration\nwould be used to show heterogeneity and not plasmids.\nR18：We sincerely appreciate the reviewer’s valuable suggestions. We have replaced the\noriginal plasmid-based fluorescent reporter system with a chromosomally integrated construct\ngenerated by homologous recombination. Specifically, a fluorescent marker gene (mCherry or\nGFP) and a gentamicin resistance cassette were inserted in-frame at the native comP or bamE\nlocus, under control of the endogenous promoter, using approximately 700 bp homologous\narms. This strategy ensures that fluorescence reflects intrinsic, physiologically relevant gene\nexpression dynamics, eliminating variability due to plasmid copy number or expression\nartifacts. To further address the reviewer's concern, we have revised Fig. 4a with new data\nobtained from these chromosomally labeled strains. Integration was verified by PCR and\nSanger sequencing, and correct expression was confirmed by fluorescence microscopy. This\nrevised approach allows us to directly visualize and quantify single-cell expression\nheterogeneity of comP and bamE in both wild-type and resistant backgrounds, as now detailed\nin the revised manuscript (Fig. 4a, lines 312-319). The updated methodology is also described\nin detail in the revised Methods section (lines 767-784). We appreciate the reviewer’s\nsuggestion, which has substantially improved the technical rigor and interpretability of our\nresults.\nC19: The differences seen in qPCR (figure 4e) appear to be statistically significant but it is not\nat all clear if such a small fold-wise effect (less than 2-fold) would actually have biological\nmeaning.\nR19: We thank the reviewer for raising this important point. As comP and bamE are\nheterogeneity-associated genes identified at the single-cell level, their upregulation occurs\npredominantly within specific bacterial clusters. Bulk RNA-seq or qPCR analyses dilute strong\nsignals from subsets of cells, underestimating their impact, even when these genes play decisive\nroles within critical functional subsets. To further clarify this issue, we repeated the qPCR\nexperiments with additional biological replicates, and the results updated in Fig. 4e consistently\nreproduced the modest but statistically significant upregulation. Importantly, to confirm\nbiological relevance, we conducted quantitative biofilm formation assays, which demonstrated\nthat even modest transcriptomic shifts in comP and bamE are associated with significant\ndifferences in resistance and biofilm phenotypes (Fig. 4g). Thus, in the context of single-cell\nheterogeneity, even modest changes in gene expression may nonetheless underpin biologically\nrelevant phenotypic adaptations, as supported by our functional validation.\nC20: To show involvement of comP and bamE they were overexpressed with a strong promoter,\nit would be nice to see the effect of deletion or a weaker promoter.\nR20: We sincerely thank the reviewer for suggesting alternative strategies to strengthen our\nfunctional validation. We have employed a chromosomal integration approach through precise\ngene deletions of comP and bamE in N. subflava. We used homologous recombination to\ngenerate knockout constructs with ~700-bp homology arms flanking each target gene. Each\nknockout cassette contained a selectable marker (gentamicin resistance). Electroporation was\nused to introduce linearized knockout constructs into electrocompetent Neisseria cells,\nfollowed by selection and confirmation by PCR and Sanger sequencing. We successfully\ngenerated ΔcomP and ΔbamE strains. Using these deletion mutants, we conducted phenotypic\nassays to delineate their contributions to ceftriaxone resistance and biofilm formation under\nbiologically relevant conditions. These comprehensive analyses are now clearly presented in\nthe revised manuscript (lines 322-326) and Methods (lines 746-766), demonstrating a rigorous\napproach toward validating the biological significance of these genes.\nReviewer #2 (Remarks to the Author):\n----------\nC21: This study investigates ceftriaxone resistance in a Neisseria commensal, N. subflava,\nwhich can act as a pathobiont. The authors evolve resistance by serial passage and antibiotic\nselection, and then characterize derived mutations and transcriptional responses of the evolved\nlines in response to ceftriaxone. The authors find resistant strains almost always harbored ataA\ngene mutations. Resistant strains also upregulated genes involved in biofilm formation, comP\nand bamE. Iron supplementation increased resistance further. The authors find heterogeneity\namong the transcriptional response of distinct colonies. As not a lot of work has been done on\nantibiotic resistance in commensal Neisseria, I see this as an important study for the field.\nCommensals can be used as a model for the pathogenic Neisseria; and findings in this study\nmay shed light on undescribed resistance mechanisms in pathogens.\nR21: We sincerely thank this reviewer for their kind appraisal and supportive comments.\nC22: Please add additional text in the discussion describing why there may be different lineages\nemerging with distinct transcriptional responses within homogeneous cell cultures? Is this a\nfunction of where cells are in the biofilm? Or antigenic variation? Why do the authors think\nthis difference exists?\nR22: We respectfully thank the reviewer for these incisive questions regarding the origins and\nsignificance of transcriptional heterogeneity and lineage diversification in our evolved\nNeisseria subflava populations, as well as the clinical relevance of key adaptive mutations. In\nthe revised Discussion, we have now explicitly integrated a mechanistic and ecological\nframework to clarify these points:\nWe clarify emphasize that the pronounced transcriptional heterogeneity and emergence\nof distinct lineages observed in our evolved populations are not experimental artifacts, but\nreflect fundamental adaptive strategies. Despite rigorous efforts to standardize experimental\nconditions, we consistently observed lineage diversification within N. subflava populations.\nThis phenomenon most likely results from the interplay of subtle spatial microgradients (in\nnutrients, oxygen, or antibiotic exposure) and intrinsic regulatory processes, such as stochastic\ngene expression, phase variation, and the action of small RNAs and transcription factors.\nDespite planktonic processing, transient micro-aggregates and surface-adherent\nsubpopulations likely create local biofilm-like niches. Enrichment of comP (competence/pilus\nassembly) and bamE (BAM-complex assembly) in resistant clusters points to micro-niche\nformation and outer-membrane fortification under β-lactam pressure, consistent with β-lactam–\nprimed biofilm programs rather than reliance on pre-formed biofilms. Adaptive differentiation\nmay arise even in macroscopically homogeneous cultures due to localized fluctuations in\nenvironmental conditions, most notably within micro-niches that resemble biofilm-like\nenvironments. While antigenic (pilin) variation is a well-documented source of phenotypic\ndiversity in Neisseria10, our results suggest that the observed transcriptional heterogeneity and\nlineage bifurcation primarily stem from adaptive regulatory mechanisms rather than large-scale\nantigenic switching (Discussion, lines 423-446).\nC23: Do the authors find any of their evolved mutations in natural Neisseria populations? In\nvitro evolved cell lines may not be fit and survive in nature. Unless these mutations (or similar)\nare found in nature, the impact of these findings may be limited.\nR23: We agree with the reviewer that that translational relevance depends on whether mutations\nalso appear in natural populations. We found that key adaptive mutations (ataA, comP, and\nbamE) were present in clinical isolates from bronchiectasis patients, especially in ST-a lineages.\nClinical strains, especially those of the ST-a lineage, displayed upregulation of comP and bamE,\nincreased biofilm formation, and reduced antibiotic susceptibility, closely mirroring our in vitro\nfindings. This convergence reinforces the clinical and translational relevance of our work.\nHowever, we also recognize that some adaptations, such as mutations in iron metabolism genes\n(bfrA, bfrB), are much less prevalent in natural populations, likely due to differing selective\npressures, ecological context, and population history11,12. Notably, the mutational signatures of\nthe three independently evolved C9 ceftriaxone-resistant strains exhibited striking similarities.\nThe C9 strains accumulated 190, 203, and 166 mutations, respectively, with 149 shared across\nall three lineages (Supplementary Fig. 3c). The majority of shared SNPs clustered in the\nadhesion-associated ataA gene (QFG95_10090), underscoring its role as a hotspot of\nconvergent evolution under ceftriaxone pressure. In contrast, the three C9 strains displayed\ndistinct structural variant (SV) signatures, including large insertions, deletions, duplications,\nand complex rearrangements. These strain-specific SVs occurred in genes beyond ataA. For\nexample, in C9-a, multiple deletions were identified at QFG95_07970 (mltA). C9-b exhibited\nlineage-specific insertions and deletions affecting QFG95_01255 (CRISPR-associated\nendonuclease Cas3), while C9-c had insertions and deletions at QFG95_05805\n(tetratricopeptide repeat protein). Unique mutational signatures of the three independently\nevolved N. subflava strains strongly suggest additional divergent genomic adaptations in\nresponse to antibiotic stress.\nOverall, our study demonstrates that transcriptional heterogeneity and lineage\ndiversification underpin robust, multi-layered adaptive responses to antibiotic pressure in N.\nsubflava. The striking convergence of laboratory and clinical evolutionary trajectories\nreinforces the ecological and clinical validity of the adaptive strategies identified here, while\nalso underscoring the importance of expanded environmental sampling to fully define the\ndistribution and impact of these resistance mechanisms.\nAll these points are now clearly and comprehensively discussed in the revised manuscript,\nwith detailed mechanistic rationale and relevant literature support (Discussion, lines 447–459).\nC24: Throughout the text subflava is sometimes called “subflora”. Please correct.\nR24: We have carefully corrected all instances of “subflora” to the correct species name\n“subflava” throughout the manuscript.\nReferences\n1. Kinnersley, M., Schwartz, K., Yang, D.-D., Sherlock, G. & Rosenzweig, F. Evolutionary dynamics and\nstructural consequences of de novo beneficial mutations and mutant lineages arising in a constant\nenvironment. BMC biology 19, 20; 10.1186/s12915-021-00954-0 (2021).\n2. Sandberg, T. E., Szubin, R., Phaneuf, P. V. & Palsson, B. O. Synthetic cross-phyla gene replacement and\nevolutionary assimilation of major enzymes. Nature ecology & evolution 4, 1402–1409; 10.1038/s41559-\n020-1271-x (2020).\n3. Couce, A. et al. Changing fitness effects of mutations through long-term bacterial evolution. Science (New\nYork, N.Y.) 383, eadd1417; 10.1126/science.add1417 (2024).\n4. Schenk, M. F. et al. Population size mediates the contribution of high-rate and large-benefit mutations to\nparallel evolution. Nature ecology & evolution 6, 439–447; 10.1038/s41559-022-01669-3 (2022).\n5. Bottagisio, M. et al. Proteomic Analysis Reveals a Biofilm-Like Behavior of Planktonic Aggregates of\nStaphylococcus epidermidis Grown Under Environmental Pressure/Stress. Front. Microbiol. 10, 1909;\n10.3389/fmicb.2019.01909 (2019).\n6. Jia, J., Xue, X., Guan, Y., Fan, X. & Wang, Z. Biofilm characteristics and transcriptomic profiling of\nAcinetobacter johnsonii defines signatures for planktonic and biofilm cells. Environmental Research 213,\n113714; 10.1016/j.envres.2022.113714 (2022).\n7. Yao, S. et al. Multispecies biofilms in fermentation: Biofilm formation, microbial interactions, and\ncommunication. Comprehensive reviews in food science and food safety 21, 3346–3375; 10.1111/1541-\n4337.12991 (2022).\n8. Setty, M. et al. Characterization of cell fate probabilities in single-cell data with Palantir. Nature\nbiotechnology 37, 451–460; 10.1038/s41587-019-0068-4 (2019).\n9. Korshoj, L. E. & Kielian, T. Bacterial single-cell RNA sequencing captures biofilm transcriptional\nheterogeneity and differential responses to immune pressure. Nature communications 15, 10184;\n10.1038/s41467-024-54581-8 (2024).\n10. Wielert, I. et al. Pilin antigenic variants impact gonococcal lifestyle and antibiotic tolerance by modulating\ninterbacterial forces. PLoS biology 23, e3003022; 10.1371/journal.pbio.3003022 (2025).\n11. Baquero, F. et al. Evolutionary Pathways and Trajectories in Antibiotic Resistance. Clinical microbiology\nreviews 34, e0005019; 10.1128/CMR.00050-19 (2021).\n12. Jangir, P. K., Ogunlana, L. & MacLean, R. C. Evolutionary constraints on the acquisition of antimicrobial\npeptide resistance in bacterial pathogens. Trends in microbiology 29, 1058–1061;\n10.1016/j.tim.2021.03.007 (2021)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The revisions by Zhang et al for their manuscript “Single-bacterial cell insights into mechanisms of ceftriaxone resistance in Neisseria subflava” nicely address the majority of my concerns but there are a few remaining items (motly minor) that could benefit from further clarification. First – there are a few instances where the authors use imprecise language or language that is not directly supported by the results. This does not take away from the general quality of the paper but it is important for the authors to be correct in their use of language and interpretation, to avoid adding confusion to the field. Second, the pseudotime analysis is still a concern. I am not sure what this analysis adds to the paper but there are reasons to question the interpretation of pseudotime in bacterial cultures. I add specifics to this below, and I generally appreciate that the authors included some cautionary words in their “limitations of the study” section that they added in this revision, but unless there is a real reason to include pseudotime I would urge the authors to drop this section and focus on the other data that they present which is clear and compelling.\nInstances of imprecise language or imprecise interpretation\n1) Line 224: The authors claim that they show “Multiple specialized resistance strategies” and that “two terminal clusters (5 and 6) represent evolutionary stable solutions to antibiotic stress.” – This statement is not supported by the results. The authors do not show which cells survive challenge and it is possible that cells in all clusters survive challenge. Additionally, since the authors do not follow dynamics of reporters, it is not clear that cells stay in these clusters for any length of time. These statements would therefore necessitate significant additional experiments if the authors wish to make them. Alternatively, the authors can remove suggestive language and focus on the other interesting results for which they provide direct evidence.\n2) Line 318: The authors claim that comP and bamE expression has overall higher intensity and heterogeneity in the evolved strain. They are absolutely correct about the intensity, but I am not convinced that the heterogeneity increased. In figure 4a we see bimodal expression in the wildtype strain in both comP and bamE. In fact, one may think that at least for bamE the WT sample appears to the naked eye more bimodal than C26 (at least you would say more cells are in an OFF/ON ratio, whereas almost all cells in the C26 culture are “ON”) – but the point is that if we are making claims about the extent of heterogeneity there must be statistical measurements and distributions need to be reported. I think that the important claim here is not about the level of heterogeneity but the overall increase in expression, but if the authors want to keep the claim on increased heterogeneity they must show distribution of single cell fluorescence. Another caveat here is that the reporter used is a stable fluorescent protein (GFP) and therefore may accumulate differently as cell growth dynamics play a role in overall intensity of stable fluorescent proteins.\n3) In line 218 the authors claim that pseudotime reveals distinct evolutionary trajectories between WT and the resistant strain. The pseudotime is of gene expression, not DNA and therefore does not measure evolutionary trajectories. In addition, please see comments below about pseudotime.\nConcerns about the use of pseudotime:\nThe addition of pseudotime analysis in the context of an evolutionary experiment raises many concerns. First, in an evolutionary experiment there are expected to be numerous coexisting phylotypes that are genetically distinct, and each can occupy distinct transcriptional landscapes. Since we do not know if these trajectories are linked temporally at the DNA level (this study did not study DNA lineages of single cells alongside transcription) it is impossible to determine what parts of the transcriptional landscape is due to temporal changes in gene expression. In addition, the term “pseudotime” is highly misleading, this is a trajectory map that aligns cells based on some sort of path alongside a continuum – or gradient. In the context of biofilm cells there is a large likelihood that the dominant gradient will be spatial (like oxygen gradients or nutrient gradients) – and not temporal. The authors attempt to link the transcriptional changes to a temporal order but the main genes are involved in biofilm production and even though the cells were grown with shaking there is a chance that aggregations produce temporal gradients. Additionally, it is not clear whether in bacterial cells pseudotime provides meaningful results in processes that are not developmental (in other words, processes that are not stepwise like sporulation in bacteria are often states in which cells switch on and then off, but not a temporal developmental process like in mammalian cells). It is true that other papers have used a highly questionable pseudotime approach to look at microbiome samples (in which much of the heterogeneity is likely due to distinct strains that co-exist in the gut) – and to look at a dataset in Staph aureus but there are several significant concerns about that dataset. In short, without validation using fluorescent reporters the authors pseudotime claims will only add confusion to the field and these results do not add any significant interpretations to the study so my suggestion to the authors is to either doing the careful dynamic experiments or dropping the pseudotime analysis altogether."}, "author_response": ""}]

Decision letter