Abstract
Space particle radiation is a major environmental factor in
spaceflight, and it is known to cause body damage and even trigger
cancer, but with unknown molecular etiologies. To examine these causes,
we developed a systems biology approach by focusing on the
co-expression network analysis of transcriptomics profiles obtained
from single high-dose (SE) and multiple low-dose (ME) α-particle
radiation exposures of BEAS-2B human bronchial epithelial cells. First,
the differential network and pathway analysis based on the global
network and the core modules showed that genes in the ME group had
higher enrichment for the extracellular matrix (ECM)-receptor
interaction pathway. Then, collagen gene COL1A1 was screened as an
important gene in the ME group assessed by network parameters and an
expression study of lung adenocarcinoma samples. COL1A1 was found to
promote the emergence of the neoplastic characteristics of BEAS-2B
cells by both in vitro experimental analyses and in vivo
immunohistochemical staining. These findings suggested that the degree
of malignant transformation of cells in the ME group was greater than
that of the SE, which may be caused by the dysregulation of the
ECM-receptor pathway.
Keywords: MT: Bioinformatics, radiation biology, transcriptomics,
bioinformatics, network biology, COL1A1
Graphical abstract
graphic file with name fx1.jpg
[41]Open in a new tab
__________________________________________________________________
Hu and colleagues have developed a differential network analysis
framework to examine molecular mechanism in human lung epithelial cells
exposed to two different α-particle radiation methods and found that
ECM-receptor pathway and COL1A1 are key factors in response to multiple
low-dose irradiations, providing insights into space radiation exposure
and protection.
Introduction
Spaceflight imposes a unique suite of environmental effects on
biology.[42]^1 With the expansion of spaceflight, the space environment
includes unique hazards such as radiation exposure, which can adversely
affect biological systems.[43]^2^,[44]^3 The main effect of radiation
exposure is damage to DNA, including base damage, single-strand breaks,
double-strand breaks, chromosomal aberrations, micronuclei, and genomic
instability.[45]^4^,[46]^5^,[47]^6^,[48]^7 At the phenotypic level,
radiation exposure leads to physical damage, including muscle and bone
loss, neurological damage, and impaired immunity.[49]^8 Space radiation
exposure can also result in defects in mitochondrial physiology,
leading to numerous pathological conditions, including bone loss[50]^9
and accelerated aging.[51]^10 The primary concern for radiation
exposure is the risk of developing tumors and/or cancer,[52]^11 and the
known increased cancer risk is caused by chronic low doses of radiation
exposure.[53]^12^,[54]^13 However, the mechanism underlying the
carcinogenic potential of radiation exposure is not fully understood,
especially after chronic low doses of radiation exposure, which is the
situation that astronauts must confront in long-term deep space
missions.[55]^14
Knowledge of spaceflight injury has grown immensely in recent years,
benefiting from large-scale “omics” technologies,[56]^15 making an
understanding of molecular mechanisms at the systems level possible.
In 2020, Gertz et al.[57]^16 integrated transcriptomic and proteomic
data to study how spaceflight affects innate immunity through
mitochondrial processes. This study, together with data from NASA’s
GeneLab ([58]https://genelab.nasa.gov/),[59]^17 opened the door to the
study of spaceflight in the big data
era.[60]^18^,[61]^19^,[62]^20^,[63]^21 Multiomics cancer models
integrating genomics, genome conformation, transcriptomics, and
epigenomics have been used to predict how mutations from spaceflight
radiation and cytoskeletal remodeling/DNA modification aberrations due
to microgravity affect astronauts.[64]^22 Due to the complexity of
spaceflight biology, including different environmental factors such as
particle radiation exposure and microgravity, systems biology
approaches leveraging network modeling are required to characterize
these components in a holistic fashion.[65]^8
Gene function is not isolated but forms protein-protein interaction
networks (PPINs) governed by universal laws, and the network effect of
genes is the driving force moving from the normal to the exposure group
in response to environmental changes.[66]^23 Thus, elucidating the
changes in individual nodes and edges, distinct modules, and the entire
network as a whole in terms of several network centrality parameters
may provide functional information to describe such biological
processes.[67]^24 Differential network analysis can capture changes
between conditions in the interplay between molecules, rather than
changes in single molecules, and has many applications, including
illuminating the molecular basis of tissue-selective processes or
disease-specific genes,[68]^25 and identifying critical signaling
pathways altered during tumor transformation and
progression,[69]^26^,[70]^27 or between different cancer
subtypes.[71]^28 On the other hand, biological networks, including
PPINs, are intrinsically modular, with proteins belonging to the same
module usually sharing common functions.[72]^29^,[73]^30 The
identification of network clusters not only helps us to reveal the
modular organization of cell signaling networks,[74]^31 but also
identifies functional modules in different conditions and in radiation
biology.
In this work, employing transcriptome data associated with two
different α-particle radiation methods (single high-dose irradiation
[SE] and multiple low-dose irradiations [ME]) on human lung epithelial
cells, we applied an unbiased multiple-level differential network
analysis framework for systematically investigating the network changes
that mimic the spaceflight environment ([75]Figure 1). First, two
weighted PPINs were constructed for two different α-particle radiation
groups. Second, the difference between these networks was investigated
through module differential analysis to identify the functional
modules. Third, the functional enrichment and topological analyses were
performed for the functional module, aiming to reveal key pathways and
genes. Finally, the biological functionality of the identified hub
genes was validated through experimental verification.
Figure 1.
[76]Figure 1
[77]Open in a new tab
The differential network analysis framework
Results
Data comparison of exposure groups with the control group and The Cancer
Genome Atlas lung samples
The expression profiles of nine SE samples, nine ME samples, and three
controls were obtained from RNA-sequencing (RNA-seq) ([78]Figure 2A).
Initially, we performed principal-component analyses (PCAs) on the
expression profiles, and the results showed differences between the
control groups and exposure groups ([79]Figure 2B). All the control
groups clustered close together, while the SE and ME samples clustered
into two groups, indicating a relatively high degree of similarity of
the exposure groups, and this expression pattern changed after
exposure, regardless of the exposure methods.
Figure 2.
[80]Figure 2
[81]Open in a new tab
Experimental design and sample information
(A) RNA-seq of BEAS-2B cells after two types of irradiation. Three
replicate controls were BEAS-2B cells without any treatment. SE:
BEAS-2B cells after a single dose of 0.5 Gy and then sub-cultured to
the 10^th passage, 30^th passage, or 50^th passage. ME: BEAS-2B cells
after 25 exposures of 0.02 Gy and then sub-cultured to the 10^th
passage, 30^th passage, and 50^th passage. (B) PCA plot depicting gene
expression profiles for the control, SE, and ME groups. (C) Expression
correlation of SE groups with normal lung samples, and ME groups with
normal lung samples. (D) Expression correlation of R1_05_50 samples
with lung cancer stage I samples and R25_05_50 samples with lung cancer
stage I samples. ^∗p < 0.05, ^∗∗p < 0.01, ^∗∗∗p < 0.001.
The risk model of lung cancer from α-particle radiation has been
primarily established based on BEAS-2B cells.[82]^32 To evaluate the
different exposure methods on BEAS-2B cells, we also performed
correlation analysis between SE groups with normal lung samples and ME
groups with normal lung samples. As shown in [83]Figure 2C, both groups
were significantly correlated with normal lung samples, but the
correlation coefficients for the former comparison were significantly
larger than those for the latter comparison, indicating that the SE
groups were more similar to the control than the ME group; in other
words, the cells under multiple exposures departed more from normality.
Moreover, correlation analysis between R1_05_50 and lung cancer stage I
samples, as well as between R25_05_50 and lung cancer stage I samples
was performed ([84]Figure 2D). The 50^th passage samples were
significantly correlated with the stage I samples, but the correlation
coefficients for the former comparison were significantly smaller than
that for the latter comparison, indicating that the ME groups were more
similar to the stage I samples than the SE group.
Multiple exposures reveal genetic information and tumorigenesis-related
pathways
Network theory provides new insights to study the disturbances of
systems from global and local interactions. Here, we employed WGCNA and
STRING to construct exposure way-specific PPINs (SENet and MENet) to
compare the effects of single and multiple exposures on the cells
([85]Figure 3A). As a result, SENet and MENet share the same network
structure and topology (15,541 nodes and 4,385,019 links) but with
different link weights. Then, we measured the consistency of the two
networks by link weight distributions and RMSIP. As shown in
[86]Figure 3B, the weight distribution patterns of links in SENet and
MENet were almost the same. However, the percentage of links with
weights larger than 0.8 in SENet was larger than in MENet. In the RMSIP
plot ([87]Figure 3C), the RMSIP achieved 0.55, which means that SENet
and MENet have a relatively high topological overlap rate. All the
results indicated that SE groups were very similar to ME groups at the
whole network level, but they still have their own respective
characteristics.
Figure 3.
[88]Figure 3
[89]Open in a new tab
SENet and MENet
(A) Construction of exposure way-specific SENet and MENet PPINs. (B)
Edge weight distribution in SENet and MENet. (C) Consistency between
SENet and MENet. (D) Degree and (E) eigenvector centrality were
significantly larger in MENet than in SENet. (F) Venn diagram showing
overlap of significantly enriched pathways for the top 5% high degree
in SENet (blue) and MENet (pink). (G) The significantly enriched
pathways specific to the SE group (blue) and ME group (red).
^∗p < 0.05, ^∗∗p < 0.01, ^∗∗∗p < 0.001.
We next calculated the topological centralities for genes in SENet and
MENet to compare exposure effects from network features. As shown in
[90]Figures 3D and 3E, although the percentage of links with high
weights in MENet was smaller than that in SENet, genes in MENet had
significantly higher K and EC than those in SENet, which means that
genes in MENet had more interactions and their neighbors also had more
interactions. The results indicated that larger network wiring occurred
with more perturbed genes and interactions induced by multiple
exposures. Then, we investigated the function of hub genes with the top
5% degree in two networks using enrichment analysis of Kyoto
Encyclopedia of Genes and Genomes (KEGG) pathways. There were 171 and
167 pathways that were significantly enriched by the hub genes in SENet
and MENet, respectively, and 153 of them were common pathways
([91]Figure 3F), that is, enriched by the genes from both networks,
which indicated that most hub genes in both networks have similar
functions. However, they also have different pathways ([92]Figure 3G).
For hub genes in SENet, 18 pathways were significantly enriched, and
two of them, the Fanconi anemia pathway and Homologous recombination,
were related to DNA damage repair. When compared with SENet, MENet hub
genes were specifically enriched in 14 pathways, which involved more
genetic information processing-related pathways including Proteasome,
Protein processing in the endoplasmic reticulum, DNA replication, as
well as pathways associated with tumorigenesis, such as ECM-receptor
interaction[93]^33 and Ferroptosis.[94]^34 All the results indicated
that multiple exposures can cause more information exchange and that
these interactions were related to tumorigenesis.
Comparison of single and multiple exposure effects on BEAS-2B cells from the
network core
Important information always flows in the key interactions, which are
always involved in the network core.[95]^35 Therefore, we constructed
the core network of SENet and MENet and compared their differences to
characterize the single and multiple exposure effects on cells. First,
we selected the interactions with weights larger than 0.40 and 0.45 in
SENet and MENet as their core network, named SENetCore and MENetCore,
respectively, to reflect the main important information flows among
genes that were affected by different exposure conditions. SENetCore
consists of 2,123 nodes and 9,847 interactions, while MENetCore
consists of 2,129 nodes and 9,395 interactions. Then, the topological
parameters were calculated for both network cores ([96]Figure 4A;
[97]Figure S1). Similar to the whole network, the eigenvector
centrality (EC) values of nodes in MENetCore were significantly larger
than those in SENetCore, and they distinguished the control groups from
both SE and ME groups ([98]Figure 4B). The results indicated that in
MENetCore, not only hub genes are important, but their neighbors also
play important roles in the network, which indicates that multiple
exposures may induce interactions among genes with a high information
flow context.
Figure 4.
[99]Figure 4
[100]Open in a new tab
Network cores for SENet and MENet
(A) Boxplot representation of eigenvector centrality for SENetCore and
MENetCore. (B) Heatmap representing the expression of genes in
SENetCore and MENetCore. The color bar on the left represents the high
and low (cutoff: median) EC values of genes that are matched, the right
colored bar represents the expression level, and the upper colored bar
represents the group. (C) Significantly enriched pathways for genes in
SENetCore and MENetCore. ^∗p < 0.05, ^∗∗p < 0.01, ^∗∗∗p < 0.001.
To explore the functional differences between SENetCore and MENetCore,
pathway enrichment analysis was performed, and the results demonstrated
that both were enriched in the neuroactive ligand-receptor interaction,
cytokine-cytokine receptor interaction, cAMP signaling pathway, and
others ([101]Figure 4C), indicating that both network cores also have
similar biological functions. However, more importantly, some specific
biological pathways were highlighted. In particular, the ECM-receptor
interaction was significantly enriched in the ME groups at the network
core level, again playing an important role in the ME group. Consistent
with the whole network findings, all the network core results indicate
that MENetCore genes were more closely associated with tumors than
SENetCore genes.
Network analysis reveals the conservation of the ECM-receptor pathway and hub
genes in the ME group
We next performed module analysis for SENetCore and MENetCore to
organize the exposure type-related changes in an unbiased manner. The
Girvan-Newman and label propagation analysis algorithms were used
together to find the robust modules in both networks. With a similar
number of genes, SENetCore and MENetCore were divided into 98 and 76
modules, respectively. As shown in [102]Figure 5A, there were more
modules whose sizes were less than 100 in SENetCore than in MENetCore,
and both networks had three modules with more than 100 genes. The
results indicated that modules in MENetCore are more robust, which may
suggest the biological importance of its functional modules.
Figure 5.
[103]Figure 5
[104]Open in a new tab
Network modules in SENetCore and MENetCore
(A) Module size distribution in the SE and ME groups. (B) The top 10
significantly enriched pathways of genes in the three largest modules
in SENetCore (blue) and MENetCore (red). (C) Venn diagram showing the
overlap of genes in the largest module with dysregulated genes in the
SE and ME groups. (D) Boxplot representation of the degree in the
largest module in SENetCore (blue) and MENetCore (red). (E) Module
genes involved in ECM-receptor pathways. (F) Degree of ECM-receptor
pathway-related genes in the module. (G) Four pathways in ME-Module 1
that COL1A1 are involved. (H) Expression of ECM-receptor
pathway-related genes. ^∗p < 0.05, ^∗∗p < 0.01, ^∗∗∗p < 0.001.
We also applied pathway enrichment analysis to analyze the functional
enrichment of the three largest modules in SENetCore and MENetCore
([105]Figure 5B). The results revealed that the largest module in both
network cores significantly enriched several of the same pathways, such
as axon guidance, basal cell carcinoma and Wnt signaling pathways.
Nevertheless, as previously mentioned, the ECM-receptor interaction was
again enriched by the largest module from the ME group, which indicated
that multiple exposures might indeed affect the ECM-receptor
interaction pathways.
To further explore the effect of exposure types on the modules, we
mapped the differentially expressed genes (DEGs) between the SE group
and control group, as well as the ME group and the control group, on
the corresponding large modules ([106]Figure 5C). In the SE group, 682
DEGs were found in the largest module, but no DEGs and only one DEG was
found in the second and third largest modules. Similarly, in the ME
group, 883 DEGs were mapped in the largest module and no DEGs were
mapped in the second and third largest modules. The results indicated
that the largest module was the key module in the exposure-affected
network, which captured most of the genes associated with expression
changes. Moreover, we also calculated the intra-module topological
parameters for the largest modules, and the results showed that the
module in MENetCore had a significantly higher average degree than the
module in SENetCore ([107]Figure 5D). Other parameters for the largest
three modules are listed in [108]Table S1.
Owing to the conservation of ECM-receptor pathways at different network
levels, we then focused on 14 genes in the largest module of MENetCore
that are involved in this pathway ([109]Figure 5E). Among them, five
genes, THBS2, COL1A1, COL4A1, COL1A2, and LAMC2, were hub genes in the
module ([110]Figure 5F; [111]Table S2). Although THBS2 has the largest
centrality measures in the network level, COL1A1 is involved in more
ECM-related pathways ([112]Figure 5G), serving as a pathway crosstalk.
In addition, most of the ECM-receptor pathway-mapped genes were
up-regulated in the exposure groups compared with the controls
([113]Figure 5H), among which COL1A1 was significantly up-regulated. Of
note, the high degree in network and pathway levels ([114]Figure S2),
as well as the high expression of COL1A1 may suggest its important
biological role, which needs further experimental investigation.
COL1A1 promotes the emergence of neoplastic characteristics in BEAS-2B cells
To further confirm which genes were driving the malignant
transformation process, we examined each key gene in the ECM-receptor
pathway and found that the mRNA level of COL1A1 was most significantly
increased after single and multiple exposures ([115]Figure S3). Then,
the exposure effects of COL1A1 were further validated by bioinformatics
through big data analysis, as well as through experimental expression
data at both the RNA and protein levels. First, the expression level of
COL1A1 was significantly higher in lung adenocarcinoma and lung
squamous cell carcinoma cancer subtypes than in normal tissues, and its
high expression was associated with poor survival outcomes
([116]Figures S4A–S4C). Moreover, its expression level was also
significantly increased in TNM stage samples, T[1-4] samples, N[0-3]
samples, and M[0-1] samples, when compared with normal samples
([117]Figures S4D–S4G). Then, we measured COL1A1 expression levels in
BEAS-2B cells and in the lung cancer cell lines Calu-1, and also in
BEAS-2B cells after single and multiple exposures using qRT-PCR. As
shown in [118]Figures S5 and [119]6A, the expression of COL1A1 was
significantly lower in BEAS-2B cells than in Calu-1 cells and was
significantly increased in BEAS-2B cells after single and multiple
exposures. The results indicated that both single and multiple
exposures can induce tumor-like high expression of COL1A1. In
particular, COL1A1 was also up-regulated in the multiple exposure group
R25_05_50 compared with the single exposure group R1_05_50, which
indicated that the effects of multiple exposure were greater than those
of single exposure on COL1A1 expression.
Figure 6.
[120]Figure 6
[121]Open in a new tab
COL1A1 promotes the emergence of the neoplastic characteristics of
BEAS-2B cells
(A) qRT-PCR of COL1A1 transcripts in SE and ME groups. (B) Western
blots of COL1A1 in the SE and ME groups. (C) Quantitative analysis
results of western blot. (D) BEAS-2B cells were fixed, stained for
COL1A1 (green) and DAPI (blue) and imaged with phase and confocal
microscopy. (E) Zoom in of stained tissue microarrays. (F) The ratio of
possible cells in immunohistochemistry of tissue microarrays. (G)
Transwell invasion assays of BEAS-2B cells, R1_05_P50 cells and
R25_05_P50 cells transfected with NC or siCOL1A1. (H) Quantitative
analysis of transwell invasion assay. I. soft agar colony formation
assay in BEAS-2B cells and COL1A1 knockdown cell lines. (J)
Quantitative analysis of soft agar colony number. ^∗p < 0.05,
^∗∗∗∗p < 0.0001.
At the protein level, we evaluated the COL1A1 expression in SE and ME
groups by western blot and found that the COL1A1 level was also
significantly higher in the exposure compared with the control group
([122]Figures 6B and 6C). In addition, COL1A1 protein was mainly
localized to the cytoplasm of BEAS-2B cell lines ([123]Figure 6D).
Next, we conducted immunohistochemical analysis using a tissue
microarray comprising 32 pairs of lung cancer and adjacent normal
tissues and found that the ratio of COL1A1-positive cells in lung
cancer tissue was significantly larger than in adjacent normal tissue
([124]Figures 6E and 6F; [125]Figure S6).
Finally, to examine the biological role of COL1A1 in the process of
irradiation-induced cell transformation, we knocked down the expression
of COL1A1 in normal, single, and multiple irradiated BEAS-2B cells by
small interfering RNA (siRNA), respectively. The cells with COL1A1
knockdown were established via transfection and screening
([126]Figure S6). We employed a loss-of-function approach using siRNA
specifically targeting COL1A1. Following transfection of
COL1A1-specific siRNA into the cells, the efficacy of gene silencing
was assessed at both mRNA and protein levels ([127]Figure S7). From the
results, cells irradiated with single and multiple exposures of
α-particles showed higher migration and invasion capacity, and in the
COL1A1 knockdown cells, each group with COL1A1 knockdown exhibited
significant decreases in migration and invasion abilities compared with
the corresponding control group, especially in the irradiated group
([128]Figures 6G and 6H; [129]Figure S8). To elucidate the functional
significance of COL1A1 in cellular transformation, we established a
stable knockdown cell line using lentiviral-mediated short hairpin RNA
(shRNA) specifically targeting COL1A1. The effectiveness of this
knockdown strategy was rigorously assessed at both transcriptional and
translational levels ([130]Figure S9), and we also find that in the
COL1A1 knockdown cell line, the number of soft agar colony formation
was significantly reduced ([131]Figures 6I and 6J).
The above results confirmed that irradiation could increase the
expression level of COL1A1 at both RNA and protein levels regardless of
the irradiation type, and the COL1A1 level was higher in the multiple
exposure group. Therefore, we proposed that COL1A1 may be the key
factor for irradiation-induced tumorigenesis, especially under multiple
exposure conditions.
Discussion
Utilizing the transcriptomics induced by α-particle radiation and
protein-protein interaction data, we have proposed a network-based
method to try to answer the fundamental question in radiation biology,
namely, which effects are induced by high- and low-dose radiation
exposure and how this affects the assessment of low-dose cancer risk.
First, two exposure-specific networks were constructed by combining
gene expression and protein-protein interaction data. Then, a
comparison analysis of exposure-specific networks at the global, core
network, and module levels was performed to elucidate the molecular
mechanisms underlying different irradiation types at the systems level.
In addition, we found surprising links between network topology and
α-particle exposure phenotypes. Although the two global PPINs of SENet
and MENet share similar organizational landscapes, two network
parameters, degree and EC, could capture special biological indicators.
Several studies suggest that nodes with a high degree are likely to be
encoded by essential genes, a phenomenon termed the
centrality-lethality rule.[132]^36 In contrast to the degree, EC takes
into account the influence of a node’s neighbors, which is also a good
indicator to identify hub biological pathways and genes.[133]^37
According to its mathematical definition, the higher EC in MENetCore
suggests that the network has more robust modules, which are involved
in more signal transduction processes. The network results also suggest
that multiple low doses of radiation exert a stronger pathogenic
effect, especially via the ECM-receptor interaction pathway to affect
cancer risk.
Research on astronaut health and model organisms has revealed
six fundamental biological features resulting from radiation
exposure, including oxidative stress, DNA damage, mitochondrial
dysregulation, epigenetic changes, telomere length alterations, and
microbiome shifts.[134]^11 Here, we found that the ECM-receptor pathway
was conserved from the whole network to the core network and was the
largest module in the ME group. The biophysical characteristics of the
ECM strongly regulate cellular responses and are used as indicators of
cancer progression and metastasis[135]^38; hence, it is vital to
research the effect of the ECM on tumor development.[136]^39 We thus
suggest that the ECM, as a novel cancer hallmark, may be an additional
biological feature for multiple low-dose radiation exposures. However,
other studies have found that ECM-related pathways are affected by
spaceflight[137]^16 and particle radiation.[138]^40 Illa-Bochaca
et al.[139]^41 found that particle radiation induces more aggressive
tumors and that its carcinogenic effect is mediated by stromal
remodeling as well as some signaling pathways in the microenvironment,
such as Notch signaling. Yet, how ECM-related pathways contribute to
the carcinogenic effect of particle radiation remains unclear. The
present work is the first to systematically investigate the
transcriptomic effects induced by particle radiation delivered by
differential exposure methods, and the results suggest that multiple
radiation exposures may affect the expression of collagens, resulting
in ECM remodeling in human bronchial epithelial cells.
Collagen is the primary constituent of lung tissue proteins and may be
responsible for a plethora of severe lung diseases.[140]^42 Among
collagens, type I is the most prominent within the ECM and has a
crucial role in sustaining the lung’s overall structure[141]^43;
furthermore, COL1A1 serves as a potential prognostic marker and
therapeutic target in lung cancer.[142]^44 In our work, the biological
role of COL1A1 in multiple low-dose radiation exposures was further
validated by the bioinformatics analyses of lung cancer public
databases, a series of biological experiments, and clinical sample
research.
Based on our data and observations, we propose possible molecular
etiologies of α-particle radiation. As shown in [143]Figure 7, in cells
after multiple exposures of low-dose radiation, the expression of
COL1A1 was higher, was associated with wider fibers, and caused
collagen protein fibers to gradually thicken and change to a linear
shape compared with cells with a single exposure of high-dose
radiation. The linearized collagen fibers are harder than their coiled
counterparts, resulting in an increase in ECM stiffness. This
stiffening can significantly enhance ECM density and disorder after
irradiation, leading to malignant cell transformation. Due to the
pathophysiological and therapeutic potential of both the ECM-related
pathway[144]^33 and collagens[145]^45 in cancer, we suggest that COL1A1
may serve as a cancer biomarker in multiple low-dose radiation
exposure.
Figure 7.
[146]Figure 7
[147]Open in a new tab
Model of the impact of α-particle radiation on human bronchial
epithelial cells
High expression of COL1A1 remodels the ECM organization, increasing ECM
density and disorder after multiple low doses of radiation.
Overall, our work highlights the power of using the differential
network method with “omics” data to understand the fundamental
biological problem in radiation biology. We found that ECM-receptor
pathway dysregulation is a central hub for the effect of low-dose
radiation. This new biological feature not only highlights the
ECM-receptor pathway together with COL1A1 as a cancer risk in
spaceflight, but also suggests a possible approach for pharmaceutical
interventions and studies in space biology.
Materials and methods
Cell line
The human bronchial epithelial cell line BEAS-2B was purchased from the
American Type Culture Collection (CRL-9609; Rockville, MD, USA) and
maintained in DMEM (Gibco, Grand Island, NY, USA) containing 10% fetal
bovine serum (FBS) (Gibco Invitrogen, Carlsbad, CA, USA), penicillin
(100 U/mL) and streptomycin (100 μg/mL) in a fully humidified incubator
with 5% CO[2] at 37°C.
α-particle irradiation
The α-irradiator that we used in this study consisted of a ^241Am
source, a rotating radiation source holder, a sample holder, and other
necessary parts.[148]^46 The ^241Am source emitted α-particles at a
dose rate of 0.14 Gy/min. The BEAS-2B cells were irradiated either with
a single dose of 0.5 Gy one time or were irradiated with 0.02 Gy (9 s
of irradiation time) once every 3 days for 25 exposures. Meanwhile, a
blank control was set up, which was sub-cultured together with the
irradiated group. Thus, the three groups of cells were labeled as the
Ctrl (blank control) group, the single exposure (SE) group (0.5 Gy),
and the multiple exposures (ME) group (25 single-dose exposures of 0.02
Gy). After all exposures were completed, the cells were continuously
sub-cultured to the 10^th passage (approximately 3 weeks), 30^th
passage (approximately 9 weeks), or 50^th passage (approximately
15 weeks) and harvested for RNA-seq at various passages. Thus, the SE
groups included R1_05_10, R1_05_30, and R1_05_50, and the ME group
included R25_05_10, R25_05_30, and R25_05_50.
RNA-seq data profiling and analysis
All cells were harvested in TRIzol reagent (Invitrogen). Total RNAs
were reverse-transcribed using the PrimeScript RT Reagent Kit (Takara,
Kusatsu, Shiga, Japan). Then, transcriptomics profiling was performed
by RNA-seq analysis using an Illumina Novaseq 6000 system by OE Biotech
(Shanghai, China). The count level RNA-seq data was used and genes with
low expression (count = 0 across more than 80% of samples) were
filtered out. rlog normalization and DEG analysis were performed on
filtered data using the R package “DESeq2.” Genes with an adjusted p
value <0.05 and |fold change| >2 were selected as significant DEGs. PCA
was performed on the expressional matrix, and all samples were
projected onto the two dimensions determined by the first two PCs.
Exposure way-specific PPIN construction
To compare the effects of different exposure methods on the cell
ecosystem, we constructed two exposure way-specific PPINs for the SE
group and ME group, by integrating the gene co-expression network
constructed by WGCNA[149]^47 and the PPIN downloaded from
STRING.[150]^48 First, the topological overlap matrix (TOM) was
calculated, and a weighted gene co-expression network was constructed
using the R package WGCNA based on the gene expression profile of the
control group and one of the exposure groups with a soft threshold of
16 for the SE group and eight for the ME group. Next, the network was
combined with the human PPIN downloaded from STRING to gain the final
exposure way-specific PPIN for the SE group and the ME group, named
SENet and MENet, respectively.
Network consistency computation
To measure the consistency of SENet and MENet networks, we calculated
the root-mean-square inner product (RMSIP), which measures similarity
in the subspace spanned by given eigenvectors, for the TOMs of the two
networks.[151]^49 As shown in [152]Equation 1, RMSIP was calculated
from the eigenvectors of TOMs.
[MATH: RMSIP=(1n(∑i=1n∑j=1n(Vi.
Uj)2)1/2 :MATH]
(Equation 1)
where n is the number of vectors of the matrix, and V[i] and U[j] are
the i^th and j^th eigenvectors of the TOM[SE] and TOM[ME],
respectively. We traversed the first 600 eigenvectors of the TOM to
calculate the corresponding RMSIP.
Exposure way-specific PPIN core construction and analysis
To reveal the efficient interactions among genes in SENet and MENet, we
filtered out the edges with low weight values to keep the high-weighted
links and construct the core networks of SENet and MENet, namely,
SENetCore and MENetCore. Here, we kept the links with weights larger
than 0.4 and 0.45 to construct SENetCore and MENetCore, respectively.
The R package “igraph” was employed to calculate the topological
parameters degree (K), betweenness (B), closeness (C), and eigenvector
centrality (EC) for each node.[153]^50 Their definitions are listed in
the supplemental methods.
Robust network module identification
Since the implementation of important biological functions is always in
the form of a “module,” we identified the robust module in the exposure
way-specific PPIN core to indicate the biological function of the
network by combining the top-down and bottom-up module discovery
methods.[154]^51 Specifically, the weighted Girvan-Newman algorithm and
label propagation analysis were used to generate the modules for the
network core, and then a hypergeometric test was performed for each
pair of modules to find the similarity in the two model results. The
overlapping modules with p values <0.05 were considered as robust
modules.
Functional enrichment and survival analysis
Functional enrichment analysis and visualization were performed using
the R package “clusterProfiler”.[155]^52^,[156]^53 Gene Ontology terms
with adjusted p values <0.01 were considered significantly enriched,
while the KEGG pathways with p values <0.05 were considered
significantly enriched pathways. The survival analysis was conducted by
the R packages “survival” and “survminer” based on The Cancer Genome
Atlas (TCGA) data.
siRNA mediated COL1A1 gene knockdown
BEAS-2B cells were seeded at 60%–70% confluence prior to transfection.
COL1A1-specific siRNA and scrambled siRNA (negative control) were
purchased from Sangon Biotech (Shanghai). Transfections were carried
out using Lipofectamine RNAiMAX (Thermo Fisher Scientific) according to
the manufacturer’s protocol. Cells were harvested 48 h
post-transfection for downstream analyses. The target sequences for
COL1A1 were as follows: siRNA-1 GGCAAGACAGTGATTGAATAC, siRNA-2
CAAAGGAGACACTGGTGCTAA, siRNA-3 AACCGGAGAACTTACAATAC.
Construction of COL1A1 knockdown cell line using shRNA
A lentiviral vector expressing shRNA targeting COL1A1 was designed and
synthesized by Sangon Biotech (Shanghai). The shRNA sequence was
selected based on its efficiency in silencing COL1A1 without off-target
effects, the TargetSeq of sh-COL1A1-1 is ACAGGGCGACAGAGGCATAAA, the
TargetSeq of sh-COL1A1-2 is CGATGGATTCCAGTTCGAGTA. The construct
contained a GFP marker for selection. Virus production and infection
HEK293T cells were transfected with the lentiviral packaging plasmids
(pMD2.G and psPAX2) along with the COL1A1-shRNA or scrambled shRNA
vector using calcium phosphate transfection. After 48 h, viral
supernatants were collected, filtered, and used to infect BEAS-2B cells
in the presence of polybrene. Stable transductants were screened by
puromycin resistance and verified for GFP expression.
qRT-PCR analysis
A 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. The PCR amplification was carried
out using a Life Technologies system (Vii7A) and initiated by 2 min at
95°C before 40 thermal cycles, each consisting of 30 s at 95°C and 40 s
at 62°C, with a final extension of 10 min at 72°C. Data were analyzed
by Ct value comparison method and normalized by control expression in
each sample. The primer sets are listed in [157]Table S3.
Western blotting
Cells were harvested and lysed using RIPA buffer (Beyotime, Shanghai,
China). Samples were sonicated and centrifuged at 12,000 × g for 15 min
at 4°C. The concentration of total protein was determined using the DC
Protein Assay Kit I (Bio-Rad, Richmond, CA, USA). The samples were then
denatured at 100°C for 5 min. Total proteins were separated by 12%
SDS-PAGE and transferred to Hybond nitrocellulose membranes (Amersham,
Piscataway, NJ, USA). The membranes were blocked with 5% nonfat milk
powder in Tris-buffered saline (pH 7.5) and hybridized overnight with
primary antibodies against COL1A1 and β-actin (Cell Signaling
Technology, Danvers, MA, USA). Finally, the membranes were incubated
with horseradish peroxidase-conjugated anti-immunoglobulin (Ig)G for
1 h at room temperature and visualized with an ECL kit (Millipore,
Billerica, MA, USA). Protein expression levels were normalized to the
loading controls based on their intensity analyzed with ImageJ
(National Institutes of Health, Bethesda, MD, USA).
Immunofluorescence
BEAS-2B cells were seeded on coverslips in 12-well plates. After
irradiation, cells were fixed with 4% formaldehyde in PBS at room
temperature for 10 min and then in methanol at −20°C for 20 min, then
were permeabilized in 0.1% Triton X-100 in PBS for 10 min and blocked
with 5% skim milk for 1 h. Cells were then incubated at room
temperature for 2 h with anti-COL1A1 before staining with Alexa Fluor
555 anti-rabbit antibody (Molecular Probes, Eugene, OR, USA) at 37°C
for 1 h. Following extensive washing in PBS, the cells were mounted on
slides using DAPI mounting medium (Molecular Probes). The stained cells
were observed under an Olympus IX71 microscope (Olympus, Tokyo, Japan)
and also under a laser scanning confocal microscope (Olympus FV1200,
Tokyo, Japan) at Soochow University (Suzhou, China).
Immunohistochemistry of clinical tissue
Non-small cell lung cancer tissue microarrays and patient pathological
information were provided by the Department of Pathology of the Second
Affiliated Hospital of Soochow University. These lung cancer tissues
and matched adjacent normal tissues were reviewed and approved by the
Department of Pathology of the Second Affiliated Hospital of Soochow
University. All staining was assessed by a quantitative imaging method,
and the percentage of positive immunostaining and the staining
intensity were recorded. The H-score was calculated using the following
formula: H-score = (percentage of cells of weak intensity
[MATH: × :MATH]
1) + (percentage of cells of moderate intensity
[MATH: × :MATH]
2) + (percentage of cells of strong intensity
[MATH: × :MATH]
3).
Cell migration and invasion assays
BEAS-2B cells were plated into 6-well plates and infection was
performed when the confluence reached 30%. Briefly, the culture medium
was replaced with a mixture of 1 mL of fresh medium and 1 mL of
lentivirus suspension. After infection for 24 h, the cells were
screened using medium containing 2 μg/mL puromycin (Invitrogen,
Carlsbad, CA, USA) for 1 week to establish the cell line with stable
COL1A1 knockdown. The lentivirus for COL1A1 knockdown was designed and
packed by RiboBio Biotechnology (Guangzhou, China). The stable COL1A1
knockdown was verified by western blot analysis.
The migratory capacity of cells was evaluated using a wound-healing
assay. Exposed cells were seeded at a density of 5 × 10^5 cells/mL/well
in a six-well plate. When the cells reached confluence, a scratch was
made using a 200-μL pipette tip. Detached cells were washed off and the
remaining cells were treated with fresh media without FBS for 48 h.
Pictures were taken under a light microscope and the gap width was
quantified by using ImageJ software. To evaluate the invasiveness of
control and irradiated cells, commercial transwell chambers (24-well
plates, 8-μm pore size; BD Biosciences, Franklin Lakes, NY, USA) were
used. Briefly, the upper surface of the chambers was coated with 10 μL
of matrigel (Corning, Corning, NY, USA). The cells, suspended in
serum-free medium, were seeded at a density of 5 × 10^4 cells/mL/well
in the upper chamber, and 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 at the bottom
surface were stained with crystal violet. Stained cells were extracted
with 0.5 mL of extraction solution (50:49:1 of ethanol:distilled
water:1 M HCl) per well for 10 min with shaking. The optical density of
each well was measured at 540 nm.
Lung cancer data collection and analysis
A dataset containing gene expression profiles of 492 patients with lung
cancer and 59 normal samples from TCGA was employed for the comparison
study. A total of 286 patients with stage I, 122 patients with stage
II, and 84 patients with stage III were included in this dataset. The
count RNA-seq data were directly used for further analysis.
Statistics
All experiments were independently repeated at least three times, and
all data are presented as means ± standard error. Student’s t tests
were employed for statistical analysis, and a probability (p) value
less than 0.05 was considered statistically significant.
Data and code availability
The α-particle irradiation-induced RNA-seq data have been deposited in
NCBI’s Gene Expression Omnibus ([158]https://www.ncbi.nlm.nih.gov/geo/)
under the accession number [159]GSE235882
; lung adenocarcinoma and lung squamous cell carcinoma TNM sample data
were downloaded from TCGA. Deposited data and all code are available at
GitHub ([160]https://github.com/CSB-SUDA/RDNA).
Acknowledgments