Screening Plasma Exosomal RNAs as Diagnostic Markers for Cervical Cancer: An Analysis of Patients Who Underwent Primary Chemoradiotherapy

This preliminary study aimed to screen non-coding RNAs (ncRNAs) from plasma exosomes as a new method for cervical cancer diagnosis. Differentially expressed RNAs were initially selected from among a group of 12 healthy individuals (normal group) and a pretreatment group of 30 patients with cervical cancer (cancer group). Then, we analyzed the association between an ncRNA-mRNA network and cancer using ingenuity pathway analysis after secondary selection according to the number and correlation of mRNAs (or ncRNAs) relative to changes in the expression of primarily selected ncRNAs (or mRNAs) before and after chemoradiotherapy. The number of RNAs selected from the initial RNAs was one from 13 miRNAs, four from 42 piRNAs, four from 28 lncRNAs, nine from 18 snoRNAs, 10 from 76 snRNAs, nine from 474 tRNAs, nine from 64 yRNAs, and five from 67 mRNAs. The combination of miRNA (miR-142-3p), mRNAs (CXCL5, KIF2A, RGS18, APL6IP5, and DAPP1), and snoRNAs (SNORD17, SCARNA12, SNORA6, SNORA12, SCRNA1, SNORD97, SNORD62, and SNORD38A) clearly distinguished the normal samples from the cancer group samples. We present a method for efficiently screening eight classes of RNAs isolated from exosomes for cervical cancer diagnosis using mRNAs (or ncRNAs) altered by chemoradiotherapy.


Introduction
The age-standardized incidence rate (ASR) in 100,000 for cervical cancer in South Korea was 8.7% in 2017, which was a significant improvement over the rate of 18.6% in 1999 [1,2]. This decrease in ASR is thought to be related to the establishment of regular health check-ups and dissemination of human papillary virus vaccines. However, unlike the overall ASR improvement, the ASR of 25-29-year-olds gradually increased from 3.6% in 2000 to 6.5% in 2011 [2]. The incidence of cervical cancer has decreased in 60-70-year-old individuals but increased in those aged between 20-30 years, and this age-dependent incidence has become clearer as the level of national development increases [3]. This pattern may be attributed to the active participation of young women in the national health check-up program, the high intake of junk food and smoking habits in young women, and sexual intercourse at younger ages [4,5]. In Korea, a cervical cytology test is currently performed every three years for women over the age of 20 years [6]. However, not only is the pap smear as a screening test a potential source of fear and shame, its sensitivity can also be as low as 50% [7]. Although this limitation can be overcome by the high sensitivity of an HPV test, 12.7% of squamous cell carcinoma (SCC) and 15-38% of adenocarcinoma (AC) globally are not related to HPV [8]. In addition, the sensitivity of the pap smear test can further decrease depending on the skill level of the physician and the examination

Blood Samples and Clinical Data
Two sets of 5-10 mL blood samples were collected from 30 patients diagnosed with stage IB-IVB cervical cancer and treated with CCRT at the Department of Radiation Oncology at Ajou University from June 2018 to March 2020. These samples were stored at the Biobank at Ajou University Hospital, which is a member of the Korea Biobank Network. Sixty samples in total were acquired before treatment and after the second week of CCRT. Blood plasma samples (3 mL) from 12 healthy individuals matched for sex and age were obtained from the Biobank at Ajou University Hospital (institutional review board approval number: BMR-EXP-20-428). Plasma exosomal RNA sequencing was conducted by Macrogen (www.macrogen.com, Supplementary Methods). The diagnosis of all patients was histologically confirmed via biopsy. Regional lymph node metastasis and distant metastasis were evaluated using magnetic resonance imaging and positron emission tomographycomputed tomography. Detailed treatment and follow-up procedures were previously reported for most of the current patients [17]. Clinical data, such as information regarding age, 2018 International Federation of Gynecology and Obstetrics stage, pathology, radiation therapy (RT) field, pretreatment hemoglobin levels, pretreatment absolute lymphocyte count (ALC), ALC after two weeks of CCRT (ALC2), and the levels of pretreatment tumor markers SCC antigen and cytokeratin fragment 21-1 were collected retrospectively from electronic medical records.

Screening Process
The screening process of seven classes of ncRNAs and mRNAs from plasma exosomes for cervical cancer diagnosis is presented in Figure 1. The process was divided into statistical and biological screening phases.
Biomolecules 2021, 11, x FOR PEER REVIEW 3 of 16 pretreatment tumor markers SCC antigen and cytokeratin fragment 21-1 were collected retrospectively from electronic medical records.

Screening Process
The screening process of seven classes of ncRNAs and mRNAs from plasma exosomes for cervical cancer diagnosis is presented in Figure 1. The process was divided into statistical and biological screening phases. Figure 1. A plasma exosomal RNA screening method for cancer diagnosis that consists of a statistical screening phase followed by a biological screening phase.

Statistical Screening
The expression of plasma exosomal ncRNAs and mRNA was compared between a group of 12 samples from healthy individuals (normal group) and a group of 30 samples from patients with cervical cancer before CCRT (cancer group). RNAs with both |log2FC| > 2 and p-values < 0.05 were selected from DEG analysis (A). The p-value of (A) was defined as the "DEG p-value". The expression of plasma exosomal ncRNAs and mRNA was compared between the normal group and a group of 24 samples from patients with cervical SCC or unclassified carcinoma before treatment (non-AC group). RNAs with |log2FC| >1.5 and p-values < 0.05 were selected from DEG analysis (B).
The expression of plasma exosomal ncRNAs and mRNA was compared between the normal group and a group of six samples from patients with cervical adenocarcinoma or adeno-squamous cell carcinoma before treatment (AC group). RNAs with |log2FC| >1.5 and p-values < 0.05 were selected from DEG analysis (C). The significant DEGs were primarily selected when the RNAs of (A) were simultaneously included in those with |log2FC (B) + log2FC (C)| > 4. These results were visualized using volcano plots.

Biological Screening
We calculated the log2FC values of plasma exosomal RNAs from 30 patient samples taken two weeks after CCRT, which were compared with those collected before CCRT to screen for secondary biological functions of the initially selected DEGs. If the number of Figure 1. A plasma exosomal RNA screening method for cancer diagnosis that consists of a statistical screening phase followed by a biological screening phase.

Statistical Screening
The expression of plasma exosomal ncRNAs and mRNA was compared between a group of 12 samples from healthy individuals (normal group) and a group of 30 samples from patients with cervical cancer before CCRT (cancer group). RNAs with both |log 2 FC| > 2 and p-values < 0.05 were selected from DEG analysis (A). The p-value of (A) was defined as the "DEG p-value". The expression of plasma exosomal ncRNAs and mRNA was compared between the normal group and a group of 24 samples from patients with cervical SCC or unclassified carcinoma before treatment (non-AC group). RNAs with |log 2 FC| >1.5 and p-values < 0.05 were selected from DEG analysis (B).
The expression of plasma exosomal ncRNAs and mRNA was compared between the normal group and a group of six samples from patients with cervical adenocarcinoma or adeno-squamous cell carcinoma before treatment (AC group). RNAs with |log 2 FC| >1.5 and p-values < 0.05 were selected from DEG analysis (C). The significant DEGs were primarily selected when the RNAs of (A) were simultaneously included in those with |log 2 FC (B) + log 2 FC (C)| > 4. These results were visualized using volcano plots.

Biological Screening
We calculated the log 2 FC values of plasma exosomal RNAs from 30 patient samples taken two weeks after CCRT, which were compared with those collected before CCRT to screen for secondary biological functions of the initially selected DEGs. If the number of DEGs was less than 20 for miRNA and snoRNA, we analyzed the extent to which the network of mRNAs relative to the log 2 FC DEG(s) was associated with cancer category as described by IPA, and calculated the percentage of RNAs relative to the cancer category in the network. If the number of DEGs was more than 20, the DEGs were divided into two groups to maximize the difference in "DEG p-values" according to number of mRNAs associated with the log 2 FCs in the selected ncRNAs (e.g., lncRNA, piRNA, snRNA, tRNA, and yRNA). When a class of selected RNAs was mRNA, we used the number of miRNAs, piRNAs, and lncRNAs relative to the identified DEGs. We performed a network analysis and IPA after screening the DEGs satisfied with many mRNAs or ncRNAs, and high −log 10 (DEG p-value) based on this grouping.

Statistical Analysis of Differential RNA Expression
Raw data (i.e., the reads for each RNA) were normalized by TMM (trimmed mean of M value) using edgeR. For pre-processing, the RNAs undetected in over 12 of 42 samples and those undetected in over 30 of 60 samples were filtered during DEG analysis for statistical and biological screening, respectively.

Multidimensional Scaling and Heatmap Construction
Classical multidimensional scaling (MDS) of all 42 samples calculated by the cmdscale function was visualized by scattered data partitioned into two or three groups using k-means clustering. A hierarchical clustering heatmap was plotted using the pheatmap function. In MDS and heatmap analysis, log 2 count per million (CPM) values from normalized read counts were used for each RNA class from the statistical screening, while those from raw read counts were applied to optimally integrate the different classes of RNAs that resulted from the screening process.

Network Analysis
We formed a network using Prim's algorithm for the minimum spanning tree in the igraph package for R. Edges are presented in red and blue, which indicate positive and negative correlations, respectively.

Ingenuity Pathway Analysis
The most significant diseases and bio-functions were analyzed using the IPA software (Qiagen, https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis) and selected based on the activated Z-scores of the downstream effects of the analysis [18]. Positive and negative Z-scores indicated the induced and inhibited functional activities, respectively.

Receiver Operative Characteristic Analysis
We visualized receiver operative characteristic (ROC) curves and calculated the area under the curves, sensitivities, and specificities using the Epi and pROC packages.

Table and Boxplots
Continuous variables between groups were compared using the t-test or Wilcoxon rank-sum test according to a Shapiro-Wilk test, while categorical variables were compared using Fisher's exact test or the chi-square test (Table 1). Statistical analysis for logarithmic values of reads per million mapped reads (RPM) or reads per kilobase per million mapped reads (RPKM) between groups was performed using the Wilcoxon rank-sum test or Kruskal-Wallis test in all boxplots. Data analysis and visualization were performed using R version 4.1 (https://www.r-project.org).  Table 1 shows the clinical characteristics of the cancer group (n = 30). Patients with cervical cancer confined within the pelvis (stage IB-IIIC1) made up 66.7% of the cohort, while those that had extra-pelvic metastases made up 33.3% of the group. Of the 24 non-AC patients, 23 (95.8%) had SCC pathology.

Statistical Screening
We selected 13 miRNAs, 43 piRNAs, 28 lncRNAs, and 67 mRNAs from three DEG analyses, as shown in Figure 2A,D,G,J (normal vs. cancer group, normal vs. non-AC group, and normal vs. AC group). The MDS scatter plots and k-means clustering of the selected RNAs divided 42 samples into two groups that had both sensitivity and specificity for cancer diagnosis in each of the four RNA classes ( Figure 2B,E,H,K). The sensitivity and specificity of miRNAs were 83% (predicted 25 of 30 patients) and 100%, piRNAs were 87% (predicted 26 of 30 patients) and 92% (predicted 11 of 12 healthy individuals), lncRNAs were 90% (predicted 27 of 30 patients) and 75% (predicted 9 of 12 healthy individuals), and mRNAs were 90% (predicted 27 of 30 patients) and 83% (predicted 10 of 12 healthy individuals), respectively. Heatmaps of DEG expression levels for the four groups showed differences between the normal and cancer groups, but did not display differences in disease stage or pathology ( Figure 2C,F,I,L). A total of 18 snoRNAs, 76 snRNAs, 474 tRNAs, and 64 yRNAs were selected from three DEG analyses, as presented in Figure 3A,D,G,J (normal vs. cancer group, normal vs. non-AC group, and normal vs. AC group). The MDS scatter plot and k-means clustering of selected snoRNAs identified three groups with a sensitivity of 70% (predicted 21 of 30 patients) and specificity of 100% for cancer diagnosis ( Figure 3B). In Figure 3C, the heatmap of snoRNA expression in 42 samples shows two groups of snoRNAs that could distinguish between the normal and cancer groups. A group of snoRNAs ( Figure 3C; red dotted line) was homogeneously overexpressed in 12 patients within the cancer group, A total of 18 snoRNAs, 76 snRNAs, 474 tRNAs, and 64 yRNAs were selected from three DEG analyses, as presented in Figure 3A,D,G,J (normal vs. cancer group, normal vs. non-AC group, and normal vs. AC group). The MDS scatter plot and k-means clustering of selected snoRNAs identified three groups with a sensitivity of 70% (predicted 21 of 30 patients) and specificity of 100% for cancer diagnosis ( Figure 3B). In Figure  3C, the heatmap of snoRNA expression in 42 samples shows two groups of snoRNAs that could distinguish between the normal and cancer groups. A group of snoRNAs ( Figure 3C; red dotted line) was homogeneously overexpressed in 12 patients within the cancer group, while a group of snoRNAs ( Figure 3C; red solid line) was nonhomogeneously overexpressed in the other 12 patients within the cancer group. We defined the former group as the DEG in snoRNA. MDS scatter plots and k-means clustering using selected snRNAs, tRNAs, and yRNAs divided 42 samples into two groups with a sensitivity of 40% (predicted 12 of 30 patients) and specificity of 100% for cancer diagnosis for all three types of RNAs ( Figure 3E,H,K). The heatmaps for these three groups showed that the 12 patients with DEG in snoRNA were classified into the same group in a column hierarchical clustering through the expression of selected RNAs ( Figure 3F,I,L). while a group of snoRNAs ( Figure 3C; red solid line) was nonhomogeneously overexpressed in the other 12 patients within the cancer group. We defined the former group as the DEG in snoRNA. MDS scatter plots and k-means clustering using selected snRNAs, tRNAs, and yRNAs divided 42 samples into two groups with a sensitivity of 40% (predicted 12 of 30 patients) and specificity of 100% for cancer diagnosis for all three types of RNAs ( Figure 3E,H,K). The heatmaps for these three groups showed that the 12 patients with DEG in snoRNA were classified into the same group in a column hierarchical clustering through the expression of selected RNAs ( Figure 3F,I,L).   The number and Pearson's correlation of mRNAs associated with the log 2 FCs of 13 miRNAs identified after CCRT are presented in Figure 4A. There were significantly more mRNAs whose expression was altered by hsa-miR-142-3p than the other 12 miRNAs studied. A network with 139 mRNAs affected by hsa-miR-142-3p (R > 0.8) was more relative to the cancer category than the 28 mRNAs that were affected by has-miR-4306 and the other seven miRNAs (R > 0.7). The log 2 (RPM+1) values of hsa-miR-142-3p were significantly lower in the cancer group than in the normal group regardless of pathology, while there was no significant change according to disease stage (Supplementary Figure S3A).

lncRNA
The number and Pearson's correlation of mRNAs associated with log 2 FCs of 28 lncRNAs after CCRT were sorted in descending order ( Figure 4D). A cut-off value of 100 was used to maximize the −log 10 (DEG p-value) between the two groups of 28 lncRNAs according to the number of related mRNAs (R > 0.7; Supplementary Figure S2A). A group of lncRNAs with related mRNAs > 100 had a significantly higher −log 10 (DEG p-value) than the group with related mRNAs ≤ 100 ( Figure 4E). This meant that the number of mRNAs relative to lncRNAs expressed after CCRT was positively correlated with the degree of statistical significance between the normal and cancer groups. Four lncRNAs with related mRNAs > 100 were selected owing to the presence of mRNAs with R > 0.9. A network indicated 76 mRNAs altered by LINC0089 and the other three lncRNAs that were relative to the cancer category ( Figure 4F,G). The log2 (RPKM+1) values for these four lncRNAs were significantly lower in the cancer group than in the normal group regardless of pathology, while there was no significant change according to disease stage, except for LOC105374768 (Supplementary Figure S3B).

mRNA
The number and Pearson's correlation of miRNAs, piRNAs, and lncRNAs associated with the log 2 FCs of 67 mRNAs expressed after CCRT were sorted in descending order ( Figure 4H). A cut-off value of 10 was used to maximize the −log 10 (DEG p-value) between two groups of 67 mRNAs according to the number of related ncRNAs (R > 0.7; Supplementary Figure S2B). A group of mRNAs with related ncRNAs > 10 had a significantly higher −log 10 (DEG p-value) than the group of mRNAs with related ncRNAs ≤ 10 ( Figure 4I). This meant that the number of ncRNAs that were relative to the mRNAs expressed after CCRT was positively related to statistical significance between the normal and cancer groups. Five mRNAs with related ncRNAs > 10 were selected because of the presence of ncRNAs with R > 0.9. A network constructed from five lncRNAs and one miRNA that were affected by five mRNAs showed a significant association with the cancer category ( Figure 4J,K). The log2 (RPKM+1) values of five mRNAs were significantly lower in the cancer group than in the normal group regardless of pathology, while there was no significant change according to disease stage, except for CXCL5 (Supplementary Figure S3C).

snoRNA
The number and Pearson's correlation of mRNAs associated with log 2 FCs for 18 snoR-NAs after CCRT are presented in Figure 4L. There were significantly more mRNAs affected by URS0000822206 or URS000067E6DC than by the other nine snoRNAs. A network of 207 mRNAs affected by URS000002084A and eight other snoRNAs (R > 0.7) was closely related to the cancer category, while a network of 13 mRNAs that were changed by URS0000822206 and URS000067E6DC (R > 0.9) was not ( Figure 4M,N). The log2 (RPM+1) values of four snoRNAs that were not included in the DEG in snoRNA group were significantly higher in the cancer group than in the normal group regardless of pathology, while there was no significant change according to disease stage. (Supplementary Figure S3D). pathology, while there was no significant change according to disease stage. (Supplementary Figure S3D). 3.3.5. piRNA, snRNA, tRNA, and yRNA The initial selection process of ncRNAs was similar to that of the lncRNAs for four classes of piRNAs, snRNAs, tRNAs, and yRNAs (Supplementary Figures S1, S2C-F, and S3E-H).  The initial selection process of ncRNAs was similar to that of the lncRNAs for four classes of piRNAs, snRNAs, tRNAs, and yRNAs (Supplementary Figures S1, S2C-F and S3E-H).

DEG in snoRNA
The 12 patients with DEG in snoRNA had significantly lower ALC2 levels and higher pretreatment tumor marker levels than the other 18 patients (Table 1). Of the five snoRNAs included in the DEG in snoRNA group ( Figure 4L, red dotted line), the log 2 (RPM+1) values of URS000002084A (SNORA12) were most relative to ALC2 according to DEG in snoRNA (R 2 = 0.19; Figure 5A). The difference in Z-score according to DEG in snoRNA was interpreted as the secondary reward activation (activation of cells↑, exocytosis↑) following weakened anticancer activity (degranulation of leukocytes↓; Figure 5B).

Discussion
To select plasma exosomal ncRNA (or mRNA) for the diagnosis of cervical cancer, we targeted RNAs that had statistical differences in expression levels between the normal and cancer patient groups. Candidate ncRNAs (or mRNAs) were selected based on the increase in the number of mRNAs (or ncRNAs) altered by log2FC of ncRNAs (or mRNAs) and the Pearson's correlation between mRNAs and ncRNAs. The association between ncRNA-mRNA networks and cancer category was confirmed using IPA. The process for additionally exploring the changes in RNA expression in response to CCRT helped in the efficient selection of many RNAs with statistical significance from various classes. The core of this process was based on the hypothesis that the number of cancer cells can be reduced by irradiation, which alters the expression of cancer-associated RNA that has a stronger biological influence on the size and correlation of the ncRNA-mRNA network.
When the RNAs present in exosomes were selected and integrated for all classes, they

Integration
The snRNAs, tRNAs, and yRNAs were excluded from integration owing to their low sensitivity in the MDS plots ( Figure 3E,H,K); however, five of the selected nine snoRNAs were DEGs between the two groups based on the DEG in snoRNA ( Figure 5C,D, red dotted line). First, the miRNAs, piRNAs, lncRNAs, mRNAs, and snoRNAs selected via the two screening process steps were visualized as a hierarchical clustering heatmap ( Figure 5C). This heatmap was divided into three groups, a group of 10 patients with increased snoRNAs (two red lines), a group of 10 patients with decreased miRNA, piRNAs, lncRNAs, and mR-NAs (blue line), and a group of 10 patients with both increased snoRNAs and four classes of decreased RNAs (purple line). Second, we selected a combination of miRNAs, mRNAs, and snoRNAs to distinguish the normal and cancer groups after MDS and hierarchical clustering heatmap analyses with various combinations that included snoRNAs ( Figure 5D,E). Third, we selected RGS18, SNORA12, and URS00003B57B1 (SNORD97), which were visualized as an MDS plot and a hierarchical clustering heatmap ( Figure 5F,G). In mRNA and lncRNA networks, RGS18 was interrelated with LOC105374768 and LINC00989, which were selected through the lncRNA screening process (Figure 4F,J). We selected SNORA12, which was at the center of the snoRNA network ( Figure 4M, lower) and closely related to ALC2 ( Figure 5A). Excluding the five snoRNAs in the DEG in snoRNA group, the snoRNA combinations that included SNORD97 were better, than those that included the other three snoRNAs, at distinguishing between the normal and cancer groups in the MDS or heatmap analyses. The ROC curves for cancer diagnosis using log 2 (RPM+1) or log 2 (RPKM+1) showed that RGS18, SNORA12, and SNORD97 had higher area under the curve (AUC), sensitivity, and specificity than RGS18 alone or both SNORA12 and SNORD97 (0.992, 96.7%, and 100% for RGS18, SNORA12, and SNORD97; 0.964, 96.7%, and 91.7% for RGS18; and 0.883, 73.3%, and 100% for SNORA12 and SNORD97; Figure 5H).

Discussion
To select plasma exosomal ncRNA (or mRNA) for the diagnosis of cervical cancer, we targeted RNAs that had statistical differences in expression levels between the normal and cancer patient groups. Candidate ncRNAs (or mRNAs) were selected based on the increase in the number of mRNAs (or ncRNAs) altered by log 2 FC of ncRNAs (or mRNAs) and the Pearson's correlation between mRNAs and ncRNAs. The association between ncRNA-mRNA networks and cancer category was confirmed using IPA. The process for additionally exploring the changes in RNA expression in response to CCRT helped in the efficient selection of many RNAs with statistical significance from various classes. The core of this process was based on the hypothesis that the number of cancer cells can be reduced by irradiation, which alters the expression of cancer-associated RNA that has a stronger biological influence on the size and correlation of the ncRNA-mRNA network.
When the RNAs present in exosomes were selected and integrated for all classes, they were broadly divided into four classes of miRNAs, piRNAs, lncRNAs, and mRNAs and a single class of snoRNAs ( Figure 5C). In particular, the combination of miRNA-mRNA-snoRNA clearly distinguished the normal group from the cancer group ( Figure 5D,E). This can imply various characteristics, including growth, invasion, metastasis, neovascularization, evasion of tumor suppression, genetic instability, inflammation, immune evasion, and alteration of energy metabolism, which are acquired by the tumor microenvironment (TME), as shown in Table 2 [19]. Table 2. Suggested biological functions of selected plasma exosomal RNAs.

CXCL5
Recruits and activates granulocytes and promotes angiogenesis, tumor growth, and metastasis in the tumor microenvironment [23] ↑(CC) [24,25] Tumors with exosome-derived CXCL5 use it to facilitate their progression through infiltration of leukocytes in the tumor microenvironment ↓ Table 2. Cont.

KIF2A
Required for cell mitosis [26] ↑(CC) [27] Rapid mitosis of cancer cells may promote the absorption of KIF2A from exosomes ↓ RGS18 Negative regulator of G protein-coupled receptors and controls platelet activation and production [28,29] ↑(OC) [30] Tumors may absorb RGS18 present in exosomes, which can promote thrombogenesis. The reduction of exosomal RGS18 by tumors may promote activated platelets around the primary tumor, which can facilitate tumor growth and invasion. Therefore, dysregulation of RGS18 can result in tumorigenesis through persistent platelet activation ↓ DAPP1

Activation of antigen-specific T cells [31] NA
This may contribute to tumorigenesis through deficiency of tumor-specific immunity ↓

LINC00989
Decreases with RGS18 in tumor-educated platelets [32] ↓(PaC) [32] The two lncRNAs may facilitate platelet activation in cancer patients via targeting RGS18 The derived RNA positively correlates with CD8 T cell infiltration in thymoma and stomach cancer [33] ↑(COC) [34] Promotion of these snoRNAs present in exosomes may be related to cancer related-lymphopenia The derived RNA negatively correlates with CD8 T cell infiltration in LGG, PC, pancreatic cancer, and HNC [33] ↑(PC) [36] ↑ SNORA12 NA ↓(CC) [37] ↑(LC) [35] ↑ SCARNA1 NA ↑(LC) [35] ↑ SNORD97 NA ↓(CC) [37] Promotion of these snoRNAs present in exosomes may be related to decreased lymphocyte activity The derived RNA negatively correlates with CD8 T cell infiltration in HNC, LC, TGCT, and PCPG [33] ↑(COC) [34]  Previous studies have suggested that expression of miR-142-3p, which suppresses tumor proliferation and invasion, is not only reduced in cervical cancer tissues, compared with that in adjacent normal cervical tissue, but is also associated with prognosis [20,21]. Cunha et al. [22] reported that ARL6IP5, which was one of the selected five mRNAs and included in the tumor suppressor gene database (https://bioinfo.uth.edu/TSGene/, accessed on 1 September 2021), was downregulated more in the soft tissue sarcoma tissues than in benign tumors. The reduced expression of miR-142-3p and ARL6IP5 in plasma exosomes from the cancer group may indicate their role as tumor suppressors during cervical tumorigenesis. The expression of CXCL5 and KIF2A was upregulated in various cancer tissues, including cervical cancer [24,25,27,38]; additionally, the expression of RGS18 in the tissues of patients with ovarian cancer was also upregulated [30]. However, each of the three mRNAs can contribute to tumorigenesis through a different biological mechanism. CXCL5 from tumor cells can promote tumor growth and angiogenesis by recruiting tumor and immune cells into the TME [23]. We might suggest a negative correlation between CXCL5 and cervical cancer stage if CXCL5 from blood plasma exosomes is absorbed and secreted into the TME by tumor cells. The expression level of KIF2A in exosomes could be lower in the cancer group than in the normal group if tumors absorbed KIF2A required for cell mitosis and proliferation [26]. RGS18 deficiency can reportedly lead to platelet activation and reduce platelet survival, whereas its abundance can inhibit platelet activation and promote platelet production [28,29]. The decrease in plasma exosomal RGS18 in the cancer group may be attributed to the absorption of RGS18 by cervical cancer cells. This may result in the promotion of tumorigenesis and tumor progression by facilitating tumorassociated platelet activation, and the presence of RGS18 in the tumor may contribute to sustained blood platelet concentration through its transfer to hematopoietic stem cells. According to a preclinical study, DAPP1 deficiency could not efficiently activate antigenspecific T cells [31], suggesting that DAPP1 downregulation is associated with a weakness in cancer-specific immunity. Four mRNAs (KIF2A, RGS18, ARL6IP5, and DAPP1) and two lncRNAs (LINC00989 and LOC105374768) were included in both the mRNA-ncRNA network and lncRNA-mRNA network ( Figure 4F,J). This study showed that the log 2 FC of RGS18 with CCRT strongly correlated with two lncRNAs in the network. A previous study on the RNA expression of tumor-educated platelets (TEPs) revealed that the expression level of both RGS18 and LINC00989 was lower in platelets from cancer patients than in those from healthy donors [32]. This supports that the two lncRNAs may be associated with tumorigenesis through RGS18-mediated platelet activation. In addition, TEPs may secrete CXCL5, which can promote the recruitment of leukocytes into the TME [39]. Therefore, plasma exosomal RGS18 may be more associated with cervical cancer development and progression than the other three mRNAs (KIF2A, ARL6IP5, and DAPP1).
The nine selected snoRNAs were different from other classes of selected RNAs in that the number of mRNAs relative to the log 2 FC of the snoRNAs was smaller than that of the mRNAs relative to the change in other selected RNAs following CCRT (Figure 4 and Supplementary Figure S1). Therefore, the two snoRNAs with many relative mRNAs ( Figure 4L) may function as a normal physiological response following CCRT. Unlike other ncRNAs, snoRNA contained seven RNAs that were specifically upregulated in 12 patients with cancer (DEG in snoRNA), and many of their snRNA, tRNA, and yRNA genes were differentially expressed in the other 30 patients (Figure 3). These RNAs were associated with both a decrease in ALC2 levels and an increase in pretreatment tumor marker expression (Table 1) and may contribute to a weakening of immunity that was associated with cancer according to the IPA ( Figure 5B). Previous studies of the changes in snoRNA content in extracellular vesicles of immune cells using immunosuppression or immunostimulatory factors and snoRNA secretions stimulated by inflammation suggested that snoRNA functions as a modulator of inflammation [40][41][42]. In addition, RNAs further processed from snoRNAs, including SNORD17 or SNORA6, among seven snoRNAs (DEG in snoRNA), were associated with the number of tumor-infiltrating CD8 T cells [33]. Therefore, we can assume that cervical cancer may be associated with snoRNA, which would reflect an impaired inflammatory or immune response followed by lymphopenia. Additionally, lymphopenia in patients with cervical cancer was associated with clinical results of CCRT [43,44], indicating that DEG in snoRNA may be related to cervical cancer prognosis. We selected SNORA12 as the most relative to ALC2 among the five selected DEGs in snoRNA ( Figure 5A,D, red dotted line), which showed that snoRNA was more suitable for diagnosing cancer than the other three RNA classes (e.g., snRNA, tRNA, and yRNA) because it had an additional four RNAs that were differentially expressed in the cancer group ( Figure 5D, red solid line). snoRNAs may also be relevant to cancer immunity, considering the association between SNORD38A (URS00003640C3 or URS000067EB9D) of the four snoRNAs and the number of tumor-infiltrating CD8 T cells [33]. Among the four selected snoRNAs, the combination of SNORD97 with both RGS18 and SNORA12 was the most appropriate for distinguishing between the normal and cancer groups using the MDS plot, heatmap, and ROC curve analyses ( Figure 5F,G,H). The downregulation of SNORA12 and SNORD97 in tissues from patients with cervical squamous cell carcinoma [37] implies that their secretion from cancer cells may help to evade host immunity against tumors through suppression of both systemic lymphopenia and lymphocyte activity in the TME. Taken together, our data suggest that the 15 plasma exosomal RNAs from the three classes of miRNA, mRNAs, and snoRNAs distinguish the normal and cancer groups by reflecting the evasion of tumor suppression (miR-142-3p and ARL6IP5), tumor proliferation (KIF2A), tumor progression through the TME (CXCL5 and RGS18), and cancer immunity (DAPP1 and 9 snoRNAs). In particular, tumor-related platelet activation by downregulation of RGS18 and an immune suppression by upregulation of both SNORA12 and SNORD97 may be essential mechanisms for cervical cancer development and progression.
We report a secondary screening method using the log 2 FC of mRNAs (or ncRNAs) before and during CCRT in patients with cervical cancer and the potential diagnostic RNAs selected by this method; however, there are several limitations to these results and require further investigation. First, biases can arise from the process of obtaining and analyzing RNA sequences, such as the heterogeneity of exosomes and their RNAs due to CCRT, the different methods used to isolate exosomes, and the different RNA profiling methods [42]. Therefore, the reproducibility of the present study is not guaranteed. Second, this is a preliminary analysis with a small number of patients, and therefore, the reliability of this study is limited. Third, additional validation of RNAs selected using the proposed screening method is still required, including biological mechanism characterization or animal studies for the selected RNAs presented in this study. We attempted to overcome the bias of the RNA sequencing analytic process by performing paired comparisons of the 60 samples from 30 patients before and during CCRT. Moreover, the 30 samples collected during CCRT were irradiated at a constant dosage per time point (approximately 18 Gy every 2 weeks) and cisplatin regimen (30-70 mg/m 2 every week), as reported in a previous study [17]. However, the reliability of the reported methodology should be validated by applying a similar method to select plasma exosomal RNAs for diagnosis using various sizes of datasets from different cancer patients treated with RT. Furthermore, the 15 RNAs selected from this pilot study should be verified in a large cohort dataset and validated with preclinical studies, despite their diagnostic potential for cervical cancer and implications as hallmarks of cancer.

Conclusions
We present the first method for efficiently screening cancer-related RNAs using mR-NAs (or ncRNAs) relative to the log2FC of ncRNAs (or mRNAs) altered by CCRT and suggest the association between the RNAs identified by this method and their dysregulation in cancer. This method may reduce the time and resources needed to develop diagnostic and therapeutic biomarkers for cervical cancer and should be validated in further preclinical and clinical studies. Informed Consent Statement: Informed consent was waived because the analyzed blood plasma samples from 12 healthy individuals were procured from the institutional Biobank, and the 58 samples from 29 cancer patients that were used for NGS had been analyzed in a previously published article by the authors. Two samples from one cancer patient were obtained with informed consent. Data Availability Statement: All data analyzed during this study are available. Original sequencing data: ArrayExpress (accession number: E-MTAB-10215, E-MTAB-10930) Coding and dataset: https: //data.mendeley.com/datasets/bnzzvk38d8/draft?a=f5ea1853-dcf8-4f74-ba61-a8e527a11a9e.