Identification of IPF specific genes expressed by lung tissue-derived fibroblasts using next generation sequencing.

Idiopathic pulmonary fibrosis (IPF) is a progressive interstitial lung disease characterized by aberrant activation of fibroblasts and excessive deposition of extracellular matrix (ECM) components such as collagens, elastin, and glycoproteins. These changes disrupt alveolar architecture, increase lung stiffness, and reduce lung volume, a key predictor of mortality ( 1 , 2 ). In normal wound healing, myofibroblasts are transiently activated and undergo apoptosis after tissue repair. The ECM is largely degraded and replaced with normal tissue, and myofibroblasts are rarely observed in healthy lungs ( 3 , 4 ). In contrast, IPF is thought to result from aberrant wound healing following repetitive alveolar epithelial injury. Rather than resolving, this process leads to persistent fibroblast activation and progressive fibrosis ( 1 , 5 ). A key component of this dysfunctional repair is the excessive proliferation of fibroblasts and their differentiation into contractile myofibroblasts, forming fibroblastic foci—a histological hallmark of IPF associated with worse survival ( 6 ). Myofibroblast accumulation is driven by both external cues and intrinsic cellular abnormalities. Alveolar epithelial damage triggers the release of profibrotic cytokines, including interleukin-1 (IL-1), tumor necrosis factor-alpha (TNF-α), interleukin-6 (IL-6), and transforming growth factor-beta (TGF-β) ( 7 ). In parallel, fibroblasts from IPF lungs exhibit phenotypic alterations, including elevated expression of contractile proteins and enhanced ECM synthesis, relative to fibroblasts from healthy tissue ( 8 ). Transcriptomic profiling of IPF-derived fibroblasts offers enables elucidation the molecular mechanisms underlying fibroblast activation and fibrosis progression. To date, three studies have characterized differentially expressed genes (DEGs) in fibroblasts isolated from IPF lung tissue: GSE180415 identified 716 DEGs (419 upregulated and 297 downregulated) ( 9 ); GSE17978 reported 1,363 DEGs ( 10 ); and a meta-analysis of GSE1724 , GSE10921 , GSE44723 , and GSE40839 identified 227 DEGs (114 upregulated, 113 downregulated) ( 11 ) as illustrated in the Venn diagram ( Figure S1a ). Nevertheless, only two genes — SECTM1 and RPB1 — exhibited consistent differential expression across all three studies ( Table S1 ), reflecting limited reproducibility. Similarly, GO enrichment analysis yielded 19 overlapping terms between pairs of studies, but none were shared across all three ( Table S2 , Figure S1b ). Such discrepancies may stem from small sample sizes and heterogeneous study designs: Hanmandlu et al analyzed nine samples (5 IPF, 4 controls), and Emblom-Callahan et al used 18 (12 IPF, 6 controls). The Plantier et al combined 40 samples from four datasets. Furthermore, demographic and anatomical differences added variability. These limitations highlight the need for transcriptomic studies with larger, clinically consistent cohorts to improve reproducibility and refine our understanding of IPF pathogenesis. Here, we performed next-generation RNA sequencing (RNA-seq) on cultured lung-derived fibroblasts from a relatively large cohort of 33 IPF patients and 10 controls. We then identified DEGs and performed functional enrichment and network analyses to characterize fibroblast-specific transcriptional changes. Comparative analysis with three previously published IPF fibroblast datasets indicated robust, reproducible ECM-related signatures and uncovered novel candidate genes not previously implicated in IPF, offering new insights into the molecular mechanisms of fibroblast activation and fibrosis progression.

Clinical characteristics of the study participants are summarized in Table 1 . A total of 43 individuals were included in the analysis, comprising 33 patients with IPF and 10 controls. The mean age was higher in the IPF group than in the control group, although this difference did not reach statistical significance ( p =0.081). Sex distribution and smoking status were comparable between the two groups. In contrast, pulmonary function parameters showed significant impairment in the IPF group. The median predicted forced vital capacity was significantly lower in IPF patients compared to controls ( p =0.010). Similarly, the median predicted diffusion capacity for carbon monoxide was markedly reduced in the IPF group relative to controls ( p =0.002). Clinical characteristics of the study subjects. Values are presented as median (interquartile range) or mean ± standard error of the mean. Abbreviations: IPF: idiopathic pulmonary fibrosis; NS/ES/SM: non-smokers, ex-smokers, smokers; FVC: forced vital capacity; DL CO : diffusion capacity of the lung for carbon monoxide. From the initial dataset comprising 46,427 genes, we excluded 31,375 genes with fewer than 10 read counts and applied TMM normalization to the remaining 15,252 genes. Based on the criteria of |FC| ≥ 2 and adjusted p (Bonferroni) < 0.05, a total of 475 genes were identified as differentially expressed, including 402 upregulated and 73 downregulated in IPF fibroblasts. ( Figure S2 ). Volcano plot and hierarchical clustering heat map of these genes are presented in Figure 1 . The complete gene list is provided in Table S3 . The top 20 genes with the highest FC values are listed in Table 2 . Among them, eight genes— BMP5, TMEM176B, C7, GDF10, ADAMTS8, EDNRB, FGFR4 , and TNFSF15 —are associated with the development of IPF ( 20 - 27 ). Five genes-including ADH1B, ADH1A, EDNRB, GDF10, and FGFR4 have been linked to liver fibrosis ( 28 - 32 ), while SCUBE1 and NKD2 are associated with kidney fibrosis ( 33 , 34 ). ADAMTS8 and ADAMTS19 are implicated in cardiac and skin fibrosis, respectively ( 35 , 36 ) and CFTR and TNFSF15 are connected to cystic and colonic fibrosis ( 37 , 38 ). In contrast, HSD17B2 , HMGCLL1 , RASL12 , HBG1 , the noncoding RNA RIPOR3-AS1 , and the predicted gene LOC105375566 have not yet been studied in the context of fibrosis. Top 20 differentially expressed genes with the highest absolute fold change between 33 IPF and 10 control fibroblasts. Values are presented as mean ± standard error of the mean. p -values shown are Bonferroni-adjusted values from the Exact test. Abbreviations: IPF: idiopathic pulmonary fibrosis; log 2 FC: fold change of IPF/Control; FDR: false discovery rate; DOI: digital object identifier; N/A: not applicable. Differential gene expression analysis of lung tissues - derived fibroblasts between 33 IPF patients and 10 controls. (a) Volcano plot displaying differential expression of 15,252 genes in lung fibroblasts, highlighting 402 significantly upregulated (red) and 73 downregulated (blue) genes in IPF compared to controls (cutoff: |fold change| > 2, Bonferroni adjusted p -value < 0.05). (b) Hierarchical clustering heatmap of 475 differentially expressed genes, demonstrating distinct expression profiles. Abbreviations: IPF: idiopathic pulmonary fibrosis; log2FC: log fold change of IPF/control. In gene ontology analysis of the 475 DEGs, 463 were recognized by the DAVID v6.8 database ( https://david.ncifcrf.gov/ ) and included in the functional enrichment analysis. FDR-adjusted p-values < 0.05 considered significant. A total of 21 GO terms were identified across the three major GO categories: 8 in BP, 8 in CC, 1 in MF and 4 in KEGG pathway. Out of 402 up-regulated DEGs, 392 were used for GO analysis. A total of 21 GO terms were identified across the three major GO categories: 7 in BP, 8 in CC, 2 in MF, 4 in KEGG, and 4 in REACTOME pathway. Out of 73 down-regulated DEGs, 71 were used for GO analysis. None of GO terms were identified. The full list of GO annotations, along with associated gene sets, is provided in Table S4 . Both 475 total DEGs and 402 up-regulated DEGs presented the top five significantly enriched GO terms as extracellular region (GO:0005576), extracellular space (GO:0005615), plasma membrane (GO:0005886), collagen-containing extracellular matrix (GO:0062023), and extracellular exosome (GO:0070062) ( Table 3 ). These GO terms primarily reflect processes and components related to ECM organization and localization, which are known to play critical roles in the pathogenesis of IPF. Moreover, there were 4 shared KEGG pathways between 475 DEGs and 402 DEGs: complement and coagulation cascades (hsa04610), cytokine-cytokine receptor interaction (hsa04060), rheumatoid arthritis (hsa05323), and staphylococcus aureus infection (hsa05150). 402 up-regulated DEGs also exhibited 4 REACTOME pathways: Immune system (R-HSA-168256), extracellular matrix organization (R-HSA-1474244), molecules associated with elastic fibers (R-HSA-2129379) and elastic fiber formation (R-HSA-1566948) as shown in Table S4 . Top five significantly enriched GO terms among DEGs from our IPF fibroblast dataset and overlap DEGs of four studies. #: GSE180415 (Hanmandlu et al. 2022) and GSE17978 (Emblom-Callahan et al. 2010) and GSE1724 , 10921, 44723, 40839 (Plantier et al. 2016) and our study. Enrichment significance values correspond to FDR-adjusted p -values (FDR < 0.05). Abbreviations: GO: gene ontology; DEG: differentially expressed gene; IPF: idiopathic pulmonary fibrosis; GSE: gene expression omnibus series; FDR: false discovery rate; CC: cellular component. We compared the 475 DEGs identified in our study with those reported in three previously published gene expression data of IPF fibroblasts ( GSE180415 ( 9 ), GSE17978 ( 10 ) and GSE1724 , 10921, 44723, and 40839 ( 11 ). A total of 85 DEGs were identified as overlap genes among four studies ( Table 4 ). Table S5 showed the GO analysis result of 85 DEGs. There were the five GO terms and one KEGG pathway such as extracellular region (GO:0005576), extracellular space (GO:0005615), collagen-containing extracellular matrix (GO:0062023), external side of plasma membrane (GO:0009897), extracellular matrix (GO:0031012), extracellular matrix organization (R-HSA-1474244), and cytokine-cytokine receptor interaction (hsa04060). Notably, the result shared the top three GO terms of our study such as extracellular region (GO:0005576), extracellular space (GO:0005615), and collagen-containing extracellular matrix (GO:0062023) as top GO terms with other two ECM related GO terms ( Table 3 ). A list of 85 differentially expressed genes shared across four studies. A table shows the genes shared with between our study and three other studies: GSE180415 (Hanmandlu et al. 2022) and GSE17978 (Emblom-Callahan et al. 2010) and GSE1724 , 10921, 44723, 40839 (Plantier et al. 2016). Abbreviation: GSE: gene expression omnibus series. STIRNG analysis of 475 DEGs identified 410 nodes and 726 edges, where MCODE made 3 clusters (score≥4, Figure 2 ). 437 DEGs within 475 were employed for the GeneClip3 analysis based on the number of functional connections (degree) within the literature- and GO-based interaction network ( 39 ). The top ~3% of nodes by degree were defined as hubs ( Figure 3 ). Three genes, SERPING1 , NR4A1 and C3 , were concurrently shown in two analyses. Gene and protein network (8 node, 14 edges, score=4) of 85 shared DEGs in 4 studies indicated that two well-known IPF related genes, DPP4 and CXCL2 ( Figure 4 ). Protein-protein interaction network analysis and cluster of 475 differentially expressed genes by STRING and MCODE. (a-c) Cluster 1=15 nodes, 54 edges, score:7.7, cluster 2=14 nodes, 28 edges, score:4.3, cluster 3=11 nodes, 20 edges, score:4. Functional gene association network of 475 differentially expressed genes via GeneClip3. Circle: gene; square: transcription factor; triangle: enzyme; sky blue mark on shapes: already known as IPF related genes. Network analysis of 85 differentially expressed genes. (a) gene network and (b) protein network analysis of overlapping genes among 4 studies. Circle: gene; square: transcription factor; triangle: enzyme; sky blue mark on shapes: already known as IPF related genes.

Our large-scale RNA-seq analysis of IPF lung fibroblasts identified 475 DEGs, with a strong enrichment of GO terms related to ECM organization and localization. Among the top five GO terms were extracellular region (GO:0005576), extracellular space (GO:0005615), and collagen-containing extracellular matrix (GO:0062023). Comparative analysis with three previous IPF fibroblast transcriptomic datasets revealed that, although only 85 DEGs and 14 GO terms overlapped with any of the prior studies and no single gene was common to all ( Figure S3 ), these ECM-related GO terms consistently emerged across datasets ( Tables S4 - 6 ). This convergence, despite differences in sample origin, cohort size, and analytical pipelines, underscores the robustness and disease relevance of ECM remodeling as a central process in IPF fibroblasts. Among the enriched pathways, the cytokine–cytokine receptor interaction (hsa04060) pathway was consistently activated across all four datasets. This underscores the central role of cytokine signaling in fibroblast activation and intercellular communication in the fibrotic lung microenvironment. Pro-fibrotic cytokines such as IL-1β, IL-6, TNF-α, and TGF-β1 sustain myofibroblast activation and ECM deposition ( 40 ). The recurrent enrichment of this pathway suggests that fibroblasts may both respond to and contribute to the local cytokine milieu, promoting persistent fibrotic remodeling. In addition to cytokine signaling, three immune- and inflammation-associated pathways—complement and coagulation cascades (hsa04610), rheumatoid arthritis (hsa05323), and Staphylococcus aureus infection (hsa05150)—were also enriched. The complement pathway has been implicated in IPF progression, with elevated levels of C3, C5, and C5–C9 complexes correlating with disease severity in plasma and bronchoalveolar lavage fluid ( 41 ). The rheumatoid arthritis (RA) pathway enrichment highlights shared molecular mechanisms between IPF and RA-associated interstitial lung disease, both of which often present with usual interstitial pneumonia patterns ( 42 ). Triggianese et al proposed that complement activation, fibroblast activation, and pro-inflammatory cytokine signaling are common mechanisms in both diseases ( 43 ). Enrichment of the Staphylococcus aureus infection pathway is notable, as colonization has been associated with acute exacerbations and poor prognosis in IPF ( 44 ). Altered ECM structure and immune dysfunction may facilitate bacterial persistence, and fibroblast responses to microbial components may further exacerbate fibrosis. Network analysis identified SERPING1 , NR4A1 , and C3 as hub genes in STRING and GeneClip3-based networks, indicating their central roles in fibroblast regulatory circuits. SERPING1 has been reported as a serum biomarker in IPF ( 45 ); NR4A1 modulates TGF-β signaling and suppresses fibrosis in several organs ( 46 ); and C3 -deficient mice are partially protected from bleomycin-induced fibrosis ( 47 ). These genes were upregulated in primary IPF fibroblasts, supporting their relevance to IPF pathology. Analysis of the 85 shared DEGs revealed a subnetwork including CXCL2 , linked to IL1R1 , ICAM1 , and OAS2 . While IL1R1 and ICAM1 are established pro-fibrotic mediators ( 48 , 49 ), OAS2 , an interferon (IFN)-induced gene with antiviral function ( 50 ), has not been studied in IPF. Emerging evidence has implicated interferon signaling in IPF pathogenesis, with single-cell RNA-seq studies showing elevated IFN-γ activity in immune cells in fibrotic lungs ( 51 , 52 ). Importantly, these studies revealed that interferon-activated immune cells are often positioned near fibroblast-rich regions. This physical closeness makes it easier for immune-derived cytokines and paracrine signals to directly influence fibroblast behavior, including activation and extracellular matrix production. Such immune–fibrotic interplay suggests that OAS2-driven interferon responses could shape fibrogenic remodeling in IPF, highlighting a potential translational link between antiviral defense programs and fibrosis ( 53 ). Notably. of the 457 DEGs identified in our study, 399 were not found in any of the three datasets. Several of these genes — BMP5, TMEM176B, C7, GDF10, ADAMTS8, EDNRB, FGFR4 , and TNFSF15 — are known to participate in fibrosis in other organs such as the liver, kidney, and heart, suggesting conserved fibrogenic pathways across tissues. We also uncovered several novel genes not previously linked to IPF, including ADH1A, ADH1B, SCUBE1, NKD2, CFTR, ADAMTS19, HSD17B2, RASL12 , HMGCLL1, HBG1 and the long non-coding RNA RIPOR3-AS and LOC105375566. These genes may represent new regulators of fibroblast activation or matrix remodeling. ADH1A, ADH1B, SCUBE1 , NKD2 , CFTR, and ADAMTS19 are associated with liver, renal, cystic fibrosis, and skin fibrosis ( 28 , 29 , 33 , 34 , 36 , 37 ) but haven’t been studied in IPF. HSD17B2 exhibiting a well-established function in hormone metabolism ( 54 ), has recently been reported to be correlated with TGF-β/Snail mediated epithelial-mesenchymal transition in endometriosis ( 55 ). This suggests that in fibroblasts HSD17B2 may contribute EMT-like or mesenchymal gene activation under profibrotic signaling through hormonal regulation. The reported significance of HMGLL1 in lipid biosynthesis and energy metabolism ( 56 ) suggests a plausible link to IPF pathogenesis. By supplying cytosolic acetyl-CoA, HMGCLL1 may promote histone acetylation and profibrotic gene expression, while simultaneously fueling lipid biosynthesis that predisposes fibroblasts to ferroptotic lipid peroxidation such as 4-hydroxy-2-nonenal, malondialdehyde accumulation). As Xin Geng et al. reported that increasing hemoglobin-oxygen affinity alleviates lung fibrosis in a bleomycin-induced model ( 57 ), HBG1 may be functionally linked to IPF pathology. RASL12 , a member of Ras-related GTPase, participates in intracellular signaling involved in fibrosis ( 58 ). Given that RIPOR3-AS1 was recently discovered as senescence-induced lncRNA ( 59 ) and that fibroblast senescence together with metabolic reprogramming are hallmarks of IPF ( 60 , 61 ), its dysregulation may contribute to disease progression. According to the NCBI Gene database (Gene ID: 105375566), LOC105375566, currently annotated as an uncharacterized lncRNA, is located on chromosome 7q36 in close proximity to TMEM176B (NCBI Gene, Gene ID: 28959, accessed September 17 2025). Since TMEM176B has recently been reported to alleviate pulmonary fibrosis by inhibiting TGF-β/SMAD signaling ( 21 ), it is tempting to speculate that LOC105375566 might exert cis-regulatory effects on nearby genes. However, genomic proximity alone does not imply functional relatedness, and direct experimental evidence is still lacking. In addition, we observed dysregulation of genes involved in organ development ( GLDN ), wound healing ( DMKN ), and ion transport ( KCNJ2 , KCND1/2/3 , KCNE3, KCNG ) that have not been characterized in IPF. Similar changes in developmental regulators and potassium channels were noted in previous methylation and transcriptome studies ( 9 ), suggesting that fibroblasts may acquire aberrant features beyond matrix production, potentially due to epigenetic reprogramming. Several limitations must be considered. First, the use of stringent criteria (Bonferroni-adjusted p 2), may have excluded genes with modest but biologically relevant expression changes. Future studies employing more permissive thresholds or machine learning–based feature selection could uncover additional candidates. Second, although primary fibroblast cultures enable cell type–specific transcriptomic profiling, in vitro expansion may lead to phenotypic drift or overrepresentation of specific subpopulations, limiting fidelity to the in vivo microenvironment. Third, control fibroblasts were derived from histologically normal lung tissue adjacent to resected tumors. This pragmatic approach is widely adopted in human studies, but the possibility of subtle transcriptional alterations related to the tumor microenvironment cannot be completely excluded. Finally, our gene ontology analysis using the DAVID database depends on existing annotations and may overlook under-characterized pathways. Nevertheless, the consistent enrichment of ECM-related terms across our and previous datasets reinforces the biological validity of our findings. Future studies should prioritize functional validation of novel DEGs identified in our analysis. Gene silencing or overexpression experiments in fibroblasts could determine their roles in collagen deposition, myofibroblast differentiation, and cell migration. Single-cell RNA-seq may reveal whether these genes are enriched in specific fibroblast subsets. Integration with epigenomic and proteomic data may reveal upstream regulatory mechanisms and therapeutic targets. Clinically, correlation of gene expression with patient survival or treatment response could inform biomarker discovery. In summary, this transcriptomic analysis of primary lung fibroblasts provides a cell-type–specific view of IPF-associated transcriptional remodeling. Our findings confirm the central role of ECM dysregulation and identify numerous novel candidate genes that expand our understanding of fibroblast heterogeneity and activation in IPF. These data offer a valuable resource for future mechanistic and translational research aimed at improving IPF diagnosis and treatment.

Lung fibroblasts from patients with IPF were obtained from the biobank of Soonchunhyang University Hospital (Bucheon, Korea), following Institutional Review Board approval (IRB numbers: SCHCA-IRB-2018-10-034 and 201910-BR-058). Informed written consent was obtained from all participants. Subjects underwent clinical evaluation, including medical history review, chest X-ray, high-resolution computed tomography, and pulmonary function testing. Patients with evidence of collagen vascular diseases were excluded. The diagnosis of IPF was established according to the 2011 and 2018 ATS/ERS/JRS/ALAT guidelines ( 12 ). Primary lung fibroblasts were isolated from surgical biopsy specimens from 33 patients with histologically confirmed usual interstitial pneumonia, and from histologically normal lung tissue of 10 patients who underwent resection for stage I or II lung cancer. Fibroblast culture protocols followed previously published methods ( 13 ). Briefly, tissues were minced and cultured in 150 cm 2 flasks containing Dulbecco’s Modified Eagle Medium (Lonza), 10% fetal bovine serum (Thermo Fisher Scientific), 2 mmol/L glutamine, and 1% penicillin-streptomycin-amphotericin B (Lonza). Cultures were maintained at 37°C in a 5% CO 2 incubator and serially passaged to obtain a morphologically homogeneous population of adherent fibroblasts. Total RNA extraction and sequencing were performed by Macrogen Inc. (Seoul, Korea), following previously described protocols ( 14 ). Total RNA was extracted from lung fibroblasts using TRIzol reagent (Invitrogen). RNA concentration was measured using Quant-iT RiboGreen RNA Assay Kit (Invitrogen, # R11490 ), and RNA integrity was assessed with the Agilent TapeStation RNA ScreenTape (Agilent, #5067-5576). Only RNA samples with an RNA integrity number > 7.0 were used for library preparation. For each sample, 0.5 µg of total RNA was used to prepare a library using the Illumina TruSeq Stranded Total RNA Library Prep Gold Kit (#20020599). Ribosomal RNA was depleted, and the remaining RNA was fragmented by divalent cations under elevated temperature. First-strand cDNA synthesis was performed using SuperScript II reverse transcriptase (Invitrogen, #18064014) and random primers, followed by second-strand synthesis using DNA polymerase I, RNase H, and dUTP. The resulting cDNA fragments underwent end repair, A-tailing, and adapter ligation. Libraries were purified and PCR-amplified to yield the final cDNA libraries. Library quality and concentration were assessed using the KAPA Library Quantification Kit (KAPA BIOSYSTEMS, #KK4854) and Agilent TapeStation D1000 ScreenTape (#5067-5582). Libraries were sequenced on the Illumina NovaSeq platform (paired-end, 2 × 100 bp). Raw sequencing data quality was assessed with FastQC ( https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ ). Adapter sequences and low-quality reads (Phred score < Q20) were trimmed using FASTX_Trimmer ( http://hannonlab.cshl.edu/fastx_toolkit/ ) and BBMap ( https://sourceforge.net/projects/bbmap/ ). Clean reads were aligned to the human reference genome using TopHat ( 15 ). Gene expression was quantified using BEDTools ( 16 ), and Cufflinks ( 17 ), and normalized using EdgeR in R ( https://www.r-project.org ). Trimmed mean of M-values (TMM) was applied, and expression was presented as log2 counts per million. Trimmed mean of M-values (TMM) was applied, and expression was presented as log2 counts per million. DEGs were identified using the Exact test in EdgeR, and statistical significance was determined using Bonferroni-adjusted p-values (< 0.05). The RNA-seq data have been deposited in NCBI Gene Expression Omnibus under accession number GSE301181. Functional enrichment of DEGs was performed using DAVID (Database for Annotation, Visualization, and Integrated Discovery) v6.8 tool ( 18 ). GO terms were classified into Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) based on the GO database ( 19 ). Enrichment was also assessed using KEGG and REACTOME pathways integrated in DAVID. Statistical significance for enrichment was evaluated using the false discovery rate (FDR) correction, with FDR-adjusted p-values < 0.05 considered significant. Interaction networks were generated using GeneClip3 ( http://cismu.net/genclip3/analysis.php ) and STRING v2.2.0 (confidence cutoff = 0.4) via Cytoscape v3.10.3. Network clustering and module detection were performed using MCODE v2.0.3 in Cytoscape. Statistical analyses were conducted using SPSS (IBM). Data normality was assessed using the Shapiro–Wilk test. Continuous variables were expressed as mean ± standard deviation or median with interquartile range. Group comparisons were performed using independent t-tests or Mann–Whitney U tests as appropriate. Categorical variables were compared using Pearson’s Chi-square test. A two-tailed p-value < 0.05 was considered statistically significant.