The wild-type G. parasuis strain MY1902 (serotype 13) was generously provided by the Swine Disease Research Center at Sichuan Agricultural University. The ΔqseC mutant and its complemented strain (c-ΔqseC) were constructed via homologous recombination as detailed in previous study (He et al., 2018), with all strains cryopreserved at −80 °C in our laboratory repository. The ΔqseC mutant was constructed by replacing the gene with a kanamycin resistance cassette using overlap PCR and homologous recombination via the suicide vector pK18mobsacB. The complementation strain was generated by cloning the qseC gene into the replicative vector pSF116 and introducing it into the mutant. Strains were cultured Tryptic Soy Broth/Agar (TSB/TSA; BD-Difco) supplemented with 5% heat-inactivated newborn bovine serum (Solarbio) and 0.1% (w/v) NAD (Sigma-Aldrich). Kanamycin (100 μg/mL) was added for mutant selection.

The effect of qseC gene deletion on the morphology of G. parasuis was assessed using TEM. The wild-type strain MY1902 and ΔqseC mutant strain were inoculated into fresh TSB at a 1:100 dilution and cultured at 80 rpm for 24 h. When the OD_{600 nm} reached approximately 1.2, the culture was centrifuged at 4,000 rpm for 10 min. The resulting bacterial pellet was gently rinsed once with phosphate-buffered saline (PBS) and re-centrifuged at 5,000 rpm for 2 min. After discarding the PBS supernatant, we added 3 mL of electron microscopy fixative (2.5% glutaraldehyde solution) to the pellet. The pellet was then gently resuspended and dispersed in the fixative. Fixation was performed at room temperature in the dark for 2 h, followed by storage at 4 °C. The fixed samples were subsequently processed and imaged using a transmission electron microscope by a commercial service provider. Images were acquired and analyzed to evaluate morphological differences.

Intracellular metabolites were extracted from exponentially growing cultures (OD600 ≈ 0.8) of both the G. parasuis MY1902 and the ΔqseC mutant strain to compare their metabolite profiles. Cells were harvested by centrifugation (4 °C, 4,000 rpm, 5 min) from TSB cultures. A representative aliquot of the resulting bacterial pellet was subjected to metabolite extraction. For this purpose, the pellet aliquot was resuspended in 300 μL of ice-cold 80% (v/v) aqueous methanol to quench metabolism and initiate extraction. Immediate flash-freezing in liquid nitrogen for 5 min was performed to further arrest metabolic activity and facilitate cell disruption. The frozen suspension was then thawed on ice and vortexed vigorously for 30 s. Complete cell lysis was achieved by ultrasonication for 6 min on ice to ensure efficient metabolite release. Cellular debris was removed by centrifugation (4 °C, 5,000 rpm, 1 min). The metabolite-containing supernatant was carefully transferred to a new tube and lyophilized to dryness to concentrate the metabolites and remove volatile solvents. The dried metabolite extract was subsequently reconstituted in a volume of 10% (v/v) methanol solution equivalent to the initial pellet aliquot volume to achieve a consistent concentration suitable for LC–MS/MS injection. Finally, the reconstituted extracts were analyzed by liquid chromatography–tandem mass spectrometry (LC–MS/MS) to characterize the metabolite composition.

Raw mass spectrometry data files (.raw format) acquired from LC–MS/MS analyses were imported into Compound Discoverer (CD) software (version 3.3) for comprehensive preprocessing. Initial data processing included alignment of chromatographic retention times and mass-to-charge ratios (m/z) across all samples to ensure consistent feature detection and matching. To correct for potential instrumental drift and enhance quantitative accuracy, peak areas were normalized relative to the first quality control (QC1) sample injected during the analytical sequence. Putative metabolite identification was then performed by matching the experimentally observed molecular ion masses ([M + H]⁺, [M−H]⁻, etc.) and characteristic fragmentation patterns against spectral libraries within the mzCloud, mzVault, and Masslist databases. Signals originating from solvents, reagents, or background contaminants were removed by subtracting features consistently detected in blank (solvent-only) samples run alongside the biological samples. For relative quantification, the raw peak area of each metabolite in each sample was normalized using the following formula to account for variations in overall metabolite abundance and instrumental response: Relative Peak Area = (Raw Peak Area of Metabolite in Sample)/(Total Sum of All Metabolite Peak Areas in Sample/Total Sum of All Metabolite Peak Areas in QC1). This normalization yielded the relative abundance of each metabolite. To ensure data quality and reproducibility, metabolites exhibiting a coefficient of variation (CV) of relative peak area exceeding 30% across the replicate QC samples were excluded from further analysis. The final curated dataset comprised the confidently identified metabolites and their corresponding relative quantitative values, providing the basis for subsequent statistical comparisons between strains. Processed metabolomics data underwent multivariate analysis using the metaX software package, where log-transformed and scaled data were subjected to Partial Least Squares-Discriminant Analysis (PLS-DA) to derive the Variable Importance in Projection (VIP) score for each metabolite. For univariate analysis, Student’s t-test was applied to calculate the statistical significance (p-value) of intergroup differences for individual metabolites, while the fold change (FC) between groups was computed to quantify relative abundance differences. Differentially abundant metabolites (DAMs) were defined using the following thresholds: VIP > 1.0, FC > 1.2 or FC < 0.833, and p < 0.05.

Identified metabolites were functionally annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG), Human Metabolome Database (HMDB), and LIPID MAPS databases to assign putative biological roles and pathway associations. For multivariate analysis, the processed metabolomics data were log-transformed (or otherwise scaled) and mean-centered using the metaX software package prior to dimensionality reduction. Principal Component Analysis (PCA) was performed to visualize inherent data structure and detect potential outliers. Partial Least Squares-Discriminant Analysis (PLS-DA) was subsequently applied to maximize the separation between the MY1902 parental strain and the ΔqseC mutant groups and to identify metabolites contributing most significantly to this discrimination. Differentially abundant metabolites (DAMs) between the two strains were selected based on the following combined criteria: a Variable Importance in Projection (VIP) score > 1 from the PLS-DA model, a p < 0.05 derived from univariate statistical testing (e.g., Student’s t-test or Wilcoxon rank-sum test), and an absolute fold change (FC) ≥ 2. Pairwise Pearson correlation coefficients among the DAMs were calculated using the cor () function in R. The statistical significance of these correlations was assessed using the cor.mtest () function in R, with a p < 0.05 considered significant. To elucidate the biological implications of the observed metabolic alterations, functional enrichment analysis of the DAMs was conducted based on KEGG pathways. A pathway was considered enriched if the proportion of DAMs mapping to it (x/n) significantly exceeded the proportion of all identified metabolites mapping to the same pathway (y/N), where: x = Number of DAMs mapped to the pathway, n = Total number of DAMs, y = Number of all identified metabolites mapped to the pathway, N = Total number of identified metabolites. The statistical significance of pathway enrichment was evaluated using the hypergeometric test, with a hypergeometric p < 0.05 indicating significantly enriched pathways.

Exponentially growing cultures (OD600 ≈ 0.8) of the G. parasuis parental strain (MY1902) and the isogenic ΔqseC mutant were harvested by centrifugation (4 °C, 4,000 rpm, 5 min) to capture transcriptomes representative of mid-logarithmic phase growth. Total RNA was isolated from the bacterial pellets using the OMEGA Bio-Tek Bacterial Total RNA Extraction Kit according to the manufacturer’s protocol to obtain high-quality, genomic DNA-free RNA suitable for sequencing. RNA integrity and concentration were rigorously assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, United States) to ensure sample quality met the stringent requirements for downstream strand-specific library construction.

Ribosomal RNA (rRNA) was depleted from qualified total RNA samples using a prokaryote-specific rRNA removal kit to enrich messenger RNA (mRNA) and non-coding RNAs. The purified rRNA-depleted RNA was fragmented under optimized conditions. First-strand cDNA synthesis was performed using random hexamer primers. To enable strand-of-origin determination, dUTP was incorporated during second-strand cDNA synthesis instead of dTTP. The resulting double-stranded cDNA underwent end repair, 3′-end adenylation (A-tailing), and ligation of Illumina-compatible adapters. Adapter-ligated fragments within the target size range (typically 250–450 bp) were size-selected. Uracil-containing second strands were enzymatically digested using USER enzyme (NEB) to generate strand-specific libraries. The libraries were then amplified with index primers via limited-cycle PCR and purified. Final library concentration was quantified using Qubit fluorometry and qPCR, while fragment size distribution was verified using an Agilent Bioanalyzer. Libraries passing quality control (appropriate size distribution, absence of adapter dimers, sufficient concentration) were pooled in equimolar ratios based on qPCR-derived effective concentrations to achieve uniform sequencing depth across samples during paired-end sequencing (2×150 bp) on an Illumina NovaSeq 6,000 platform.

The method for standardizing differences in transcriptomic data in this study is TMM. The differential analysis software is edgeR, and the filtering threshold is: padj≤0.005, |log2FoldChange| ≥1.0. Bioinformatic analysis commenced with alignment of filtered sequencing reads to the G. parasuis reference genome using Bowtie 2 to facilitate transcript reconstruction. Rockhopper software assembled these alignments, identifying novel transcript regions through comparison with annotated gene models; putative functions for novel transcripts were assigned via Blastx against the non-redundant (nr) database. Gene expression quantification (FPKM metric) was performed using HTSeq v0.6.1. Differential expression analysis between MY1902 and ΔqseC strains was conducted with DESeq2, defining differentially expressed genes (DEGs) as those with Benjamini-Hochberg adjusted p < 0.05 to control false discovery. Functional enrichment analysis of DEGs elucidated the biological impact of qseC deletion: GOseqR identified enriched Gene Ontology terms (p < 0.05, length-bias corrected), KOBAS revealed significant KEGG pathways, and protein–protein interaction networks were reconstructed from DEGs using the STRING database to uncover functional associations.

The qRT-PCR was performed to validate RNA-seq results by assessing transcription levels of eight randomly selected differentially expressed genes in G. parasuis strains MY1902 and ΔqseC. Bacterial cultures were harvested at mid-log phase (OD₆₀₀ = 0.8), followed by total RNA extraction and reverse transcription into cDNA. Using 16S rRNA as the endogenous control (primers in Table 1), reactions were conducted in technical triplicates on a LightCycler^® 96 system (Roche) with SYBR Premix Ex Taq II (Takara Bio) under standardized cycling conditions: 95 °C/30s; 45 cycles of 95 °C/5 s, 54 °C/30s, 72 °C/30s; and melting curve analysis. Gene expression fold-changes were calculated via the 2⁻ΔΔCT method to confirm transcriptomic differences between strains.

Integrated transcriptomic-metabolomic analysis was conducted to identify key post-transcriptional regulators (e.g., RNA-binding proteins, non-coding RNAs) governing metabolic flux in G. parasuis. This approach correlated RNA-level expression patterns with small-molecule metabolite dynamics to elucidate system-level adaptations resulting from qseC deletion.

All experimental data were analyzed using GraphPad Prism 6.0, with continuous variables assessed by Student’s t-test (two-group) or one-way ANOVA (multi-group), survival data analyzed via log-rank (Mantel-Cox) test, and statistical significance defined as p < 0.05 (denoted by asterisks).

As QseC is a membrane-localized sensor kinase, we first interrogated whether its deletion compromises cellular integrity using TEM-the gold standard for ultrastructural analysis. Wild-type MY1902 exhibited intact cell walls and membranes without cytoplasmic separation, with electron-dense cytoplasmic inclusions observed (Figure 1A, red arrows). In contrast, ΔqseC mutants showed focal wall/membrane dissolution (Figure 1B, red arrows), cytoplasmic separation (Figure 1B, yellow arrows), and protoplast disintegration. These defects indicate QseC was essential for maintaining envelope integrity and cytoplasmic organization in G. parasuis, validating our structural stability hypothesis.

Given the structural defects observed, we next performed untargeted metabolomics to test whether QseC deletion disrupts metabolic networks supporting envelope biogenesis and stress adaptation. We performed untargeted metabolomics using UHPLC-Q-TOF MS. Data were deposited in MetaboLights (MTBLS12504).

Based on QseC’s predicted function as a signal transduction kinase, we hypothesized it directly regulates transcriptional programs for envelope maintenance and stress adaptation. This hypothesis was tested through comparative transcriptomics of WT and ΔqseC strains. The raw transcriptome data have been deposited in the National Center for Biotechnology Information (NCBI) database under the accession number PRJNA1297586.

We integrated transcriptomic and metabolomic data to identify regulatory networks affected by qseC deletion. Pearson correlation analysis revealed significant associations between differentially expressed genes and metabolites. The top 50 genes and metabolites with strongest correlations were shown in Figures 7A,B. Key relationships between the top 10 DEGs and top 5 metabolites were displayed in Figures 7C,D (Supplementary Tables S14, S15). To map shared biological pathways, we analyzed KEGG enrichment. In positive ion mode, co-enriched pathways included amino acid and purine metabolism. In negative ion mode, enriched pathways involved cell wall synthesis and lipid metabolism (Figures 8A,B, Supplementary Tables S16, S17). Integrated visualizations showed coordinated changes in these pathways (Figures 8C,D). These results demonstrated that QseC synchronizes transcriptional and metabolic programs to maintain structural-metabolic coupling.

Environmental stressors activated the histidine kinase QseC, triggering phosphorylation of response regulator QseB. This phosphorylation drove transcriptional reprogramming with upregulation of membrane genes plsB and wecA, and downregulation of hrpA and pilW. Concurrent metabolic changes showed increased epigallocatechin and KQH alongside decreased thromboxane B2. Transmission electron microscopy confirmed membrane damage in ΔqseC mutants, linking this network to cell integrity maintenance. Loss of QseC function caused impaired biofilm formation, reduced stress resistance and compromised osmotic tolerance. Figure 9 illustrated this QseC-dependent regulatory pathway, confirming its role in bacterial adaptability through coordinated transcriptomic-metabolic control. The cell signal transduction pathway map was created with BioGDP.com (Jiang et al., 2024).