Animal modeling and sample collection

C57BL/6 mice (8 weeks, 18∼20 g, male) were obtained from Shanghai SLAC Laboratory Animal Co., Ltd (Shanghai, China) and randomly divided into four groups (15 for control, 11 for BLM, 11 for Agomir-NC, and 11 for Agomir). Animals in the BLM group were intratracheally administered a 1 mg/kg dose of bleomycin (BLM) (MedChemExpress LLC, China) dissolved in 50 µL sterile saline. The control group received 50 µL sterile saline intratracheally in the same manner. Animals in the Agomir group received a tail vein injection of Agomir-125b-5p or Agomir-NC (Shanghai GenePharma Co., Ltd) at a dose of 4 mg/kg dissolved in 100 µL sterile saline on day 7 of BLM treatment. The mice were euthanized on day 14 after BLM administration, after blood oxygen monitoring and pulmonary function testing [38]. Lung tissues were collected for fixation, embedding, and sectioning. Subsequently, H&E, Masson staining, and immunohistochemical assays were performed. Fresh lung tissues were also required for protein extraction and Western blot (WB) assays. All animal procedures were approved by the Research Ethics Committee of the Second Affiliated Hospital of Fujian Medical University (2019 − 103).

Lung histology, trichromatic staining, immunohistochemistry (IHC) and tyramine signal amplification (TSA)

Mouse lung tissue was soaked in a 10% neutral buffered formalin solution overnight, dehydrated according to standard procedures, degreased in xylene embedded in paraffin, and cut into 5 μm- thick slices. The slices were then stained with hematoxylin and eosin (H&E). The fractional area of pathological alterations in the lung parenchyma was calculated using the Image J software (v1.8.0). To observe the degree of lung fibrosis, the slices were stained using the Masson Trichrome Staining Kit (Cat # 77148, Guangzhou Wexis Biotech Ltd.), and the volume fraction (%) of collagen fibers was analyzed using Image J. Fibrosis markers (FN1 and α-SMA) and BAK1 (12105 S, CST) expression were evaluated using IHC [39]. After dewaxing, the slices were incubated in a 0.01 mol/L citric acid buffer, pH 6.0,  in a microwave oven for 15 min for antigen repair. The slices were then cooled and incubated in 3% H2O2 for 10 min to inactivate endogenous peroxidase, then blocked with a non-immune goat serum at room temperature for 10 min. Next, the primary antibody was added overnight at 4 ℃. Thereafter, on adding the MaxVisionTM (KIT-5005 MaxVisionTM HRP-Polymer anti-rabbit IHC Kit) reagent, the tissue slices were incubated at room temperature for 15 min. Finally, the freshly prepared DAB color developing solution was added. After 5 min, the stained cells were mounted and viewed with a light microscope.

For TSA experiments on paraffin sections of mouse lungs, the kit (G1236, Servicebio® TSAPLus, Wuhan servicebio Co., Ltd, China) was used, and before incubation of the antibodies, the tissue slices were autoclaved with sodium citrate for 3 min, endogenous peroxidase was eliminated with 3% H2O2, and blocked with goat serum for half an hour. The staining schemes were iF488-Tyramide-SFTPC (AB3786, Millipore), iF555-Tyramide-BAK1 (12105 S, CST) and iF647-Tyramide-CC10 (ab40873, Abcam). It is worth noting that a sealing step with goat serum is required after each stripping of antibodies before subsequent antibody incubation can take place.

MicroRNA sequencing

Fresh mouse lung tissue was cut into 5 mm small tissue blocks and immediately frozen in liquid nitrogen. The total RNA was extracted using the TRIzol reagent (Invitrogen, 15596018, USA). The microRNA sequencing was conducted by Haplox (Jiang Xi, China). The concentration, quality, and integrity of the total RNA were measured using a Qubit RNA HS Assay Kit in the Qubit ® 3.0 Fluorometer (Life Technologies, USA) and an Agilent 2100 Bioanalyser (Agilent Technologies, USA). Next, for each qualified sample, a 3 µg RNA was used to establish the small RNA library. Following the manufacturer’s instructions, NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB, USA) and index codes were assigned to attribute sequences to each sample, and clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE ClusterKit v3-cBot-HS (Illumina, USA). After cluster generation, the prepared library was sequenced on an Illumina platform. Raw sequencing data were deposited on the NCBI SRA database (accession number PRJNA1165890).

Cell culture and cell transfection

The alveolar epithelial type II cell line A549 was obtained from the National Collection of Authenticated Cell Cultures (serial number SCSP-503, Shanghai, China). Beas-2B, the human bronchial epithelial cell line, was purchased from IMMOCELL (Xiamen, Fujian, China). Human embryonic kidney cell line 293T was purchased from the Shanghai Institute of Cell Biology. All cells were cultured in DMEM (Gibco) supplemented with 10% FBS. The cells were maintained in a 37 ℃, 5% CO2 incubator in a medium of DMEM or 1640 supplemented with 10% FBS (A5669701, Gibco), 100 U/mL penicillin G, and 100 µg/mL streptomycin (Solarbio, P1400). In the cell experiments section, three biological replicates were utilized for the experiments, and for each cell plating, six technical replicates were included to ensure data reliability.

MiRNA mimic and inhibitor and plasmid transfection

The miR-125b-5p mimic and inhibitor were purchased from GenePharma Co., Ltd. (Shanghai, China). Sangon Biotech (Shanghai) Co., Ltd. constructed BAK1 overexpression plasmids by synthesizing and cloning the sequences into the pcDNA3.1(+) vector. Transfection of A549 and Beas-2B cells was performed under optimal growth conditions at 50–70% confluence using Lipofectamine® 3000 (Cat # L3000015; Invitrogen Thermo Fisher Scientific, Inc.). A 150 pM transfection of miR-125b-5p mimics and inhibitors was performed at 37 ℃ for 48 h. The cells were then harvested for subsequent experimentation. The BAK1 3’ UTR fragment containing the miR-125b-5p binding site was inserted into the pmirGLO vector (Promega Corporation) and named pmirGLO-BAK1-WT. Meanwhile, the mutation vector was constructed and named pmirGLO-BAK1-MUT. The BAK1 3’ UTR WT and MUT sequences are listed in Table 1.

MiRNA target prediction

To accurately identify the core target of microRNA has-miR-125b-5p, we used the following five databases: ENCORI, miRDB, RNAInter, TargetScan, and miRWalk. The Encyclopedia of RNA Interactomes (ENCORI, formerly known as starBase, https://starbase.sysu.edu.cn/) is a public platform for RNA-RNA and RNA-protein interactions identified from CLIP-seq and other high-throughput sequencing data [40]. The miRNA regulators were selected using seven target-predicting programs (PITA, RNA22, miRmap, DIANA-microT, miRanda, PicTar, and TargetScan) of the miRNA-Target “miRNA-mRNA” database in ENCORI. MiRDB (http://www.mirdb.org/) uses a bioinformatics tool to predict miRNA-target interactions, including 3.5 million predicted targets regulated by 7000 miRNAs [41]. RNAInter v4.0 (http://www.rnainter.org/), a complete resource of RNA interactome data from the literature and other databases, contains over 47 million RNA-related interactions [42]. For further analyses, we used the mRNA–miRNA interactions targeted by > 4 programs, score > 90, and score > 0.5 in ENCORI, miRDB, and RNAInter, respectively. MiRNA–mRNA interactions were additionally retrieved from TargetScan 8.0 (https://www.targetscan.org/vert_80/) [43]. MiRWalk (http://mirwalk.umm.uni-heidelberg.de/) holds data predicted by a machine learning algorithm, including experimentally validated miRNA-target interactions [44]. MiRNAs interacting with mRNAs commonly determined by > 5 databases (ENCORI, miRDB, RNAInter, TargetScan, miRWalk) were considered as mRNA regulators. Venogram was used to visualize the targets predicted by different databases.

Luciferase reporter assays

The cells were seeded into 96-well plates (2 × 104 cells/well). Co-transfection of reporter plasmids pmirGLO-BAK1-WT and pmirGLO-BAK1-MUT (200 ng each) with miR-125b-5p mimic or mimic-NC was performed using Lipofectamine 3000 at 37 ℃ for 48 h according to the manufacturer’s protocol. Firefly luciferase activity was measured according to the manufacturer’s instructions (Promega, Beijing, China) and normalized to Renilla luciferase activity.

Quantitative RT-PCR (qRT-PCR)

Total RNA was extracted from mouse tissues and cell lines using the TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.). 500ng of total RNA were reversely transcribed with specific primers using the TAKARA TB Green® Premix Ex Taq™ (Code No. RR420A) for miRNA and mRNA. The AceQ® qPCR SYBR® Green Master Mix kit (Code No. Q131, Vazyme Biotech Co., Ltd., Nanjing, China) was used for cDNA amplification that was conducted as follows: 94 ℃ for 5 min, then 35 cycles of 94 ℃ for 1 min, 60 ℃ for 30 s, and 72 ℃ for 1 min, followed by a final cycle of 72 ℃ for 10 min. The qRT-PCR was performed and analyzed in an ABI QuantStudio5 Q5 (Applied Biosystems, USA). The relative expression of miR-125b-5p was quantified using the 2−ΔΔCT method. U6 was selected as an internal control for miRNA. GAPDH or β-actin was used as an internal control for mRNA. The sequences of the miRNA mimics (Shanghai GenePharma Co., Ltd.) are listed in Table 1. The type of negative control used was non‑targeting. The specific primers (Shanghai Sangon Biotech Co., Ltd.) are listed in Table 2.

Differentially expressed gene analysis

Differentially expressed genes (DEGs) between the IPF (n = 30) and health control (n = 49) groups from RNA-seq dataset GSE124685 were determined using the limma software (v3.5.1) [45]. The genes expressed in less than 75% of samples were filtered. DEGs were defined as genes with p-value < 0.05 and |Fold Change| > 2. Expression of the target DEG was visualized using the ggboxplot function in the ggpubr R package.

Table 1 Sequences used for miR-125b-5p mimic, inhibitor, and BAK1 3’ UTR WT/ MUTTable 2 Primers used for reverse transcription-quantitative PCRWestern blot (WB)

A solution of RIPA: PMSF (Solarbio, Beijing, China) was used to extract total protein. Protein concentrations were measured using the BCA Assay Kit (Solarbio, Beijing, China). After separating on 11% SDS-PAGE, the proteins were transferred onto a PVDF membrane. Next, the membranes were blocked with 5% nonfat milk for 2 h, which was followed by incubation with primary antibodies (Vimentin, 1:1000; FN1, 1:100; E-cadherin, 1:1000; N-cadherin, 1:1000; α-SMA,1:20000; 1:1000; BAK1, 1:1000; β-actin, 1:1000) at 4 ℃ overnight. The used antibodies are listed in the Supplementary Table S1. Next, the primary antibodies were washed with the TBST buffer, and the membranes were incubated with the secondary antibodies for 1 h. The chemiluminescent signal analysis was performed using an imaging system (Invitrogen™, iBright™ CL1500 Imaging System) and quantified by the Image J software (v1.8.0).

Measurement of blood oxygen, pulmonary inspiratory resistance, and vital capacity

The mice were deeply anesthetized with sodium pentobarbital (40 mg/mL) and tested for respiratory health using the Small Animal Vital Signs Monitor, MouseOX® Small Animal Vital Signs Monitor (STRR, USA). Then, the animals were tracheostomized and placed into the body cavity of the system. The trachea was connected to the system tubing and ventilated for pulmonary function testing. Respiratory function indices, including relative pulmonary inspiratory resistance (RI) and pulmonary capacity (FVC), were measured using the Buxco Pulmonary function test system (DSI Buxco, PFT Controller, St. Paul, USA) [38]. The mice that died from respiratory failure or before the end of the experiment were not analyzed.

Analysis of BAK1 expression in different cell types from single-cell transcriptomic data of IPF clinical samples

The scRNA-seq data were downloaded from the NCBI GEO database under accession numbers GSE1437061 (IPF: HC = 7:5) and GSE1329152 (IPF: HC = 3:3), and included 43,691 cells from 10 IPF patient samples and 8 healthy controls. Figure S1 shows the range of number of genes detected in each sample, the depth of sequencing, and the percentage of mitochondrial, ribosomal, and erythroid content. We excluded cells with less than 500 expressed genes and more than 40,000 UMI. In addition, we excluded cells with more than 20% mitochondrial content, > 50% ribosomal content, and > 2% erythrocyte content. After these rigorous quality control screens to remove low-quality cells, we obtained 40,153 cells for downstream analysis. Based on the highly consistent results between patients and the poor results of using the Harmony algorithm to remove heterogeneity between individual patient IPF cells, we chose the batch effect correction algorithm.

The genes were normalized by the functions Normalize Data and Scale Data. Variable genes were then selected and subjected to principal component analysis using the FindVariableFeatures function. The most appropriate number of PCs (PCs = 30) was determined using Dim Plot and Elbow Plot (Figure S2). The first 30 principal components and Resolution 1.0 were used with the Find Clusters function to generate cell clusters. Clusters were manually annotated by publicly known biomarkers for each cell type (Figure S3). On this basis, we demonstrated the expression values of BAK1 gene in different cell populations in the IPF and HC groups.

Statistical analyses

The data were expressed as means ± SEM and analyzed using the GraphPad Prism software (v8.0). One-way ANOVA was used for multiple group comparisons. Student’s t-test was used only for two-group comparisons. A p-value < 0.05 was considered statistically significant.