Low-Dose Radiation Affects Cardiovascular Disease Risk in Human Aortic Endothelial Cells by Altering Gene Expression under Normal and Diabetic Conditions

High doses of ionizing radiation can cause cardiovascular diseases (CVDs); however, the effects of <100 mGy radiation on CVD remain underreported. Endothelial cells (ECs) play major roles in cardiovascular health and disease, and their function is reduced by stimuli such as chronic disease, metabolic disorders, and smoking. However, whether exposure to low-dose radiation results in the disruption of similar molecular mechanisms in ECs under diabetic and non-diabetic states remains largely unknown; we aimed to address this gap in knowledge through the molecular and functional characterization of primary human aortic endothelial cells (HAECs) derived from patients with type 2 diabetes (T2D-HAECs) and normal HAECs in response to low-dose radiation. To address these limitations, we performed RNA sequencing on HAECs and T2D-HAECs following exposure to 100 mGy of ionizing radiation and examined the transcriptome changes associated with the low-dose radiation. Compared with that in the non-irradiation group, low-dose irradiation induced 243 differentially expressed genes (DEGs) (133 down-regulated and 110 up-regulated) in HAECs and 378 DEGs (195 down-regulated and 183 up-regulated) in T2D-HAECs. We also discovered a significant association between the DEGs and the interferon (IFN)-I signaling pathway, which is associated with CVD by Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, protein–protein network analysis, and module analysis. Our findings demonstrate the potential impact of low-dose radiation on EC functions that are related to the risk of CVD.


Introduction
With industrial and scientific developments, the risk of radiation exposure has increased in all parts of our lives. Issues such as occupational radiation exposure, the future of nuclear power, manned space flights, and the threat of radiological terrorism call for a thorough understanding of the health risks associated with low-dose radiation exposure [1]. People also experience low-dose radiation exposure from the medical use of radiation for diagnostic purposes. According to the International Commission on Radiological Protection (ICRP) recommendation, those at risk for repeated radiation exposure include health care and nuclear industry workers, who are typically monitored and restricted to effective doses of 100 mSv every five years (i.e., 20 mSv per year), with a maximum of 50 mSv allowed in any given year [2,3]. In contrast, radiation exposure in patients who commonly undergo multiple medical imaging procedures is not typically monitored [4]. While experimental and epidemiologic evidence has linked exposure to low-dose ionizing radiation with the development of solid cancers and leukemia, the association between long-term risk of cardiovascular disease (CVD) and low-dose radiation exposure is unknown [5]. Since the biological mechanisms underlying CVD after low-dose are unclear and results from epidemiological studies are inconsistent, only a weak relationship between CVD and low-dose radiation exposure has been reported. Therefore, molecular studies should be conducted to improve our understanding of the pathogenesis and risk estimation of radiation-induced CVD at low doses.
Endothelial cells (ECs) are a key component of the cardiovascular system, and their function is diminished by the presence of cardiovascular risk factors. Several stimuli, including chronic disease states, metabolic conditions (e.g., type 2 diabetes mellitus (T2DM), obesity, dyslipidemia), smoking, and disturbed blood flow drive EC dysfunction [6]. There is a particularly close link between T2DM and CVD. Progression of T2DM eventually involves the development of chronic vascular impairment, which results in CVD. These cardiovascular complications are the main causes of death in patients with diabetes mellitus (DM) worldwide [7,8]. It has been reported that diabetes emerged as a late effect of radiation therapy among childhood cancer survivors who underwent total body irradiation or radiation to the head and the abdomen [9][10][11][12][13][14][15]. In addition, based on epidemiological studies of A-bomb survivors, an association between radiation and the incidence of diabetes was hypothesized [16]. Interestingly, the therapeutic effect of repeated low-dose radiation (25 or 50 mGy) on diabetes-induced cardiopathy has been reported; low-dose radiation ameliorated DM-induced cardiopathy caused by anti-oxidant exposure [17] and prevented cardiac inflammation and pathological remodeling [18]. The major mechanisms of action and therapeutic potential of low-dose radiation on chemical-or radiation-induced tissue damage were hypothesized to include the induction of hormesis and adaptive responses [19], suggesting that low-dose radiation may promote different responses in healthy or diseased conditions. The risk of disease rises over time as the number of metabolic syndrome characteristics increases; therefore, early intervention is warranted.
In this study, we investigated the effect of low-dose radiation on EC function and pathophysiology using RNA sequencing (RNA-seq) to identify its molecular mechanisms using primary HAECs derived from healthy or T2DM donors. This study provides valuable insights into the effects of combining low-dose radiation with other disease risk factors, including metabolic syndrome and DM, on CVD.
We performed RNA-seq analysis using total RNA isolated from HAECs and T2D-HAECs treated with or without 100 mGy of ionizing radiation. Whole transcriptome analysis showed that T2D affected the normal EC transcriptome and that this T2D-induced transcriptional alteration was greater than that in ECs treated with 100 mGy of ionizing radiation. To confirm the EC gene signatures affected by T2DM, we analyzed the gene profiles of T2D-HAECs compared to those of HAECs. A total of 27,685 genes were analyzed; 4701 differentially expressed genes (DEGs) with p-values < 0.05 and fold changes > 1.2 were found (Table S1). These DEGs were represented in a heatmap with hierarchical clustering ( Figure 2A). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that T2DM in ECs most significantly altered the pathways associated with the ribosome, DNA replication, cell cycle, oocyte meiosis, and purine metabolism ( Figure 2C, Table S2). Gene ontology (GO) analysis of biological processes showed that the DEGs were most significantly involved in G-protein-coupled receptor activity, sensory perception, ribonucleoprotein complex biogenesis, ribosome biogenesis, and the mitotic cell cycle in T2D-HAECs ( Figure 2D, Table S3). . Quantification of each protein is shown in the graph. All data represent th mean ± standard deviation (SD) of each independent experiment; * p < 0.05; ** p < 0.01; *** p < 0.005 ns = no significance.
We performed RNA-seq analysis using total RNA isolated from HAECs and T2D HAECs treated with or without 100 mGy of ionizing radiation. Whole transcriptome anal ysis showed that T2D affected the normal EC transcriptome and that this T2D-induced transcriptional alteration was greater than that in ECs treated with 100 mGy of ionizing radiation. To confirm the EC gene signatures affected by T2DM, we analyzed the gen profiles of T2D-HAECs compared to those of HAECs. A total of 27,685 genes were ana lyzed; 4701 differentially expressed genes (DEGs) with p-values < 0.05 and fold changes . Quantification of each protein is shown in the graph. All data represent the mean ± standard deviation (SD) of each independent experiment; * p < 0.05; ** p < 0.01; *** p < 0.005; ns = no significance.

Gene Profiling of HAECs after 100 mGy Ionizing Radiation Treatment
Next, we analyzed low-dose radiation-responsive genes in normal HAECs; results revealed 243 DEGs with a significance threshold of p-value < 0.05 and fold change >1.2. When compared with that in untreated cells, 110 DEGs were up-regulated and 133 were down-regulated after irradiation (Table S4). The 20 most altered genes by fold change are listed in Table 1. A heatmap analysis with hierarchical clustering of DEGs responsive to low-dose radiation in HAECs is shown in Figure 3A. KEGG pathway analysis showed that DEGs involved in the RIG-I-like receptor signaling pathway and taste transduction were most significantly altered in HAECs exposed to 100 mGy of radiation ( Figure 3B) and involved in GO terms associated with viral infection and innate immune-related biological processes ( Figure 3C). Protein-protein interaction (PPI) network analysis indicated proteins from 100 mGy-irradiated normal HAECs, represented by 186 nodes and 293 edges between 243 DEGs ( Figure 3D). To investigate the interactions between PPI networks, functional enrichment analysis was performed; results revealed that PPI networks from 100 mGy-irradiated normal HAECs involved interferon (IFN) alpha/beta signaling, negative regulators of DDX58/IFIH1 signaling, and ISG15-protein conjugation. Together, these results showed that low-dose radiation could influence transcriptomes and HAEC signaling pathways.   (B) Significant KEGG functional pathways associated with the identified DEGs in HAECs treated with and without low-dose radiation (p < 0.05, FDR q-value < 0.25). (C) Significant GO terms of biological processes associated with the DEGs identified in HAECs treated with and without low-dose radiation (p < 0.05). (D) The protein-protein interaction (PPI) networks in HAECs treated with and without low-dose radiation were visualized using the STRING plug-in of the Cytoscape program (version 3.9.1). Node colors represent the expression. The gradual change in color from blue to red corresponds to the change in the expression level (from down-regulation to up-regulation) in low-dose radiation-treated HAECs versus non-treated HAECs.

Gene Profiling in T2D-HAECs after Treatment with 100 mGy of Ionizing Radiation
We also analyzed low-dose radiation-induced genes in human ECs in a diabetic state. We identified 378 DEGs with a significance threshold of p-value < 0.05 and fold change > 1.2, as shown in the heatmap in Figure 4A. When T2D-HAECs were treated with 100 mGy of ionizing radiation, 183 DEGs were up-regulated and 195 were down-regulated when compared to that in untreated cells (Table S5). The 20 most significantly modulated genes by fold change are listed in Table 2. KEGG pathway analysis of T2D-HAECs showed that low-dose radiation altered a pathway involved in antigen processing and presentation ( Figure 4B), while the GO term for post-translational protein folding was also identified in these cells ( Figure 4C). PPI networks identified possible contacts between proteins encoded by 378 DEGs that contained 284 nodes and 291 edges in 100 mGy-irradiated T2D-HAECs ( Figure 4D).

Regulation of Cell Function Following Response to Low-Dose Radiation in ECs under Normal and Diabetic Conditions
We found 15 common DEGs that were significantly altered in both HAECs and T2D-HAECs after exposure to low-dose radiation (Table 3). No KEGG pathways were identified for these DEGs; however, 15 DEGs were related to two GO terms associated with the The PPI networks associated with low-dose radiation responsiveness in T2D-HAECs visualized using the STRING plug-in of the Cytoscape program (version 3.9.1). Node colors represent the expression. The gradual change in color from blue to red corresponds to the change in the expression level (from down-regulation to up-regulation) in low-dose radiation-treated HAECs versus non-treated HAECs.

Regulation of Cell Function Following Response to Low-Dose Radiation in ECs under Normal and Diabetic Conditions
We found 15 common DEGs that were significantly altered in both HAECs and T2D-HAECs after exposure to low-dose radiation (Table 3). No KEGG pathways were identified for these DEGs; however, 15 DEGs were related to two GO terms associated with the negative regulation of viral genome replication and cellular response to type I IFN ( Figure 5A). The PPI network analysis showed that 12 of these 15 DEGs closely interacted ( Figure 5B). We validated gene sets in both HAECs and T2D-HAECs treated with 100 mGy irradiation by qRT-PCR ( Figure 6). Consistent with the RNA-seq results, ACKR4, IFIH1, and LAP3 were up-regulated in HAECs but down-regulated in T2D-HAECs. However, the mRNA expression of CMPK2, CXCL10, IFI35, IFT1, ISG15, RSAD2, and USP18 was up-regulated in HAECs after low-dose radiation, while that of CMPK2, CXCL10, IFT1, ISG15, RSAD2, and USP18 was up-regulated in T2D-HAEC, which was not consistent with the RNA-seq results. The other genes (ASTE1, IFIT3, HERC6, TASR13, and TAS2R20) were not altered. We examined the cell proliferation rate, cellular senescence, and tube formation activity in HAECs and T2D-HAECs treated with 100 mGy and 2 Gy of radiation and found that low-dose radiation increased the number of β-galactose-positive cells and reduced tube formation activity, similar to that seen in non-irradiated T2D-HAECs (Figure 7). However, both cell types showed the same responses to high-dose radiation exposure, indicated by their lower proliferative rates, increased cellular senescence, and reduced tube formation compared with non-irradiated HAECs or T2D-HAECs, respectively. and found that low-dose radiation increased the number of β-galactose-positive cells and reduced tube formation activity, similar to that seen in non-irradiated T2D-HAECs (Figure 7). However, both cell types showed the same responses to high-dose radiation exposure, indicated by their lower proliferative rates, increased cellular senescence, and reduced tube formation compared with non-irradiated HAECs or T2D-HAECs, respectively.

Discussion
Although public concern about low-dose radiation exposure has increased with technological development and the increased application of radiation technology in the fields of medicine and industry, clinically relevant tissue damage does not occur below an absorbed dose of 100 mGy; this forms the basis of current radioprotection systems against non-cancer effects. However, recent epidemiological findings point to a risk of non-cancerous disease after exposure to lower doses of ionizing radiation than previously thought, especially for CVD and cataracts [25]. Owing to limited statistical support, the relationship between dose and risk is undecided below 0.5 Gy [26]. However, emerging evidence suggests that doses under 0.5 Gy may increase the long-term risk of CVD [27]. Moreover, environmental factors, including radiation, can accelerate the onset and progression of diabetes with more serious side effects [28]. Similar to diabetes, radiation causes damage to blood vessels and tissues and affects aging [29]. Extensive clinical and experimental evidence suggests that endothelial dysfunction occurs in T2DM and prediabetes patients [20]. Endothelial dysfunction has received increased attention as a potential contributor to the pathogenesis of vascular disease in T2DM [30].
In this study, we used gene profiling analysis to investigate the possible impact of low-dose radiation on the induction of endothelial dysfunction. Impairment of endothelial cell function is a critical risk factor for both cardiovascular pathology and DM. Thus, we analyzed low-dose radiation-responsive gene profiles in irradiated ECs under normal and DM conditions. We performed RNA-seq using total RNA from HAECs and T2D-HAECs exposed to 100 mGy of low-dose radiation and profiled genes that may be associated with endothelial damage and dysfunction that cause CVD. We confirmed that T2D-HAECs exhibited impaired EC characteristics and gene profiles that were strikingly different from those of HAECs. We observed that low-dose radiation of normal ECs promoted a phenotype similar to that of non-irradiated diabetes, which showed increased numbers of β-galactose-positive cells and reduced tube formation, further indicating impaired EC function. Additionally, 2 Gy of high-dose radiation induced more significant EC dysfunction, including reduced proliferation rate, cellular senescence, and tube formation in both HAECs and T2D-HAECs when compared to those of the controls. However, in T2D-HAECs, low-dose radiation did not induce significant changes in EC function. Indeed, similar to the EC functional assessment, HAECs showed more consistent results in qRT-PCR validation assays of the RNA-seq results compared with T2D-HAECs. qRT-PCR analysis showed that low-dose irradiation of T2D-HAECs resulted in only three reproducible gene products compared to that with RNA-seq results. Similar to the results of EC functional assessments, low-dose radiation had no significant effect on the qRT-PCR outcome. These different responses between HAECs and T2D-HAECs to low-dose radiation may be caused by inadequate radiation to phenotypically or mechanistically alter HAECs under diabetic conditions. Our findings were similar to those of Vieria Dias et al., who investigated the response of HAECs to different low-dose (6 mGy/h) and high-dose (1 Gy/min) rates of 0.05, 0.5, 1.0, and 2 Gy [31]. They reported that the threshold for angiogenic capacity loss in cells treated with both low-and high-dose rates of radiation was between 0.5 and 1.0 Gy [31]. However, the effect of 100 mGy under diabetic conditions has not been investigated. Zhang et al. tested whether 14 repeated exposures of 12.5, 25, or 50 mGy every two days for two weeks could protect streptozotocin-induced diabetic C57BL/6J mice against diabetesinduced cardiopathy and reported that only repeated 25 and 50 mGy exposures exerted protective effects in the form of not only anti-apoptotic and anti-oxidant effects but also in preventing cardiomyocyte hypertrophy and fibrosis [17]. However, Zhang et al. reported that irradiating mice with doses higher than 100 mGy did not show any significant effects on diabetes-induced cardiopathy, including anti-oxidant enzyme levels or histopathological alterations to the heart [17]. Takahashi et al. reported that among acute exposures to doses of 0.25, 0.5, 1.0, and 2.0 Gy, only exposure to acute doses of 0.5 Gy provided significant protection against ALX-induced diabetes, while 0.25 Gy was not adequate to exert a protective effect against the development of diabetes [32]. While it has been reported that low-dose radiation induced damage in normal ECs in the form of morphological changes, proliferation, and DNA damage [31,33], no studies have examined the effect of low-dose radiation (<100 mGy) on ECs under diabetic conditions; thus, further studies are needed to estimate the risk of CVD by providing evidence of the biological mechanism involved in its pathogenesis.
CVD and T2DM share several common pathophysiological features, including insulin resistance, inflammation, oxidative stress, hypercoagulability, high blood pressure, dyslipidemia, and obesity [34]. We investigated the radiation-induced pathogenic mechanism linked to CVD and diabetes using bioinformatic analysis. Examining the impact of low-dose radiation on EC dysfunction indicated potential CVD risk, while PPI network analysis revealed a hub of 12 genes (ACKR4, CMPK2, CXCL10, HERO6, IFI35, IFIH1, IFT1, IFT3, ISG15, LAP3, RSAD2, and USP18) out of 15 DEGs that were associated with IFN-I signaling. IFNs constitute a family of cytokines that is divided into three main subtypes: type I, II, and III IFNs. The type I IFNs (IFN-α and IFN-β) are reported to be important mediators in atherosclerosis [35]. Using qRT-PCR, we confirmed that low-dose radiation significantly altered mRNA levels of ACKR4, IFIH1, CMPK2, CXCL10, IFI35, IFT1, ISG15, LAP3, RSAD2, and USP18, which are involved in IFN-I signaling. It is known that IFIH1 (interferon-induced helicase-1), IFI35 (interferon-induced protein 35), and USP18 (ubiquitinspecific peptidases 18) activate the IFN signaling pathway; IFIH1 induces pro-inflammatory cytokines, and type I IFNs respond to viral infections in which they act as innate immune receptors [36,37]. It has been suggested that viral infections can cause autoimmune diseases such as type 1 diabetes (T1D) through β-cell disruption [38]; however, this remains debatable. IFIT3 is a novel biomarker for human ischemic cardiomyopathy that can inhibit the nuclear factor-kappa B (NF-κB) pathway in injured arteries, resulting in the inhibition of endothelial cell proliferation, migration, and re-endothelialization [39]. Therefore, IFI35 may induce diverse diseases, including hyperplasia, atherosclerosis, and CVD after endothelial injury [40].
USP18 is expressed in the heart, where it modulates inflammation and apoptosis [39]. USP18 is stimulated by IFNs, which form a negative feedback loop and act as isopeptidases for specific substrates with ISG15, a ubiquitin-like protein [41]. ISG15 is also induced by IFNs and can modify proteins in a ubiquitin-like manner. This modification in cardiomyocytes has been found to contribute to the suppression of viral replication in a CVB-3-infected mouse model; therefore, ISG15 modification may be a critical part of the innate immune system in cardiomyocytes [42]. CMPK2 (cytidine/uridine monophosphate kinase 2) is localized in the mitochondria, where it engages in macrophage differentiation that is involved in atherosclerosis processes. Deletion of CMPK2 has been reported to reduce mtROS production induced by IFN-α, suggesting that it exerts pro-atherogenic effects [43]. CXCL10 (C-X-C motif chemokine ligand 10) and RSAD2 (radical S-adenosyl methionine domain containing 2) are highly expressed in IFN-γ and LPS-treated vascular smooth muscle cells, which can induce pro-inflammatory and pro-atherogenic processes [44]. ACKR4 (atypical chemokine receptor 4) belongs to the ACKR subfamily and acts as a chemokine receptor for CCL2, CCL8, CCL13, CCL19, CCL21, and CCL25 [45]. It has been reported that the combination of chemokines and ACKR4 can lead to the recruitment of β-arrestin [46][47][48], which mediates cardioprotective signaling that results in ligand internalization [49]. LAP3 (leucine aminopeptidase 3) has been reported to play a vital role in the pathogenesis of nonalcoholic fatty liver disease (NAFLD), which is a known risk factor for T2DM [50] and is associated with cardiovascular and cerebrovascular diseases [51].
Our study has some limitations. Although we discovered unique gene profiles pertaining to key genes involved in IFN-I signaling by low-dose radiation, the up-or downregulated expression patterns of some key genes were inconsistent with the RNA-seq results in T2D-HAECs. Thus, it is our assumption that the pathophysiological features of diabetes may be more potent than the biological effects exerted by low-dose radiation. Therefore, further studies are required to uncover the biological mechanisms active in normal conditions as well as pathological conditions. Furthermore, we did not perform in vivo experiments to clarify the roles of key molecules via the modulation of gene or protein levels and their corresponding physiological features. Moreover, cultured cells are known to undergo genetic variations, and these modifications may alter the responses to irradiation, which would not be evident in primary tissues. Thus, future studies are required to further validate the available data.
In this study, we investigated the effect of low-dose radiation on EC gene profiles in HAECs derived from healthy or T2DM donors. We demonstrated that low-dose radiation could reduce the proliferative rate, increase cellular senescence, and reduce tube formation in normal HAECs, which is consistent with results from cells treated with 2 Gy high-dose radiation. Additionally, low-dose radiation induced changes in EC phenotypes compared to those observed in non-irradiated diabetic conditions, indicating impaired EC function. Moreover, we uncovered key molecules involved in the IFN-I signaling pathway that are associated with CVD in cultured primary ECs under both normal and diabetic conditions. Therefore, our study suggests a potential effect of low-dose radiation on vascular physiology in both normal and diseased states.

MTT Assay
A 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) assay was used to measure cell proliferation after low-dose radiation exposure. Cells were seeded at 5 × 10 3 cells/well in 4-well cell culture dishes (SPL Life Sciences, Gyeonggi-do, Korea) and incubated at 37 • C and 5% CO 2 . After 24 h, cells were irradiated with 100 mGy and 2 Gy and then incubated for another 3 d. A total of 5 mg MTT powder (Duchefa Biochemie, Haarlem, The Netherlands) was dissolved in 1 mL Dulbecco's phosphatebuffered saline (DPBS, Welgene, Gyeongsangbuk-do, Korea) and then mixed into the cell culture medium in a 1:10 ratio, followed by incubation for 2 h. The MTT-supplemented culture medium was removed, and dimethyl sulfoxide (DMSO, Sigma-Aldrich, St. Louis, MO, USA) was added to the cultures, followed by incubation for 30 min. A 100 µL aliquot of the solution was pipetted into 96-well cell culture plates (SPL Life Science) and measured at a wavelength of 570 nm using a Multiskan™ FC microplate photometer (Thermo Fisher Scientific, Waltham, MA, USA)

Tube Formation Assay
Matrigel (354230, Corning, MA, USA) was coated onto 24-well plates at 37 • C for 2 h before cell detachment. HAECs and T2D-HAECs were irradiated with 0, 100 mGy, and 2 Gy and incubated for an additional 24 h. The next day, cells were detached, seeded into culture plates at 1.5 × 10 4 cells, and incubated for 16 h at 37 • C. Tube formation was analyzed and quantified using a Nikon ECLIPSE microscope (Nikon, Tokyo, Japan) and its associated software.

RNA Isolation
Total RNA was isolated using TRIsure solution (Bioline, London, UK) following the manufacturer's instructions. RNA quality (expressed as an RNA integrity number) was assessed with an Agilent 2100 bioanalyzer using an RNA 6000 Nano Chip (Agilent Technologies, Amstelveen, The Netherlands). Total RNA was quantified using a NanoDrop 2000 spectrophotometer (ND-2000; Thermo Fisher Scientific, Inc.).

RNA-seq
RNA-seq was performed with high-quality RNA samples (RNA integrity number > 7) isolated from each cell type. Four separate samples were multiplexed into each lane and sequenced on a HiSeq 4000 system (Illumina, San Diego, CA, USA).

Identification of DEGs and Data Analysis
The RNA-seq reads were aligned using Bowtie2 [53]. Bowtie2 indices were generated from the genome assembly sequences or the representative transcript sequences for alignment to the genome. The alignment file was used to assemble transcripts, estimate abundance, and detect DEGs. Data mining and graphic visualization were performed using the Excel-based Differentially Expressed Gene Analysis (ExDEGA; Ebiogen Inc., Seoul, Korea). Probe sets without corresponding gene symbols were removed, and sand differences with a p-value < 0.05 and absolute fold change > 1.2 were considered statistically significant.

Heatmap Visualization and Hierarchical Clustering Analysis
The expressions of DEGs were visualized as a heatmap using the open-source program MeV (version 4.9.0; https://sourceforge.net/projects/mev-tm4/, accessed on 22 November 2021). Gene expression levels of each sample were indicated by different colors using log 2 ratio adjusted z-score values, and hierarchical clustering of gene and sample trees was analyzed using the Euclidean distance metric.

Pathway Analysis
DEG lists were annotated to KEGG pathways to identify significant pathways. KEGG pathways were analyzed using GSEA 4.0.3 software. Pathways with an FDR-adjusted p-value < 0.05 were filtered, and a Bonferroni correction was applied to adjust the p-values. Pathways were considered significantly enriched if the p-value < 0.05 and FDR q-value < 0.25.

Gene Ontology (GO) Analysis
All DEGs were analyzed by the GO database (http://www.geneontology.org/, accessed on 29 November 2021). Biological processes with an FDR threshold <0.05 were filtered, and biological process-related categories were selected and grouped hierarchically. Cytoscape's ClueGO plug-in was used to analyze the functional groups of GO terms related to biological processes. Significant interrelated GO terms were determined as adjusted p-value < 0.05.

Protein-Protein Interaction Network and Module Analysis
Protein-protein interaction (PPI) networks were mapped and analyzed using Cytoscape (version 3.9.1; https://cytoscape.org/, Gladstone, San Francisco, CA, USA, 22 November 2021). Significant modules in the PPI networks were identified using Molecular Complex Detection (MCODE), a plug-in app of Cytoscape designed to analyze densely connected regions by clustering a given network.
4.14. Total RNA Extraction and Quantitative Reverse Transcription-Polymerase Chain Reaction (qRT-PCR) Total RNA was extracted using the TRIsure solution and reverse-transcribed into cDNA using the SensiFAST TM cDNA Synthesis kit (Bioline) from the isolated RNA. PCR was performed with a Mic Real-Time PCR system (Bio Molecular Systems, QLD, Australia), the Qualitative PCR SYBR 2X Master Mix Kit (Mbiotech), and primer pairs specific for target genes (Table 4). Table 4. Primers used in qRT-PCR.

Gene Symbol
Forward Sequence Reverse Sequence CTGTGCCATGGAGAGTAGCA AGGTGGATTGTCAGGGTCTG

Statistical Analysis
Data are expressed as mean ± standard deviation (SD). One-way ANOVA followed by the Tukey post hoc analysis or Student's unpaired two-tailed t-test were performed using the GraphPad Prism program (version 8.0, GraphPad Software, San Diego, CA, USA) and are indicated in the figure legends. Values with p-values < 0.05 were considered statistically significant.