Based on the background, we conducted the first systematic study on the effects of long-term low dose-rate radiation, using human bronchial epithelial cells (BEAS-2B) as a model to establish a long-term low dose-rate radiation-induced cell model. Our findings reveal that, compared to an equivalent total dose of single acute radiation, long-term low dose-rate radiation does not increase the incidence of lung cancer but significantly enhances the malignancy of cells. Furthermore, scRNA-seq analysis results provide emerging yet convincing systems biology evidence at both the molecular and cellular levels. Through this radiation-induced tumorigenesis model, we have, for the first time, uncovered the potential biological role of the ANGPTL4-SDC4 ligand-receptor pair in early-stage lung cancer. This study not only provides the first evidence and mechanistic explanation that long-term low dose-rate radiation exposure leads to increased cellular malignancy but also systematically delineates the dynamic process of malignant transformation, offering significant theoretical support and fresh perspectives for the fields of tumor biology and radiation biology.

The human bronchial epithelial cell line BEAS-2B was purchased from ATCC (CRL-9609) and maintained in DMEM (Gibco, Grand Island, NY, USA) containing 10% fetal bovine serum (FBS) (Gibco Invitrogen), penicillin (100 U/mL) and streptomycin (100 µg/mL) in a fully humidified incubator with 5% CO at 37°C. Cell lines were identified by STR profiling, and were tested mycoplasma contamination.

The α-irradiator used in this study comprises a Am source, a rotating radiation source holder, a sample holder, and other necessary components [15]. The irradiation source holder and the sample holder were 10 mm apart and parallel to each other. The Am source was positioned face down on the irradiation source holder and emitted α-particles at a dose rate of 0.14 Gy/min. Coverslips with BEAS-2B cells were placed on the sample holder. The Am source was rotated over the coverslip to irradiate the cells. The rotation of the Am source and the irradiation time were controlled automatically and accurately by a computer. The BEAS-2B cells at 70% confluence were irradiated either with a single dose of 0.2 Gy, 0.4 Gy, or 0.5 Gy for one time or irradiated with 0.02 Gy (9 s of irradiation time) once every three days for 10, 20, or 25 exposures, respectively (Fig. S1). For single exposures of 0.2 Gy (irradiated for 86 s), 0.4 Gy (171 s), and 0.5 Gy (214 s), the slides with cells were immersed in a culture medium once every minute in order to maintain humidity, and cell subculture was consistent with the long-term low dose-rate exposure group. For long-term low dose-rate exposures, all the irradiated cells were placed in a new dish for further culture, and all of them were planted on a new cell slide the next day, with irradiation repeated two days later. The sham control was similarly treated without irradiation. The seven groups were designated as following, Ctrl, sham control; R1_0.2, single exposure to 0.2 Gy; R1_0.4, single exposure to 0.4 Gy; R1_0.5, single exposure to 0.5 Gy; R10_0.2, 10 exposures of 0.02 Gy each time in 30 days; R20_0.4, 20 exposures of 0.02 Gy each time in 60 days; R25_0.5, 25 exposures of 0.02 Gy each time in 75 days. After all exposures were completed, the cells were continuously subcultured and harvested on 10 passage (Day 30), 30 passage (Day 90), and 50 passage (Day 150) for further analysis.

Mice were purchased from Shanghai Experimental Animal Center of the Chinese Academy of Sciences (Shanghai, China). Cells were harvested at a concentration of 4×10 cells/mL and mixed with the same volume of Matrigel. 0.1 mL of each mix was subcutaneously injected into either side of the posterior flank of male NOD/SCID mice (4-5 weeks old). Injection points in each mouse are shown in Fig. 1a. Each mouse had two injection points: Ctrl, R1_0.2, R10_0.2, R1_0.4, R20_0.4, R1_0.5, and R25_0.5 cells were injected into the hind leg armpits. A total of 93 mice were used in the study, random grouping was conducted. Four months after injection, the mice were sacrificed, and the subcutaneous growth of each tumor was examined. First, the tumor-like nodules were weighed individually. Then, partial tumor-like nodules from each irradiated sample were immersed in liquid nitrogen to preserve RNA. Next, the size of the remaining tumor-like nodules is documented through photographs. Tumor-like nodules were excised and subjected to pathological analysis and gene expression assays. For immunohistochemical staining, histological sections of formalin-fixed and paraffin-embedded tumor-like nodules were submitted to routine H&E staining and immunohistological assay of cancer-related proteins, including CAM5.2, CK7, and Ki67, etc. After deparaffinization and antigen retrieval using an autoclave oven, sections (5 μm) were incubated at 4°C overnight with CAM5.2, CK7, or Ki67, etc Rabbit Polyclonal Antibody (Abcam, UK). Antigen-antibody complexes were detected using the cobalt-3, 3'-diaminobenzidine reaction. The nodules from at least three mice were resected and processed for immunohistochemical staining. Positive expressions of CAM5.2, CK7 and Ki67, etc were shown in brown. For the microPET assay, F-FDG was synthesized by Jiangsu Huayi Technology Co., Ltd. Five NOD/SCID mice were fasted for 8 h in advance and anesthetized with isoflurane, the heating pads were used to keep the nude mice warm; after the animal's righting reflex disappeared, the imaging agent was injected into the tail vein, each 100 ± 20 µCi, and the volume was 0.1 mL. After the injection, the animals were kept under anesthesia for 1 h, fixed to the scanning bed, subjected to a PET static scan for 10 min, and then to a CT scan. PET scans and image analysis were performed using an Inveon microPET (Siemens Healthineers, Erlangen, Germany). The mouse experiments were carried out in strict accordance with the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health, as well as reviewed and approved by the Soochow University Institutional Animal Care and Use Committee (SUDA20230619A04).

After irradiation and subculture for up to the 10, 30 and 50 passages, cell proliferation was determined using the Cell Counting Assay Kit-8 (Dojindo Laboratories, Kumamoto. Japan). The BEAS-2B cells were seeded into 96-well plates (5×10 cells/well). After 24 h, 10 μL of CCK-8 was added to each well. Two hours post-incubation, the absorbance at 450 nm was measured using a spectrophotometer (Bio-Rad Laboratories, Inc.). The experiment was independently repeated three times. Commercial transwell chambers (24-well plate, 8 μm pore size, BD Biosciences, MA, USA) were used to assess cellular invasion ability. The upper surface of the chambers was coated with 10 μL of Matrigel (Corning, MA, USA). The cells, suspended in serum free medium, were seeded at a density of 5 × 10 cells/mL/well in the upper chamber, while the lower chamber contained 10% FBS-DMEM. After 48 h of incubation, non-invading cells were scraped from the upper surface with a cotton swab, and the invading cells on the bottom surface were stained with crystal violet. Stained cells were extracted with 0.5 mL of extraction solution (50:49:1 ethanol: distilled water: 1 M HCl) per well for 10 min under shaking. The optical density of each well was measured at 540 nm.

All cells were adjusted to a density of 5 × 10 cells/mL cell culture medium, and 200 μL of cell suspension was plated (i.e., 10,000 cells) into the wells of the sterile 96-well U-bottom cell culture plate. The plate was centrifuged at a speed of 500 rpm for 5 min and then placed in the incubator for seven days at 37°C and 5% CO. Afterward, 150 μL of the medium from each well was removed, and 25 μL of Matrigel was added into the remaining 50 μL of medium. After a short mixing with a pipette, the plate was left on ice for 5 min and then incubated for 30 min at 37°C before gently adding another 100 μL of culture medium into each well. Images were taken on days 1, 3, 5, and 7.

Initially, a quality inspection of the cells was performed for sequencing. Single-cell suspensions were mixed with an equal volume of 0.4% trypan blue staining solution, and live cell concentration was adjusted to 1,000-2,000 cells/μL using the Counters II Automated Cell Counter. Subsequently, cDNA fragments were labeled using 10X technology. Gel beads with barcode information were combined with cell and enzyme mixtures, separated by oil in the reservoir to form GEMs (Gel Beads-In-Emulsions). The gel beads dissolved, releasing the captured barcode-containing sequence. These cDNA fragments were then reverse transcribed and labeled. After breaking the gel beads and oil droplets, the cDNA served as the template for PCR amplification. The GEM products were pooled to construct a standard sequencing library. The library was established by cleaving cDNA into 200-300 bp fragments, followed by standard second-generation sequencing library construction steps, including end repair, A-tail addition, and sequencing connectors (P5, P7) and sample indices. PCR amplification was performed to obtain the DNA library, which was sequenced using the PE150 mode on the Illumina platform for high-throughput sequencing.

For data analysis, alignment (GRCh38), filtering, barcode counting, and UMI counting were performed using the Cell Ranger count module to generate the filtered gene-cell matrix. The matrix was then loaded into R software, use DoubletFinder (v2.0.4) to calculate the multi cell filtering threshold for each sample [16], and perform multi cell filtering in sequence. Then use the Seurat package (v4.2.1) to filter out low-quality cells or doublets (<200 or >7500 genes/cell, >10% mitochondria genes, >50,000 transcripts/cell) and genes (<3 cells/gene) [17]. Seven samples were integrated using data anchors with the FindIntegrationAnchors and IntegrateData functions. Gene expression levels were normalized (LogNormalize) using the NormalizeData function. A total of 2000 highly variable genes were obtained for PCA dimension reduction. Uniform manifold approximation and projection (UMAP) was performed for cell visualization. Except for the control sample, six samples were exposed to a total dose of 0.5 Gy. R1-10, 10 passage after one acute exposure; R1-30, 30 passage after one acute exposure; R1-50, 50 passage after one acute exposure; R25-10, 10 passage after low dose-rate exposure; R25-30, 30 passage after low dose-rate exposure; R25-50, 50 passage after low dose-rate exposure. The log2FC between two groups (R1-10, R1-30, R1-50, R25-10, R25-30, R25-50 versus sham control) was manually calculated using the Seurat package's FindMarkers function. The significance of the difference was determined using a two-sided Student's t test with Bonferroni correction. Signature genes were required to be expressed in >25% of cells within either of the two cell groups (marked as PCT, percentage of cells). Genes were selected as signatures based on the statistical threshold (logFC  >  0.25, two-sided Student's t test p-value < 0.01, and adjusted p-value (Bonferroni) < 0.01). Determine significantly enriched KEGG pathways of DEGs using the clusterProfiler package (v3.8.1) [18]. Transcription factor activity was analyzed using the DoRothEA package (v1.7.3) [19], with all parameters set to default values. Dotplots were generated using the DotPlot function of the Seurat package.

Somatic large-scale chromosomal copy number variation (CNV) score of exposed sample cell was calculated using the inferCNV package (v1.14.0) [20]. A raw counts matrix, annotation file, and gene/chromosome position file were prepared according to data requirements. Sham control cells without irradiation were selected as reference normal cells. The default parameters were applied. To quantify the CNV level of each cell, the CNV score of each gene on each observation set cell is subtracted by 1 and the quantitative sum is calculated, which is used as the CNV level of the cell. Next, classify and compare the CNV level of cells based on their sample and cluster attributes.

The GSEA function of the clusterProfiler package for was used for gene set enrichment analysis (GSEA) to calculate the biological functional differences between high_CNV and low_CNV cells. The input data is a gene list sorted in descending fold change order between two types of cells, calculated by FindMarkers function of Seurat package. Pathway analyses were predominantly performed on the Hallmark pathways and KEGG pathways in the MsigDB database [21]. For single-sample gene set enrichment analysis (ssGSEA), we use the gsva function of the GSVA package (v1.46.0) [22], set the kcdf parameter to 'Poisson', and performed statistical analysis on the results.

The cells from four samples, R1-30, R1-50, R25-30, and R25-50, were labeled with fluorescently conjugated CD44 (397503, BioLegend) and CD142 (365205, BioLegend) antibodies, respectively. An SH800S multi-application cell sorter (SONY Biotechnology) was used to separate viable high_CNV cells (CD44CD142) and low_CNV cells (CD44CD142).

Cells were harvested in TRIzol reagent (Invitrogen, CA, USA). Total RNAs were reverse-transcribed using the PrimeScript RT Reagent Kit (Takara, Kusatsu, Shiga, Japan). PCR analysis was performed using PowerUp™ SYBR Green Master Mix (Life Technologies, Grand Island, NY, USA), and amplified PCR products were quantified and normalized with GAPDH. PCR amplification was carried out using a Life Technologies system (Vii7A, NY, USA) and initiated by at 95°C for 2 min, followed by 40 cycles of 95°C for 30 s, 62°C for 40 s, and a final extension at 72°C for 10 min. Data were analyzed by C(t) value comparison method and normalized to control expression in each sample. The primer sets were listed in Supplemental Table S5.

The cells cultured on cover slides were fixed in 4% paraformaldehyde for 10 min, and washed with PBS for 30 min, treated with 0.5% Triton for 10 min, blocked with 5% skim milk for 1 h, and then, stained with anti-EpCAM (1:500, CST), anti-Vimentin (1:500, CST) or anti-Fibronectin antibody (1:500, BD) for 2 h. After washed three times with PBST, the cells were incubated with secondary antibodies for 1.5 h, and cell nuclei were counterstained with 4',6-diamidino-2-phenylindole (Invitrogen). Phase contrast images were captured using an Olympus IX70 microscope (Olympus). More than five fields of view from at least three independent experiments were randomly chosen.

Monocle 2 (v2.24.0) was utilized to infer cellular trajectory of cell evolution following two different exposure modes with the assumption that one-dimensional 'time' can describe the high-dimensional expression values, so called pseudotime analysis of single cells [23]. The general pipeline and default parameters were implemented for high_CNV cells in both exposure groups. The cell trajectory was visualized as tree structure in two-dimension plot after log normalization and DDRTree reduction dimension. Subsequently, heatmap was used to display the genes that significantly changed with the increase of pseudotime (cell passage), followed by unsupervised clustering of these genes. KEGG enrichment analysis was conducted on gene clusters showing increased expression levels.

CellChat was used for cell-cell communication calculation and analysis [24]. The CellChat results are divided into cell-cell communication count and cell-cell communication weight, used to reveal both incoming communication patterns of target cells and outgoing communication patterns of the source cells. We calculated the cell-cell communication information for 6 exposure samples and statistically compared their communication counts and weights. To explore the molecular processes of cell evolution after different exposures, we used a systems biology hierarchical analysis model of cell clusters to ligand-receptor pathways to ligand-receptor pairs. Specifically, at the cell cluster level, we used heatmap to illustrate the difference in the number/weight of cell-cell communication between two samples. For clusters with significant changes, we used bubble plot and heatmap to demonstrate their specific changes in ligand-receptor pathway signaling. Next, we identified representative ligand-receptor pair based on their contributions to the pathway.

Pre-cooled cell lysis buffer for Western blotting and immunoprecipitation (Beyotime, P0013) was added to a 60 mm culture dish of R25-50 sample. The lysates were then transferred into 1.5 mL EP tubes and gently shaken at 4°C for 10 min. Protein A/G Agarose (Beyotime, Shanghai, China, P2055) was used to incubate with anti-ANGPTL4 (Proteintech, 18374-1-AP) at 4°C for 8 h on a rotating incubator. The immunocomplexes were washed four times with cell lysis buffer for western blotting and immunoprecipitation and then incubated with cell lysates on a rotating incubator for another 8 h. After four additional washes, the immunocomplexes were dilute with loaded buffer. Finally, the presence of SDC4 protein in the immunocomplexes was detected using the immunoblotting method.

The protein structures were obtained from AlphaFold Protein Structure Database [25]. Molecular docking simulations were conducted using the online software GRAMM [26]. Initially, the crystal structure of the target receptor protein was prepared, by removing water and ligand molecules. Subsequently, the three-dimensional structures of multiple ligand molecules were prepared and optimized for conformation. Docking simulations were then performed using the software to explore the most probable ligand-receptor binding conformations. Finally, the binding affinity of the ligand-receptor complex was evaluated based on docking energy and binding sites using Prodigy [27]. The docking structure was visualized by PyMol (v2.2.0).

Multiplexed immunohistochemistry staining of ANGPTL4 (Rabbit, 1:2000, Proteintech) and SDC4 (Rabbit, 1:2000, Proteintech) was conducted according to the manufacturer's instructions (AiFang biological, 6-Color Multiple Fluorescence Kit, China). The images were captured using the multispectral microscope. More than 5 fields of view from at least three independent experiments were randomly chosen. All subjects were informed and agreed. The use of human tissues was approved by the Ethics Committee of Soochow University (SUDA20230619H05).

From public scRNA-seq resources, we selected 11 early LUAD samples and 11 distant normal samples and processed the data according to the reported method [28]. Automated annotation was performed on the obtained cell clusters using the SingleR package (v2.0.0) [29], with all parameters set as default. We confirmed the specific expression of ANGPTL4 and SDC4 in LUAD tissues. Additionally, we isolated epithelial cells from tumor and normal tissues and validated the interaction between ANGPTL4 and SDC4 using CellChat.

RNA-seq and clinical data from patients with LUAD and LUSC samples (stage I-II) were obtained from the Cancer Genome Atlas (TCGA) to assess the prognostic effects of ANGPTL4-SDC4 ligand-receptor pair on early-stage lung cancer. Moreover, we used the cor.test function to identify genes co-expressed with ANGPTL4 and SDC4, and determined their overlapping genes.

For survival analysis, we extracted three expression information from RNA-seq data: individual expression values of ANGPTL4 and SDC4, as well as the average expression value of both. Subsequently, we categorized the samples into high and low expression groups based on the median expression levels of these genes. The survival curves were generated using the Kaplan-Meier formula in the survival package (v3.4.0), and visualized using the ggsurvplot function of the survminer package (v0.4.9).

Store Matrigel matrix, 200 μL pipette tips, and a 96-well plate at 4°C to pre-chill. Once the Matrigel is fully thawed, transfer it to an ice box on a sterile clean bench for handling. Using pre-cooled pipette tips, carefully pipette an appropriate amount of Matrigel and slowly distribute it along the sidewalls of the wells in the 96-well plate, avoiding the formation of air bubbles. Let the plate sit undisturbed for 10 minutes to ensure uniform distribution of the Matrigel. Subsequently, transfer the plate to a 37°C incubator and allow the Matrigel to solidify for at least 30 minutes, providing a stable matrix for subsequent experiments. Digest HUVEC cells to prepare a single-cell suspension, and perform cell counting. Resuspend approximately 2×10 cells in 100 μL of conditioned medium from both the control, R25-30 and R25-50 cells, respectively. Add the cell suspension to the Matrigel-coated wells of the 96-well plate. Allow the cells to settle completely for 3 minutes before proceeding. Observe the wells under a microscope, and capture images of the vascular structures formed by HUVECs every 6 hours. Analyze the vascular structures using the ImageJ Angiogenesis Analyzer plugin to compare the effects [30].

All experiments were independently repeated at least three times, and the data were presented as the mean ± standard deviation. Student's t tests were used for statistical analysis. Probability (p) values less than 0.05 were considered as statistically significant.