Human bronchial epithelial cell line (BEAS-2B, CRL-3588) was purchased from the American Type Culture Collection, and 293 TN cells (LV900A-1) were obtained from the System Bioscience. Cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM, Gibco, United States) supplemented with 10% fetal bovine serum (Gibco, United States). All cells were routinely tested and free from microbial contamination.

Cells were seeded at a density of 1 × 10⁵ per well onto six-well culture plates and maintained in the DMEM in a humidified atmosphere of 5% CO₂ at 37°C. On the following day, BEAS-2B cells were treated with human recombinant TGF-β1 (2 μg/mL, rTGF-β1, PHG9204, Gibco, United States) at the concentration of 5 ng/mL, 10 ng/mL, 100 ng/mL, and 1,000 ng/mL for 48 h, respectively. Cells cultured in the absence of rTGF-β1 served as the control. After the treatment, cell extracts were collected for further analyses.

A human genome-wide CRISPR knockout pooled library (GeCKO v2) targeting 19,050 genes and 1,864 microRNAs (miRNAs) was purchased from the Addgene (#1000000048, United States). To produce the lentivirus for the GeCKO v2 sgRNA library, 3 × 10⁶ 293 TN cells were transduced with 9 μg of GeCKO v2 plasmids, 3 μg of pCMV-VSV-G, 3 μg of pMDLg, 3 μg of pRRE, and 3 μg of pRSV-Rev in a 150 mm dish containing 90 μL polyethylenimine. After 6 h, cells were washed with PBS, and fresh complete DMEM was added. The lentiviral supernatant was harvested after 48 h of incubation at 37°C and 5% CO₂. The supernatant was then centrifuged at 3,000 rpm for 15 min at 4°C and stored for further use.

Genome-Cas9-edited BEAS-2B (BEAS-2B-GeCKO) cells were generated as previously described (Wei et al., 2019). Briefly, BEAS-2B cells were transduced with GeCKO v2 sgRNA library at a low MOI (∼0.3). After incubating with puromycin for 14 days to eliminate untransduced cells, viable cells were collected to evaluate the transduction efficiency. (With regard to the knockout of individual candidate genes, corresponding gRNAs were introduced into the lentiCRISPR v2.1 vector (Supplementary Table S1), followed by the cell transduction.) Thereafter, 2 × 10⁷ BEAS-2B-GeCKO cells were treated with rTGF-β1 (10 ng/mL). We performed four consecutive rounds of TGF-β1 selection according to the EMT readout assay. After four rounds of rTGF-β1 treatment, viable cells were collected for subsequent analyses. Genomic DNA was extracted from TGF-β1-resistant BEAS-2B-GeCKO cells using the HiPure Tissue DNA Kit (Magen, China). The single guide RNA (sgRNA) sequences were amplified by polymerase chain reaction (PCR). The sequence of forward primer is 5′-TCTTGTGGAAAGGACGAAACACCG-3’; The sequence of reverse primer is 5′-ACCTTCTCTAGGCACCGGAT-3’. The gRNA sequences were amplified using 2 × PrimStar MasterMix (Takara, Japan) through the PCR. The conditions of PCR were as follows: 98°C for 30 s; 35 cycles (98°C for 10 s; 60°C for 10 s; 72°C for 50 s); 72°C for 1 min; 4°C. The distribution of sgRNA was assessed by Illumina next-generation sequencing (NGS).

Following the NGS, the top-scoring gene hits were identified by applying a threshold of total reads of sgRNAs>300 and the sgRNA diversity (the counts of sgRNAs targeting each gene) > 1. In addition, a dataset (GSE104908) was downloaded from the Gene Expression Omnibus database and differentially expressed genes (LogFC>2, q < 0.05) were chosen for analyses. A Venn diagram was constructed using the R software v4.3.3 to compare the overlap between the top-scoring genes and GSE104908. To investigate the expression levels of the identified key profibrotic regulators in patients with pulmonary fibrosis, the GSE40839 and GSE24206 datasets were also downloaded from the Gene Expression Omnibus database. Prior to the statistical analysis, the gene chip data were subjected to deduplication, log2 transformation, and quality control. The Mann-Whitney U Test was applied to analyze the expression levels of key profibrotic regulators between the pulmonary fibrosis patients and control groups. Data were visualized in violin plots using the R software v4.3.3. To identify significantly enriched signaling pathways, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were conducted using the R software v4.3.3. To characterize the biological process, cellular component, and molecular function, the Gene ontology (GO) analysis was performed using the WEB-based Gene Set Analysis Toolkit (https://www.webgestalt.org/). The STRING was used to visualize the protein-protein interaction (PPI) network.

First, the 1 mL of Trizol was added to cells (5 × 10⁶) and vortexed briefly. Then 200 μL of chloroform was added to the cell lysate, followed by incubation (5 min at room temperature) and centrifugation (10,000 rpm, 5 min at room temperature). After centrifugation, 500 μL isopropanol was added to the aqueous phase and mixed gently. Samples were then centrifuged, and the supernatant was discarded. The pellet was subsequently washed with 70% ethanol, and the supernatant was removed after centrifugation. The pellet was allowed to air dry and resuspended in 20 μL of RNase-free water. The purity and concentration of RNA was measured using the NanoDrop 2000 (Thermo Fisher Scientific; United States). The cDNA was reverse transcribed using a High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific; United States) according to the manufacturer’s instructions. The RT-qPCR was carried out using the PowerUp SYBR Green master mix (Applied Biosystems; United States) and target-specific primers (shown below) on a CFX96 Touch Real-Time PCR Detection System (BIO-RAD; United States). The β-actin was used as the reference gene for normalization. Relative mRNA expression of the targets was calculated as the fold change relative to PBS-treated cells using the delta-delta Ct method (2^−ΔΔCT). Sequences of primers for target genes are listed below:

β-actin FW: TGGCACCAGCACAATGAA.

β-actin RV: CTAAGTCATAGTCCGCCTAGAAGC.

Fibronectin FW: AGCAAGCCCGGTTGTTATGA.

Fibronectin RV: CCCACTCGGTAAGTGTTCCC.

Collagen-I FW: ATTCCAGTTCGAGTATGGCGG.

Collagen-I RV: CTTGAGGTTGCCAGTCTGCT.

The pLKO.1-shRNA constructs targeting COL20A1, COL27A1, and WNT11 were purchased from IGEbio (China). The stable knockdown of these three genes was conducted through the lentiviral transduction of shRNA expressing virus, followed by puromycin selection. The pLKO.1-scramble shRNA was used as the control. The sequences of shRNAs are the following:

shRNA-COL20A1: GCTACACCTTGCAGATCTTCG.

shRNA-COL27A1: GCAGTGGCTATTCGATCTTCC.

shRNA-WNT11: GCCTGTGAAGGACTCGGAACT.

shRNA-scrambled: CCTAAGGTTAAGTCGCCCTCG.

After the rTGF-β1 treatment, cells were washed twice with ice-cold PBS. Cells were then lysed by adding radioimmunoprecipitation assay buffer (Beyotime, China) with the protease inhibitor (Beyotime, China), followed by centrifugation (12,000 rpm, 15 min at 4°C). After discarding the supernatant, the loading buffer (EpiZyme, United States) was added, and cells were heated at 95°C for 10 min. Samples were electrophoresed in the 10% Sodium Dodecyl Sulfate polyacrylamide gels, and proteins were subsequently transferred to polyvinylidene difluoride membranes. Membranes were blocked with skimmed milk for 2 h at room temperature and incubated with primary antibodies overnight at 4°C. On the next day, polyvinylidene difluoride membranes were washed with TBST (100 mM NaCl, 10 mM Tris, 0.1% Tween 20), and blots were incubated with horseradish peroxidase-conjugated secondary antibodies for 2 h at room temperature. After washing the TBST, membranes were developed with ECL chemiluminescent substrate (Tanon, China). Data quantification was conducted using the Fiji software (ImageJ, United States). Relative protein levels were calculated after normalization to GAPDH.

Antibodies used are as follows: rabbit anti-Fibronectin (proteintech #15613-1-AP, 1:2000), mouse anti-Collagen Type I (proteintech #67288-1-Ig, 1:5000), and mouse anti-GAPDH (proteintech #60004-1-Ig, 1:10,000). HRP-conjugated secondary antibodies include goat anti-mouse (proteintech #RGAM001, 1:5000) and goat anti-rabbit (proteintech # RGAR001, 1:5000).

The cell viability was assessed using the cell counting kit-8 (CCK-8, Beyotime, China). Cells (3 × 10³ per well) were seeded onto 96-well plates in triplicate and incubated at 37% with 5% CO₂ for 24 h. Cells were then treated with 10 ng/mL of rTGF-β1 for 48 h. Subsequently, 10 μL of CCK-8 solution was added, followed by an additional 2 h of incubation. The absorbance at 450 nm was detected using a microplate reader (BioTeke, United States).

Cells (1 × 10⁴) were seeded onto the transwell inserts (Corning, United States) containing 200 μL of serum-free DMEM. The DMEM containing 20% fetal bovine serum (500 μL) was added into the lower chamber. After 48 h of rTGF-β1 treatment (10 ng/mL), cells that migrated/invaded into the lower surface were washed with sterile PBS, fixed with 4% paraformaldehyde (20 min), and stained with 0.1% crystal violet (15 min). Excess crystal violet and unfixed cells were then washed away with PBS, followed by the observation under a bright-field light microscopy (Olympus, Japan) at ×100 magnification. Relative cell migration was calculated as the ratio of migrated cells in experimental groups to the background migration in the control group.

Cells (1 × 10⁵ per well) were seeded onto six-well culture plates and cultured overnight. Cells were then treated with rTGF-β1 (10 ng/mL) for 48 h. Subsequently, cells were stained with 5 μL of Annexin V-APC (BioLegend, United States) and 10 μL of 7-amino-actinomycin D (eBioscience, United States) for 15 min at room temperature. The apoptotic cells were determined and quantified by flow cytometry, and data were analyzed using the FlowJo software (Tree Star, United States).

All data were expressed as mean ± SEM. Data analysis was performed using the GraphPad Prism 9.30 (GraphPad Software, United States). Data were analyzed using unpaired two-tailed t-test or Mann-Whitney U Test for comparison between two groups. The one-way analysis of variance (ANOVA) with Tukey’s post-test was applied for multiple comparisons. The p value < 0.05 was considered statistically significant.

To explore the successful induction of EMT in human bronchial epithelial cells and its potential association with ECM synthesis and deposition, we aimed to establish a TGF-β1-induced EMT in BEAS-2B cells and measured the level of Collagen-I and Fibronectin, both of which are essential components of ECM (Zhang et al., 2018). To identify the optimal concentration of rTGF-β1, BEAS-2B cells were treated with different doses of rTGF-β1 (0, 5, 10, 100, and 1,000 ng/mL) for 48 h, followed by detecting the mRNA expression of Collagen-I and Fibronectin. As shown in Figure 1A, the mRNA expression levels of Collagen-I and Fibronectin gradually increased with the rise of rTGF-β1 concentration, peaking at 10 ng/mL. This result demonstrates the presence of an ECM-generating phenotype in BEAS-2B cells and the successful induction of EMT. However, a significant decrease in the expression of Collagen-I and Fibronectin was observed at 100 ng/mL and 1,000 ng/mL (Figure 1A), possibly due to the detrimental effect of higher TGF-β1 concentrations on cell viability. We therefore confirmed the use of 10 ng/mL of rTGF-β1 for subsequent experiments. After four consecutive rounds of rTGF-β1 treatment, higher levels of Collagen-I and Fibronectin mRNA expression were found in BEAS-2B-GeCKO cells compared with control. This result confirms the successful induction of EMT in BEAS-2B-GeCKO cells, paving the way for the establishment of the GeCKO screening (Figure 1B).

To uncover new molecular regulators of the TGF-β1-induced EMT, BEAS-2B cells were transduced with GeCKO v2 sgRNA library to generate the BEAS-2B-GeCKO cells. The Figure 2A illustrates the schematic flow of EMT induction by treating BEAS-2B-GeCKO cells with rTGF-β1 and subsequent high-throughput sequencing. We analyzed two samples in the testing set using the NGS. For each sample, a minimum of 5.1 GB of data were generated. A total of 19,050 genes and 1,864 miRNAs were detected through the NGS, in line with the number of genes and miRNAs targeted by the GeCKO pooled library. Amplified sgRNA fragments were able to be detected in BEAS-2B-NC cells and BEAS-2B-GeCKO cells as compared to BEAS-2B cells, proving that the GeCKO v2 sgRNA library was successfully transduced into BEAS-2B cells (Figure 2B). By analyzing 14M Illumina NGS reads from 3.6×10⁷ BEAS-2B-GeCKO cells post the rTGF-β1 selection, we observed the enrichment of 70,255 sgRNAs, representing 57% of total 123,411 sgRNAs in the GeCKO library. In addition, we detected an enrichment of 20,368 genes, including miRNAs, with 828 or 886 genes represented by six independent sgRNAs in the library, 3,247 or 3,132 genes represented by five independent sgRNAs, 5,553 or 5,607 genes represented by four independent sgRNAs, 5,755 or 5,861 genes represented by three independent sgRNAs, and 3,546 or 3,436 genes represented by two independent sgRNAs (Figures 2C, D).

Among the enriched genes, we identified a total of 5,643 top-scoring genes by applying the threshold of total reads of sgRNAs>300 and sgRNA diversity>1. This result suggests the potential profibrotic properties of these top-scoring genes and demonstrates that our GeCKO screening can accurately screen essential and functional molecular regulators involved in the TGF-β1-induced EMT. To identify crucial molecular regulators for further study, we compared our GeCKO screening data with the differentially expressed genes from GSE104908 (Figure 2E). Through the analysis of Venn diagram, 76 top molecular regulators were confirmed (Figure 2E).

To explore the PPI network among the identified top molecular regulators, we conducted the STRING analysis. As illustrated in Figure 3A, the PPI network comprising 76 nodes and 53 edges was constructed by applying a threshold of p < 0.05 and interaction score>0.4 (Disconnected nodes were hidden). Moreover, we performed KEGG pathway analysis for the top molecular regulators, and results exhibited enrichment in the TGF-β signaling pathway as well as other pathways associated with ECM dysregulation and EMT (e.g., focal adhesion and ECM-receptor interaction) (Figure 3B). Notably, we observed the enrichment of INHBA and found that the INHBA is involved in the TGF-β signaling pathway (Figure 3B and Supplementary Table S2). This result is in line with prior studies demonstrating the crosstalk between INHBA, EMT initiation and TGF-β pathway activation (Yu et al., 2021). In addition, the enrichment of COL6A3 was found, and it participates in ECM-receptor interaction, focal adhesion, and the PI3K-Akt signaling pathway (Figure 3B and Supplementary Table S2), all of which are associated with the pathogenesis of pulmonary fibrosis. The WNT11 and Wnt signaling pathway were also found to be enriched (Supplementary Table S2). Furthermore, the top-ranked biological process, cellular component, and molecular function of those top molecular regulators were analyzed through the GO analysis. As shown in Figure 3C and Supplementary Table S3, other molecular regulators (e.g., COL20A1, COL27A1, COL6A3, SMAD7, and WNT11) are involved in ECM structural constituent, collagen trimer, ECM component, adherens junction, collagen-containing ECM, as well as ECM organization. (Detailed information of KEGG and GO analyses is provided in Supplementary Tables S2, S3) Altogether, as these enriched genes are closely associated with ECM accumulation and TGF-β pathway activation, we therefore deduced that they may play a central role in regulating the TGF-β1-induced EMT and pulmonary fibrosis.

By integrating results of the Venn diagram with findings from pathway and functional enrichment analysis, the candidate genes (COL6A3, MMP3, WNT11, SMAD7, COL27A1, WNT5B, SERPINE2, INHBA, COL20A1, and COL7A1) were finally confirmed. The sgRNA counts for each candidate gene are shown in Figure 4A.

To assess whether these candidate genes play a key role in regulating the TGF-β1- induced EMT, we generated 10 BEAS-2B cell lines through the specific knockout of every single candidate gene, followed by the rTGF-β1 treatment. As presented in Figure 4B, the knockout of candidate genes (COL20A1, COL6A3, COL7A1, SERPINE2, SMAD7, and WNT11) greatly decreased the mRNA expression levels of Collagen-I. Moreover, although higher mRNA levels of Fibronectin were detected after knocking out the COL7A1, SERPINE2, WNT5B, and SMAD7, the substantial reduction in Fibronectin mRNA was observed upon the knockout of COL20A1 and COL27A1 (Figure 4B). In addition, post the rTGF-β1 treatment, a significant reduction in cell viability was observed when knocking out the candidate genes (MMP3, COL7A1, SERPINE2, and WNT11) compared with the control (Figure 4C). In contrast, the knockout of candidate genes (COL6A3, WNT5B, and INHBA) in BEAS-2B cells promoted cell proliferation (Figure 4C). Furthermore, Figures 4D, E show that the relative migration/invasion abilities of BEAS-2B cells following rTGF-β1 treatment were inhibited upon the knockout of MMP3, INHBA, and WNT11, whereas the knockout of WNT5B and SMAD7 facilitated the migration/invasion of BEAS-2B cells.

Based on the ten BEAS-2B cell lines with specific knockout of every single candidate gene, we performed the flow cytometry to assess whether the loss-of-function of candidate genes affected the cell apoptosis. As shown in Figures 5A, B, only the knockout of COL20A1, COL27A1, and COL7A1 greatly enhanced the rate of apoptosis in BEAS-2B cells after rTGF-β1 treatment, whereas knocking out the COL6A3 led to a decrease in the rate of apoptosis. Moreover, after suppressing the expression of COL20A1, COL27A1, and WNT11, the relative protein levels of Fibronectin and Collagen-I (key fibrosis-regulating proteins) were diminished (Figure 5C). These results prove that the COL20A1, COL27A1, and WNT11 may serve as key profibrotic regulators and contribute to the ECM deposition and TGF-β1-induced EMT.

To further consolidate the findings of our study, we analyzed the expression levels of COL20A1, COL27A1, and WNT11 in patients with pulmonary fibrosis. As presented in Figures 6A, B, although there were no differences in the expression of COL20A1 between IPF patients and the control in GSE24206 dataset, a substantial increase in COL27A1 expression was observed in IPF patients. In addition, the WNT11 expression was significantly higher in patients with pulmonary fibrosis as compared with the control group in GSE40839 dataset (Figure 6C).

Altogether, all of our findings demonstrate that the COL20A1, COL27A1, and WNT11 serve as key profibrotic regulators and may play a crucial role in regulating the TGF-β1-induced EMT and the pathogenesis of pulmonary fibrosis.