Integrative Multi-omics Analysis of Buti Huatan Tang in Chronic Obstructive Pulmonary Disease

Summary

Multi-omics profiling of Buti Huatan Tang (BTHTT) in chronic obstructive pulmonary disease reveals its potential to modulate the IL-1R2/IL-1β-thyroid hormone metabolic axis. This work provides a framework for further mechanistic validation of the bioactive components of BTHTT, establishing a foundation for future causal intervention studies required to confirm this axis.

Abstract

This study utilized a multi-omics and computational biology framework to investigate the therapeutic potential of the Traditional Chinese Medicine (TCM) formula Buti Huatan Tang (BTHTT) against chronic obstructive pulmonary disease (COPD). Significant physiological improvements were observed in a rat model following BTHTT intervention. Histological analysis showed a reversal of lung pathological damage, while biochemical assays, and transcriptomics confirmed the normalization of IL-1β and IL-1R2 levels. Additionally, metabolic profiling revealed that BTHTT corrected disruptions in T3 and T4 thyroid hormone levels. A negative correlation was observed between the IL-1β/IL-1R2 axis and these thyroid hormones, indicating that their regulation is associated with the formula’s therapeutic effect. Beyond direct measurements, machine learning algorithms identified ten COPD signature genes from clinical databases. Pathway enrichment analysis suggests that BTHTT may act through cytokine-cytokine-receptor interactions and thyroid hormone synthesis pathways. Furthermore, while 283 components were identified in vivo, compounds such as tanshinone IIA and cryptotanshinone are currently considered candidate active substances. Their role as primary drivers is supported by a model in which they stably bind to IL-1R2; this inference is based on molecular docking and molecular dynamics (MD) simulations rather than direct experimental isolation. Overall, the data support a model in which BTHTT exerts a multi-target effect on COPD by modulating inflammation and metabolic homeostasis. This integrated approach provides a refined scientific basis for the clinical application of BTHTT and highlights specific pathways for future experimental validation.

Introduction

Chronic obstructive pulmonary disease (COPD), one of the most common respiratory diseases in clinical practice, has become the third leading cause of death globally due to its extremely high disability rate and recurrent acute exacerbations¹. COPD was a significant factor leading to the loss of labor capacity in the population and a sharp increase in family care costs, with social and economic losses far exceeding those of most chronic diseases². In the “Healthy China Initiative 2030,” COPD has been listed as a key disease for prevention and treatment. Due to an incomplete understanding of the mechanisms underlying COPD, the combined use of bronchodilators and corticosteroids remained the primary treatment. While this approach provided relatively obvious symptom relief, it cannot prevent the progressive decline in lung function. Moreover, the long-term use of corticosteroids inevitably increases the incidence of adverse reactions, further affecting the long-term efficacy of the drugs and the patients' tolerance³. Studies have shown that even with triple therapy, 28%–31% of patients still experience acute exacerbations or worsening of symptoms. Therefore, the development of effective drugs for COPD remains a key and primary issue in current COPD research⁴.

To address these needs, classic TCM formulas such as Buti Huatan Tang (BTHTT), Shengmai Yin, and Bu Fei Tang have been developed. Among them, BTHTT was particularly notable for its significant effects in tonifying the lung, boosting energy, and resolving phlegm to relieve cough^5,6,7,8. BTHTT was widely used in the clinical treatment of COPD, not only for its remarkable efficacy and minimal side effects, but also for its potential to reverse the progressive decline in lung function. It held great promise in clinical and social contexts and had the potential for “secondary development” as a major drug, with significant scientific and economic value^9,10,11,12. However, due to the lack of systematic, long-term scientific research, studies on BTHTT have generally suffered from an unclear understanding of its component basis and unidentified core targets of action. Fundamentally, this was due to significant deficiencies in research on the main mechanism and potential therapeutic substance basis of BTHTT.

In response to this characteristic, integrating multiple omics technologies such as metabolomics and transcriptomics, as well as multi-algorithm mining methods, focusing on the main correlations between characteristic key genes and key metabolic pathways in the body's microenvironment^13,14,15 cannot only explore the essence of treatment but also evaluate the synergistic efficacy of the active components. A particular focus is placed on the thyroid hormone axis, as emerging research indicates that thyroid dysfunction is a frequent comorbidity in COPD, often correlating with increased systemic inflammation and impaired respiratory muscle function. Modulating this axis may represent a novel pathway for restoring metabolic homeostasis in pulmonary disease^16,17.

In summary, this study adopted a three-pronged approach combining multi-omics integration with computational biology, based on the premise of syndrome and formula correspondence. On the one hand, it utilized TCM serum pharmacochemistry techniques to identify the blood-borne components of BTHTT under effective conditions. On the other hand, focusing on key characteristic genes of COPD, it used comprehensive tissue metabolomics to discover disease biomarkers from endogenous small molecules and to construct a precise evaluation of the overall effects of the formula. Meanwhile, combined with multi-algorithm mining and transcriptomic identification analysis, it identified the key characteristic genes regulated by BTHTT under effective conditions and conducted correlation analysis between the key characteristic genes and key metabolic pathways. The potential therapeutic substance basis was screened by molecular docking and molecular dynamics analysis.

Based on our multi-omics framework and the pharmacological profile of BTHTT, this study proposes the following testable predictions. Antiinflammatory modulation: BTHTT intervention will significantly reduce systemic and pulmonary inflammation markers, specifically by normalizing the IL-1β/IL-1R2 expression axis to mitigate tissue damage. Metabolic restoration: BTHTT will alter signals in the thyroid hormone synthesis pathway, leading to a measurable correction of T3/T4 levels and a subsequent improvement in the metabolic state of the lung microenvironment. Metabolic Restoration: BTHTT will alter signals in the thyroid hormone synthesis pathway, leading to a measurable correction of T3/T4 levels and a subsequent improvement in the metabolic state of the lung microenvironment.

Protocol

NOTE: The reagents and the equipment used in this study are listed in the Table of Materials.

Clinical data acquisition and preprocessing
To investigate the molecular signatures of COPD, transcriptomic datasets were retrieved from the Gene Expression Omnibus (GEO) database. A total of five datasets were selected: GSE11784, GSE12472, GSE16972, GSE38974, and GSE222965. The raw data and platform annotation files were downloaded for probe-to-gene mapping. When multiple probes targeted a single gene, the maximum expression value was retained. The resulting gene expression matrices (rows as genes, columns as samples) were merged into a single discovery dataset. To account for technical variations across different microarray platforms and study cohorts, batch-effect correction was performed using the ComBat algorithm from the sva R package. The efficacy of the correction was validated via Principal Component Analysis (PCA) plots. Following batch correction, the merged discovery dataset was utilized for differential analysis. Differentially expressed genes (DEGs) between COPD patients and healthy controls were identified using the limma package. The significance thresholds were set at |logFC| ≥ 1 and p ≤ 0.05.

To identify the most robust characteristic genes, three independent machine learning algorithms were integrated. To ensure the reliability of the models and prevent data leakage, the feature selection process was nested within the cross-validation loops where applicable, and the discovery dataset was strictly separated from the independent validation sets. LASSO model was applied to the DEGs using the glmnet package. We employed 10-fold cross-validation to determine the optimal penalty parameter. The optimal penalty parameter value corresponding to the minimum cross-validation error was selected as the threshold to identify core feature genes. SVM was utilized to rank genes based on their discriminative power. A 10-fold cross-validation strategy was implemented to identify the point of minimum generalization error, thereby determining the optimal number of characteristic genes. The randomForest package was used to rank DEGs based on their Mean Decrease Accuracy and Gini index. Genes with the highest importance scores were selected as disease-related features.

The intersection of the features identified by LASSO, SVM, and randomForest was taken to define the final core characteristic genes. The diagnostic performance of these genes was evaluated using Receiver Operating Characteristic (ROC) curve analysis. The Area Under the Curve (AUC) and associated 95% Confidence Intervals (CIs) were calculated using the pROC package. A gene was considered to have high diagnostic value if the AUC > 0.70. Finally, the expression levels and diagnostic accuracy of these genes were further validated in an independent validation set to ensure the generalizability of the findings.

Using the PubChem (https://pubchem.ncbi.nlm.nih.gov/) and drugbank (https://go.drugbank.com/), a total of 12 specific drug profiles were obtained, such as LABAs, LAMAs, SABAs, SAMAs, Fluticasone propionate, Budesonide, Beclomethasone, Fluticasone, Salmeterol, Umeclidinium, Vilanterol, Theophylline. The compiled targets were calibrated using the Uniprot database (https://www.uniprot.org/), during which non-human genes were removed, and invalid duplicate targets were deleted to obtain standardized gene names. By entering the keyword “chronic obstructive pulmonary disease”, “COPD” in the GeneCards (https://www.genecards.org/), OMIM (https://www.omim.org/), and DisGENET (https://www.disgenet.org/) databases, disease-related targets were retrieved. All targets from the three databases were consolidated into an Excel file, duplicate genes were removed, and the data was calibrated using the Uniprot database to obtain the final disease target gene information.

Machine learning model development
A multi-algorithmic machine learning framework was constructed by sequentially utilizing R packages glmnet, e1071, and randomForest. Specifically, LASSO regression was performed for penalty-based dimensionality reduction, Support Vector Machine (SVM) analysis was employed to evaluate validation errors based on sample grouping, and Random Forest was used to filter features according to importance scores. This process yielded corresponding diagnostic visualizations, including cross-validation curves and gene importance bubble plots. Following model construction, a Venn diagram analysis was performed on the gene sets identified by these multiple algorithms to extract the overlapping "intersection" feature genes, thereby enhancing the reliability of the potential biomarkers. The expression matrix of these intersection genes was then extracted to visualize inter-group expression differences via violin plots. Finally, ROC curves were generated through iterative loops for each gene to calculate the Area Under the Curve (AUC), validating their diagnostic value as candidate biomarkers.

Animal model induction
The experimental protocol was approved by the Animal Ethics Committee of Guizhou Medical University (2303411) and followed the ARRIVE guidelines and animal welfare regulations. A total of 72 male Sprague-Dawley (SD) rats of SPF grade (230 ± 20 g). After 1 week of adaptive housing under standard conditions (25 ± 1 ℃, 50 ± 5% humidity, 12 h light-dark cycle), the rats were randomly divided into the Control group (C group) (n = 24) and the Model group (M group) (n = 48).

The M group was subjected to dual-factor modeling: Intermittent cigarette smoke exposure (9 weeks, 6 days per week, 3 research-grade cigarettes per day divided into 2 sessions, 30 min per session) and intratracheal instillation of LPS (200 µg per instillation, once every other week)^18,19. The C group received an equivalent volume of normal saline. After successful modeling, the groups were randomly divided into the M group (n = 12), the high-dose BTHTT group (H group) (High, 1× clinical dose, n = 12), and the low-dose BTHTT group (L group) (Low, 1/2× clinical dose, n = 12), with continuous gavage intervention for 2 weeks. The C and M groups were given distilled water synchronously.

Body weight and spontaneous activity parameters were recorded weekly. At the end of modeling and treatment, lung tissue, serum, and serum and bronchoalveolar lavage fluid (BALF) were collected. Levels of inflammatory markers (i.e., TNF-α, IL-1β, IL-6, IL-8, OPN, and MCP-1) were measured using ELISA according to the kit instructions. Tissue pathology and multi-omics analyses were also performed.

Treatment administration
Buti Huatan Tang (BTHTT) consists of nine traditional Chinese medicinal herbs: Astragalus membranaceus (Fisch.) Bge. (Astragali Radix, 15 g), Pseudostellaria heterophylla (Miq.) Pax (Pseudostellariae Radix, 15 g), Cinnamomum cassia (L.) J. Presl (Cinnamomi Cortex, 15 g), Angelica sinensis (Oliv.) Diels (Angelicae Sinensis Radix, 10 g), Salvia miltiorrhiza Bge. (Salviae Miltiorrhizae Radix et Rhizoma, 15 g), Perilla frutescens (L.) Britt. (Perillae Folium, 10 g), Raphanus sativus L. (Raphani Semen, 10 g), Lepidium apetalum Willd. (Descurainiae Semen, 10 g), and Mahonia fortunei (Lindl.) Fedde (Mahoniae Folium, 10 g). The herbal materials were soaked in water at 10x their combined weight for 30 min and subsequently decocted for 1 h. The decoction was filtered, and the filtrate was collected and divided into three equal portions for oral administration.

Based on clinical safety and efficacy, the standard dose of BTHTT for adults was 1.57 g∙kg^-1∙day^-1 according to the clinical medication guidelines. Considering the dose conversion factor for rats, which indicated that the standard rat drug dose was 6.3× the human standard dose, the corresponding H group for rats was set at 9.9 g∙kg-1∙day-1, and the L group was set at half this dose (4.95 g∙kg^-1∙day^-1). Given the final concentrated volume of BTHTT was 50 mL, the administered volume for the H group in rats was approximately 1.6 mL, and for the L group, approximately 0.8 mL. The drug was administered once daily via oral gavage.

Tissue and BALF collection
At the experimental endpoints (Week 9 and Week 11), rats were anesthetized via intraperitoneal injection of 5% sodium pentobarbital (1 mL/100 g). Blood was collected from the portal vein, allowed to stand for 30 min at room temperature, and then centrifuged at 13,000 × g for 15 min at 4 ℃. The supernatant was stored at -80 ℃. Lung tissues were rapidly frozen in liquid nitrogen and stored at -80 ℃. BALF was collected by three consecutive lavages with cold PBS.

Transcriptomics analysis
Total RNA was extracted from lung tissue and quality-controlled using (A₂₆₀/A₂₈₀ > 1.8) (RIN ≥ 7.0. RNA sequencing libraries were constructed: mRNA was enriched with oligo(dT), fragmented, and used to synthesize double-stranded cDNA, which was then ligated to adapter oligonucleotides and amplified by PCR. Libraries were quantified and quality-checked before sequencing. Including content distribution inspection, FPKM density distribution analysis of each sample, and overall quality assessment analysis of RNA-seq²⁰.

Raw transcriptomics data generated by the sequencing platform were processed using Perl scripts to remove adapter sequences and low-quality reads (reads with Q ≤ 25 bases accounting for > 60% or N rate > 5%). Clean reads were obtained after this filtering process. The clean reads were aligned to the reference genome using HISAT2, and gene expression levels were quantified to calculate FPKM values. Differential expression analysis was performed (with screening criteria of |log2FC| > 1 and p-adj < 0.05). Transcription factor annotation was based on Animal TFDB or Pfam/DBD databases, matching gene IDs and protein domain information²¹.

An orthologous gene mapping strategy was employed to ensure methodological rigor for cross-species validation. This strategy involved retrieving orthologs in rats for core human genes (e.g., SRPX2, IL-1R2, TFF3) using the NCBI HomoloGene and Ensembl BioMart databases. The selection was restricted to gene pairs exhibiting a clear "one-to-one" mapping relationship and high protein sequence identity. In instances where multiple candidates were present, preference was given to orthologous pairs certified by the HGNC. To ensure high detection precision, specific RT-qPCR primers were designed based on the mRNA sequences of the identified rat orthologs. Validation criteria were defined by consistency in expression direction and functional verification. Consistency in Expression Direction: In a smoke-induced rat COPD model, RT-qPCR showed that the expression trends of target genes in rat lung tissue were fully consistent with those observed in human GEO clinical datasets. Genes exhibiting the same polarity of change across species were considered to possess conservation as disease biomarkers. Pathological Verification: Upon expression consistency, correlation analysis was performed to verify the involvement of these genes in the pathological evolution of COPD.

Metabolomics workflow
Lung tissue (20–50 mg) was homogenized in prechilled methanol-acetonitrile-water (2:2:1, v/v) and sonicated in an ice bath. The homogenate was centrifuged 13,000 × g for 20 min at 4 ℃. The supernatant was vacuum-concentrated, re-dissolved in acetonitrile-water (1:1, v/v), and filtered through a 0.22 µm membrane for LC-MS analysis²².

Chromatographic separation was achieved using an amide-based UPLC column (1.7 µm, 2.1 mm × 100 mm). Column temperature was maintained at 25 ℃. The mobile phase consisted of (A) water containing 25 mM ammonium acetate and 25 mM ammonia, and (B) acetonitrile. Flow rate was set at 0.5 mL/min, and the injection volume was 2 µL. Gradient elution program was as follows: 0–0.5 min, 95% B; 0.5–7 min, linear decrease of B from 95% to 65%; 7–8 min, linear decrease of B from 65% to 40%; 8–9 min, B held at 40%; 9–9.1 min, linear increase of B from 40% to 95%; 9.1–12 min, B held at 95%. During the entire analysis, samples were kept at 4 ℃ in the autosampler. To ensure system stability and the reliability of the experimental data, samples were analyzed in a random sequence, with quality control (QC) samples interspersed in the queue. Mass spectrometry analysis was performed using an ultra-high-performance liquid chromatography (UHPLC) system coupled to a mass spectrometer. Samples were ionized using electrospray ionization (ESI) in both positive and negative ion modes. ESI source and MS settings were as follows: nebulizer gas (Gas 1) was set to 50, auxiliary gas (Gas 2) to 2, ion source temperature to 350 ℃, and spray voltage (ISVF) to 3,500 V in positive ion mode and 2,800 V in negative ion mode. The mass range for MS1 was set from 70 to 1,200 Da, with a resolution of 60,000 and a scan accumulation time of 100 ms. For MS2, data-dependent acquisition (DDA) with stepped collision energy was used. The mass range for MS2 was also set from 70 to 1,200 Da, with a resolution of 60,000 and a scan accumulation time of 100 ms. The dynamic exclusion time was set at 4 s.

Raw metabolomics data were converted to the mzXML format and then processed for peak alignment, retention time correction, and peak area extraction. Data preprocessing workflow included the following steps: First, ion peaks with a missing rate > 50% were removed. Second, the remaining missing values were imputed using the KNN algorithm. Third, metabolic features with a relative standard deviation (RSD) >50% were discarded. The quality of the experimental data was assessed using principal component analysis (PCA) and clustering of QC samples. Subsequent analyses included univariate statistics (e.g., t-tests), multivariate statistics (PLS-DA), differential metabolite screening (VIP > 1 and p < 0.05), and KEGG pathway enrichment analysis (hypergeometric test)^23,24.

Establishment of analytical methods for in vivo and in vitro components
Preparation of BTHTT samples for in vitro testing: BTHTT was extracted by decoction in water (2 x 30 min), concentrated to 1.1–1.2 g/mL, and then freeze-dried. Before analysis, 600 µL of the freeze-dried powder solution was mixed with 400 µL of methanol, re-dissolved in 40% methanol, and centrifuged to collect the supernatant.

Preparation of BTHTT samples for in vivo testing: Serum was deproteinized by mixing with methanol (1:1) and precipitating at −20 ℃ for 30 min, followed by centrifugation at for 20 min. The supernatant was vacuum-dried and re-dissolved in 40% methanol to obtain the final sample. For the preparation of blank serum + BTHTT samples, an appropriate amount of blank serum was spiked with the in vitro BTHTT supernatant, and the remaining steps were performed as described.

Samples were separated using a UHPLC system equipped with a reverse-phase UPLC column (2.1 mm × 100 mm, 1.8 µm). The column temperature was maintained at 35 °C, and the flow rate was set at 0.3 mL/min. The mobile phase consisted of (A) 0.1% formic acid in water and (B) 0.1% formic acid in acetonitrile. Gradient elution was performed as shown in Table 1.

A mass spectrometer was used for the acquisition of both MS1 and MS2 spectra. The mass spectrometer was coupled with the UHPLC system and operated in both positive and negative ESI modes. ESI parameters were as follows: spray voltage 3,800 V (ESI+) / 3500 V (ESI-), sheath gas pressure 45 arb, auxiliary gas pressure 20 arb, ion transfer tube temperature 320 °C, and vaporizer temperature 350 °C. The detection mode was set to full scan/data-dependent MS2 (Full-MS/dd-MS2) with resolutions of 60,000 for MS1 and 15,000 for MS2. The top 10 MS1 ions were selected for MS/MS fragmentation with stepped normalized collision energies of 20, 40, and 60. Mass range for MS1 was set from 90 to 1,300 Da.

For in vivo analysis, including blank group samples, dosed group samples, and blank group + BTHTT samples, 6 µL of each sample was precisely injected. For in vitro analysis of BTHTT, 2 µL of the sample was injected. Each batch of blank and dosed group samples was injected once, while the blank group + BTHTT samples were injected in triplicate, and the BTHTT samples were injected in quintuplicate.

Data in the mzXML format were processed and compounds were identified based on a local high-resolution commercial TCM mass spectrometry database. The criteria for identification were set as follows: a mass error of <25 ppm for MS1 and a match score > 0.7 for MS2 (where the Score reflected the similarity of fragment ions, with ≥0.7 being a reliable threshold)^25,26. Statistical analysis included compound counting and classification (e.g., flavonoids, alkaloids), which was performed in conjunction with annotations from the mass spectrometry database²⁷.

Molecular docking and MD simulation
To investigate potential binding modes between the identified characteristic proteins and their corresponding ligands, in silico molecular docking was performed. The three-dimensional structures of the small molecules were retrieved from the PubChem database, and their geometric configurations optimized. The crystal structures of the target proteins were obtained from the RCSB Protein Data Bank (PDB). Using PyMOL, water molecules and heteroatoms were removed, and co-crystallized ligands were extracted to define the active site coordinates. Hydrogen atoms were added, and Gasteiger charges were assigned using software. Docking simulations were executed, generating 15 independent conformations per run. The conformation with the lowest binding energy was selected for further analysis. To rigorously characterize the non-covalent interactions, the receptor-ligand complexes were analyzed using the Protein-Ligand Interaction Profiler (PLIP). NOTE: It is important to emphasize that these docking results provide structural support for potential molecular interactions and serve as a foundation for further dynamic refinement; however, they do not constitute standalone evidence of biological efficacy.

To evaluate the stability and conformational evolution of the predicted protein-ligand complexes under physiologically relevant conditions, molecular dynamics simulations were performed using the GROMACS software package. Topology files for both the proteins and ligands were generated based on the GROMOS96 43a1 force field. Each complex was positioned at the center of a dodecahedral box, maintaining a minimum distance of 1.0 nm from the box edges, and solvated using the SPC water model. To ensure electrical neutrality, sodium or chloride ions were added to the system as required. Energy minimization was performed using the steepest descent algorithm until the maximum force was less than 1,000.0 kJ∙mol^-1∙nm^-1. The system was then equilibrated in two stages: first, an NVT ensemble was employed to heat the system to 300 K over 100 ps using a V-rescale thermostat; second, an NPT ensemble was used to stabilize the pressure at 1 bar over 100 ps using a Parrinello-Rahman barostat. Production simulations were carried out for a total duration of 10 ns with a time step of 2 fs. Long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method, while short-range van der Waals and electrostatic interactions were managed with a cutoff radius of 1.2 nm. To ensure the reliability of the simulations, three independent runs were performed where feasible. The stability of the complexes was quantitatively assessed by calculating the Root Mean Square Deviation (RMSD) and Root Mean Square Fluctuation (RMSF) of the protein backbone atoms relative to the initial structure. The attainment of a plateau in the RMSD profile was utilized as the primary criterion for system equilibration and structural stability.

General statistical analysis
In this experiment, group calculations were performed using t-tests or one-way analysis of variance (ANOVA).

Results

Prediction of key characteristic genes in COPD
Following batch-effect correction of clinical data from multiple microarray platforms (Supplemental File 1–Supplemental Figure S1), a differential expression analysis was performed. 34 candidate genes (32 upregulated and 2 downregulated) were identified, including SRPX2, CLDN10, and TFF3 (Figure 1A).

To refine these candidates into a robust predi...

Discussion

The synthesis and maintenance of thyroid hormone levels are significant for COPD patients and extend far beyond mere endocrine function, profoundly influencing disease progression, nutritional status, pulmonary function, and overall prognosis^31,32,33. Particularly during acute exacerbations, COPD patients often exhibited a hypothyroidism-like state referred to as “low T3 syndrome”. It was reported that decreased thyr...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Acetonitrile (chromatographic-grade)	Merck, Germany	1499230-935
Agilent 4150 bioanalyzer	Agilent Technologies, USA
Ammonium acetate	SIGMA, Germany	SRE0035-20230321
Angelica sinensis (Oliv.) Diels (Angelicae Sinensis Radix),	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240807
Astragalus membranaceus (Fisch.) Bge. (Astragali Radix),	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240601
ChemOffice		chemical structure drawing and molecular modeling software
Cinnamomum cassia (L.) J. Presl (Cinnamomi Cortex),	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240901
Eppendorf Centrifuge 5430 R	Eppendorf, Germany
Experimental cigarettes	Harbin Taihua, China	20230208
GraphPad Prism 8		statistical analysis and graphing software
Huangshan cigarettes		research-grade cigarettes for smoke exposure experiments
IL-1β ELISA kit	Shanghai Enzyme-Linked Biotechnology, China	E-EL-R0012c-20230522
Illumina Novaseq 6000 sequencing platform	Illumina, USA
Lepidium apetalum Willd. (Descurainiae Semen)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240201
Lepidium sativum L. (Raphani Semen)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240301
Lipopolysaccharide (LPS)	Shandong Huino Pharmaceutical, China	20230315
Mahonia fortunei (Lindl.) Fedde (Folium Mahoniae)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240301
Methanol (chromatographic-grade)	Fisher, USA	DK745764
Nanodrop ND-2000 spectrophotometer	Thermo Scientific, USA
Perilla frutescens (L.) Britt. (Perillae Folium)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240510
Pseudostellaria heterophylla (Miq.) Pax (Pseudostellariae Radix),	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20250201
Qubit 4.0 fluorometer	Thermo Fisher Scientific, USA
RNA library construction kit	ABclonal, USA	RK20301-20230415
SMINA		molecular docking software
Salvia miltiorrhiza Bge. (Salviae Miltiorrhizae Radix et Rhizoma)	Guizhou Shanhaitongjia Pharmaceutical Co., Ltd., China	20240902
Small animal spontaneous activity recorder (KW-ZF)	Nanjing Karlvin Biotechnology, China
Sodium pentobarbital	Shanghai Experimental Reagent Procurement Company, China	20230412
TNF-α ELISA kit	Shanghai Enzyme-Linked Biotechnology, China	E-EL-R2856c-20230522
Traditional Chinese Medicine MS database (Shanghai Applied Protein Technology Co., Ltd.)		commercial traditional Chinese medicine mass spectrometry spectral database
Trizol reagent	Magen, China	20230508
UPLC BEH Amide column		amide-based UHPLC analytical column
UPLC HSS T3 column		reversed-phase UHPLC analytical column (T3-type)
Vanquish Neo UHPLC-Orbitrap Exploris 480 MS system	Thermo Scientific, USA
Vanquish UHPLC-Q-Exactive HF-X MS system	Thermo Scientific, USA