This multicenter longitudinal study was conducted at three tertiary Grade A hospitals (A, B, and C) located in an NPC endemic area in China. All participants were informed by a doctor or nurse about the purpose and significance of this study, and their rights and obligations, and they all provided informed consent. Patients were enrolled consecutively upon admission between July 2024 and February 2025. Baseline data and an electronic mental health survey were collected during their hospitalization (The patient selects the paper or electronic questionnaire to fill in the questionnaire). Participants completed a questionnaire at each time point and provided blood samples for immune evaluation at admission and after the 30th radiotherapy session.

The inclusion criteria were: i) age ≥18 years; ii) pathological diagnosis of WHO type II-IV NPC within 1 week; iii) having stage I-IVB disease [UICC/AJCC Classification (8th edition)]; iv) NPC patients first receiving radiotherapy (intensity-modulated radiation therapy) with chemotherapy and complete it within 90 days; v) sufficient cognitive and language abilities to finish the electronic questionnaire independently. The exclusion criteria were: i) participating or participated in psychological intervention experimental studies; ii) having heart, brain, lung, liver, kidney, or other organ diseases. This study was carried out following the Declaration of Helsinki and approved by the Ethics Review Committee of Guangxi Medical University (Project No. KY0172, 2024).

The sample size of the longitudinal study was calculated by GLIMMPSE software. Hotelling Lawley Trace (23) was selected as the statistical test for the LCGM, with a significance level of 0.05 and statistical power of 0.90 (β = 0.10). Based on preliminary investigation in this study, NPC patients typically undergo 30 to 33 radiotherapy sessions over a period of 45 to 60 days. As radiotherapy-related adverse effects may accumulate over time, the total number of sessions was evenly divided to assess patients’ PD at different stages, namely at baseline (before treatment), after 15 sessions (mid-treatment), and after 30 sessions (post-treatment).

Repeated measures were specified and four time points were established: at admission (T1), after the first radiotherapy session (T2), after the 15th session (T3), and after the 30th session (T4). The estimated sample size was 190. Considering the potential for attrition or dropout in longitudinal studies, a 20% attrition rate was allowed in this study. Therefore, the required sample size increased to 190/(100%-20%) =238. Therefore, at least 238 NPC patients undergoing radiotherapy should be enrolled in this study. A sampling ratio of 3:1:1 was determined based on the annual admission numbers for NPC patients receiving radiotherapy at the three hospitals. At least, it was expected to collect 136, 46, and 46 cases in hospitals A, B, and C respectively.

Paired-end reads were obtained through sequencing, and quality control was conducted via Q30. After quality control, high-quality reads were identified by filtering raw data with fastp software to remove adapters and low-quality reads. Pruned data were mapped to the reference genome. Expression profiles, differentially expressed transcripts, and differentially expressed genes (DEGs) were calculated. Volcano plots were generated to visualize DEGs. All data processing was carried out by CloudSeq Inc., and results were verified by two researchers to ensure accuracy for subsequent analysis.

Regarding the top 10 differentially expressed circRNAs (DEcircRNAs), their target miRNAs were predicted via RNAhybrid, circBank, CircInteractome, and TargetScan. Overlapping miRNAs were identified by intersecting the four databases and the differentially expressed miRNAs (DEmiRNAs). Target mRNAs for the intersecting miRNAs were predicted via MiRTarBase. A circRNA–miRNA–mRNA (ceRNA) network was built based on the intersecting miRNAs, the top 10 dysregulated DEcircRNAs, and the predicted mRNAs. The ceRNA network was visualized through Cytoscape 3.9.1. Supplementary Table 1 illustrates the top 10 upregulated and downregulated circRNAs.

The random three upregulated circRNAs identified in the ceRNA network were validated via RT-qPCR. Peripheral blood mononuclear cells were isolated from blood samples collected from Class 4 and Class 1 patients at T1 and T4 to extract total RNAs via Trizol (Invitrogen, USA). Gel electrophoresis and Spark 10 M (Tecan, Switzerland) were utilized to appraise RNA quality. Based on the manufacturer’s instructions, reverse transcription was carried out via PrimeScript™ RT Master Mix Kit (Takara, Japan). TB Green^® Premix Ex Taq™ (Takara, Japan) was utilized to amplify cDNA on the ABI StepOnePlus system (Applied Biosystems, USA). CircRNA primers were formed through backsplicing on the CircInteractome and National Center for Biotechnology Information website and were then synthesized by Sangon Biotech (Shanghai, China). The 2^−ΔCT method was utilized to calculate the data, with GAPDH was used as internal reference. Every group consisted of three technical replicates. The specificity of primers was verified by qRT-PCR experiments, and the results showed that the primer amplification efficiency was good and the specificity was strong. Supplementary Table 2 lists all primer sequences.

The biological processes, cellular components, and molecular functions of genes can be explored through GO analysis. KEGG analysis reveals enriched pathways related to human diseases, metabolism, and cellular processes. Therefore, KEGG analysis can support constructing a regulatory network for exploring the development and progression of a disease. GO and KEGG analyses were carried out to reveal the impacts of circRNAs and miRNAs on the rise and decline distress group. The top 10 pathways based on the Enrichment Score (ES) were selected to explore the pathways related to the maternal genes of circRNAs and identified mRNAs in the ceRNA network. R version 4.4.1 was applied to conduct KEGG and GO analyses.

Hub genes in NPC patients with PD were identified to further examine the ceRNA network. The STRING database (https://string-db.org/) was utilized to build a PPI network. Moreover, the PPI network was visualized by the CytoHubba plugin in Cytoscape (version 3.9.1) to obtain the top 10 genes within the ceRNA network. The PPI network included interactions with a combined score greater than 0.7.

Target RBPs for the top 10 dysregulated DEcircRNAs were predicted via circAtlas (28) and RBPsuite databases. The intersecting RBPs were screened by Venny 2.0. Furthermore, the predicted RBPs and dysregulated circRNAs were utilized to build a circRNA–RBP network. Cytoscape 3.9.1 was then applied to visualize the network.

R version 4.4.1 was applied to carry out statistical analysis. Categorical data were presented as absolute counts and percentages, while continuous data were expressed as mean ± standard deviation for normally distributed variables, and median with interquartile ranges for non-normally distributed variables. Chi-square tests were applied to compare categorical data, and independent-sample t-tests were employed to compare normally distributed continuous data between two groups. Analysis of variance was used for multiple group comparisons, while the LSD-t test was applied to pairwise comparisons. Non-normally distributed data between two groups were compared via the non-parametric Wilcoxon test, while those among multiple groups were compared through the Kruskal–Wallis H test.

LCGM was analyzed via Mplus version 8.3. The model was tested with 1 to 5 classes sequentially. The average posterior probability were calculated for each participant in the model. This value represented the posterior probability of each individual being assigned to the corresponding trajectory group. An average posterior probability greater than 70% was considered acceptable. Model fit was evaluated by the sample size-adjusted BIC (aBIC), Bayesian Information Criterion (BIC), and Akaike Information Criterion (AIC). Smaller values indicated a better fit. The bootstrapped likelihood ratio test (BLRT) and Lo-Mendell-Rubin likelihood ratio test (LMR-LRT) were applied to compare k-class models with k−1 class models, with P-values < 0.05 indicating that the k-class model performed notably better. Entropy suggested classification accuracy, with values > 0.8 representing approximately 90% accuracy. A minimum proportion of 5% per class was required. Finally, logistic regression was carried out to explore factors influencing PD in NPC patients. Sociodemographic variables with statistical significance were included as independent variables, and latent classes of PD as dependent variables. A two-sided significance level of α = 0.05 was applied.

300 eligible patients were initially enrolled. During follow-up, nine patients were lost at T2, four at T3, and five at T4. Ultimately, 282 patients were enrolled in the study (170 from Hospital A, 56 from Hospital B, and 56 from Hospital C). Table 1 illustrates the sociodemographic and clinical features of this population.

LCGM on 282 NPC patients undergoing radiotherapy at T1, T2, T3, and T4 was built to identify the trajectories of PD. The fit results for the heterogeneous trajectories of PD are shown in Table 2 . The five-class model did not meet the criteria (P > 0.05 in the likelihood ratio test). The two-, three-, and four-class models fit well as they all had entropy values >0.8, and the LMR-LRT and BLRT reached a significant level (P < 0.05) with each class representing more than 5% of the total. The four-class model had the lowest AIC(4232.533), BIC(4287.108), and aBIC values(4239.544), suggesting that it was the best-fitting model. In this model, 32 patients (11.0%) were assigned to Class 1, 163 (58.0%) to Class 2, 30 (10.7%) to Class 3, and 57 (20.3%) to Class 4. The PD trajectories of these four classes are illustrated in Figure 1 .

Each class was designated according to its trajectory. The DT cutoff score for PD was based on the NCCN recommendations (DT score ≥4). Class 1 (C1, green) patients exhibited relatively high PD levels at T1 and then showed an accelerated decline. Thus, C1 was defined as a decline distress group; Class 2 (C2, blue) patients showed low PD levels at T1 and remained stable across T2–T4. Thus, C2 was defined as a stable distress group; Class 3 (C3, yellow) patients had low PD levels at T1 and remained consistent with DT scores below 4 at T2–T4. Thus, C3 was defined as a no distress group; Class 4 (C4, red) patients exhibited the lowest PD levels at T1 but experienced a sharp increase during T2–T4, peaking at T4. Thus, C4 was defined as a rise distress group.

Among the four trajectories identified through LCGM, we found that when four types of trajectories were selected for single-factor analysis, the only meaningful influencing factor obtained was the PSS score, which would affect the subsequent analysis of influencing factors. Therefore, in order to eliminate potential selection biases, we consulted relevant experts and ultimately the rise, stable and decline distress groups were selected for univariable and multivariable logistic regression analyses to explore potential factors influencing PD trajectories. Univariable analysis ( Table 3 ) indicated that household monthly per capita income, the coping and yielding strategies in medical coping modes, and PSS were significantly associated with the differences among the three PD trajectories in NPC patients undergoing radiotherapy (P < 0.05).

A binary logistic regression analysis was performed using the rise, stable, and decline distress groups as dependent variables, while household income, coping and yielding strategies, and PSS were considered as independent variables ( Table 4 ). Firstly, to comprehensively compare different PD trajectories, the rise distress group was used as the reference group. The results showed that, compared to the stable distress group, patients with a monthly income ≤ 1000 were more likely to be assigned to the rise distress group (OR = 0.448). Patients with higher coping strategy scores had a greater likelihood of being assigned to the stable distress group (OR = 1.535). Compared with the decline distress group, patients with a monthly income ≤ 1000 (OR = 0.316), a monthly income of 1001-3000 (OR = 0.275), and a higher score for yielding strategies (OR = 0.316) were more likely to be categorized into the rise distress group. Patients with higher coping strategy scores (OR = 1.825) and increased PSS (OR = 1.530) had a greater likelihood of being classified into the decline distress group. Furthermore, when the stable distress group was used as the reference, the results indicated that higher coping strategy scores (OR = 1.734) and increased PSS (OR = 1.108) were associated with a greater likelihood of being categorized in the decline distress group.

After identifying four trajectories of PD, the genetic factors influencing the upward or downward PD changes during the radiotherapy process in NPC patients were further investigated. Therefore, three patients from Class 1 (decline distress group) and three patients from Class 4 (rise distress group) were randomly selected. WTS was then conducted based on the 12 blood samples collected at T1 and T4. 34444 circRNAs were obtained by intersecting the sequencing data collected from the Class 1 and Class 4 patients at T1 and T4. DEcircRNAs were identified based on the criteria as follows: adjusted p ≤ 0.05, log (FC) ≥ 1.0, and average counts per million ≥ 100. Finally, 600 DEcircRNAs were identified, with 189 being upregulated and 411 downregulated between the two groups ( Supplementary Figure 1 ). As illustrated in Supplementary Figure 1 , the expression patterns of circRNA between Class 4 patients and Class 1 patients could be clearly distinguished.

1257 miRNAs were identified by intersecting miRNA-sequencing data collected from three Class 1 and three Class 4 patients at T1 and T4. The screening criteria were as follows: log (FC) ≥ 1.0 and adjusted p ≤ 0.05. 123 DEmiRNAs were selected, with 61 significantly upregulated and 62 significantly downregulated between the two groups. Hierarchical clustering demonstrated that the miRNA expression profiles could be effectively distinguished between Class 4 and Class 1 patients ( Supplementary Figure 2 ).

The upregulated and downregulated genes were used to carry out KEGG and GO analyses to investigate the functions of the maternal genes of the DEcircRNAs in NPC patients with PD. KEGG analysis unraveled that the upregulated genes were enriched in the polycomb repressive complex, neurotrophin signaling pathway, and MAPK signaling pathway. GO analysis demonstrated that the upregulated genes were enriched in positive regulation of inflammasome-mediated signaling pathways, chromosome segregation, and positive regulation of NLRP3 inflammasome complex assembly. Additionally, KEGG analysis demonstrated that the downregulated genes were enriched in the viral life cycle, RNA degradation, and the polycomb repressive complex. While GO analysis indicated that the downregulated genes were enriched in lymphocyte differentiation, mononuclear cell differentiation, and positive regulation of post-translational protein modification ( Supplementary Figure 3 ).

A ceRNA network based on the top 10 dysregulated DEcircRNAs was built to better understand the impact of the ceRNA network on the dynamic changes in gene expression regulation linked to PD in NPC patients. For the top 10 upregulated DEcircRNAs, circBank, TargetScan, CircInteractome, and RNAhybrid were applied to predict 73 overlapping miRNAs (four circRNAs without predicted miRNAs were removed). Nine miRNAs were obtained by intersecting the predicted miRNAs and DEmiRNAs. 1062 target mRNAs for these miRNAs were identified via MiRTarBase. 6 circRNAs (including hsa_circ_0004277, hsa_circ_0003684, hsa_circ_0000721, hsa_circ_0106601, hsa_circ_0023249, hsa_circ_0031814), 9 miRNAs (including hsa-miR-1287-5p, hsa-miR-1343-3p, hsa-miR-148b-5p, hsa-miR-181a-3p, hsa-miR-191-5p, hsa-miR-193b-3p, hsa-miR-330-5p, hsa-miR-582-3p, hsa-miR-92b-3p), and 1062 mRNAs were obtained ( Figure 2A ).

Regarding the top 10 downregulated DEcircRNAs, the four databases were utilized to predict 171 overlapping miRNAs (four circRNAs without predicted miRNAs were removed). 11 miRNAs were selected by intersecting the predicted miRNAs and DEmiRNAs. 1753 target mRNAs for these miRNAs were identified via MiRTarBase. 6 circRNAs (including hsa_circ_0003923, hsa_circ_0027702, hsa_circ_0000657, hsa_circ_0008833, hsa_circ_0000042, hsa_circ_0069613), 11 miRNAs (including hsa-miR-493-5p, hsa-miR-132-3p, hsa-miR-21-5p, hsa-miR-24-3p, hsa-miR-22-3p, hsa-miR-556-3p, hsa-miR-655-3p, hsa-miR-31-5p, hsa-miR-369-3p, hsa-miR-889-3p, hsa-miR-495-3p), and 1753 mRNAs were obtained ( Figure 2B ).

The potential roles of the target mRNAs were investigated via GO and KEGG analyses. KEGG analysis revealed that target mRNAs for the upregulated DEcircRNAs in the ceRNA network were enriched in oocyte meiosis, cell cycle pathways, and protein processing in the endoplasmic reticulum. GO analysis demonstrated that these target mRNAs were enriched in nuclear division, nuclear chromosome segregation, and chromosome segregation in NPC patients with PD ( Supplementary Figure 4 ).

Regarding the downregulated DEcircRNAs in the ceRNA network, the KEGG analysis demonstrated that their target mRNAs were enriched in Th1 and Th2 cell differentiation, natural killer (NK) cell-mediated cytotoxicity, and axon guidance in NPC patients with PD. GO analysis indicated that these target mRNAs were enriched in antigen receptor-mediated signaling pathways and the signaling pathways of immune response-regulating and response-activating cell surface receptors ( Supplementary Figure 4 ).

Blood samples at T1 and T4 from 20 Class 1 patients and 20 Class 4 patients were randomly selected (The T1 and T4 samples from the six patients previously sent for sequencing were both included). The random three upregulated circRNAs (hsa_circ_0004277, hsa_circ_0003684, and hsa_circ_0000721) identified in the ceRNA network were validated via RT-qPCR. The results of the RT-qPCR aligned with the WTS data. No significant differences in the expression levels of hsa_circ_0003684 and hsa_circ_0000721 between Class 1 and Class 4 patients at T1 and T4 were found (P > 0.05) ( Figures 3B, C, E, F ). However, as shown in Figures 3A, D , significant differences in the expression levels of hsa_circ_0004277 between Class 1 patients and Class 4 patients at T1 and T4 were observed (P < 0.01). The results all showed that the expression level of hsa_circ_0004277 in the ascending group was higher than that in the descending group. Among them, the expression level of the T1 rising group was 1.305 times that of the falling group, and the expression level of the T4 rising group was 1.605 times that of the falling group.

Hub genes in the ceRNA network were selected via PPI networks ( Supplementary Figure 5 ). Hub genes including RPS27A, H3C13, TP53, UBA52, H4C6, BRCA1, CDK1, EZH2, PLK1, and CCNA2 were found to play vital roles in the ceRNA network based on the upregulated circRNAs. Additionally, hub genes including TP53, MYC, CTNNB1, AKT1, EGFR, EP300, RPS27A, HDAC1, ACTB, and SRC were identified in the ceRNA network based on the downregulated circRNAs.

Through a synchronous prediction, 35 RBPs were selected from the upregulated DEcircRNAs via the circAtlas and RBPsuite databases. 19 RBPs were identified from the downregulated circRNAs. The circRNA–RBP networks are visualized in Supplementary Figure 6 .

Most included NPC patients received 30 to 33 sessions of radiotherapy over a period of 45 to 60 days. Diagnosis was typically confirmed within three to five days of admission, with radiotherapy starting around day seven to ten and concluding by day 60 to 90. During this period, our research team was able to closely monitor changes in PD among NPC patients undergoing radiotherapy. This study identified four PD trajectories of NPC patients undergoing radiotherapy via LCGM: the decline distress group (11.0%), Class 1; the low-stable distress group (58.0%), Class 2; the no distress group (10.7%), Class 3; and the rise distress group (20.3%), Class 4. Patients with PD were mainly from Classes 1, 2, and 4, indicating the lasting PD during NPC radiotherapy. This finding is aligned with other longitudinal studies in cancer patients (15). Notably, Class 1 patients experienced high PD levels before receiving radiotherapy, which then declined steadily through the course of treatment. We hypothesized that these patients may have utilized more active coping strategies and received greater social support to help them cope with the physical and emotional changes caused by cancer. Our results indicated that compared to Class 4 and Class2 patients, those in Class 1 had higher scores for coping strategies and PSS, suggesting that adaptive coping may prevent the deterioration of psychological health (29).

Homeostasis and normal physiological activities in the body are maintained by accurately regulating gene expression. In addition to investigating the sociodemographic data of PD in NPC patients undergoing radiotherapy, this study explored the influence of genetic factors. The ceRNA and circRNA-RBP networks were built to provide insights into the pathophysiology of PD in NPC patients undergoing radiotherapy. KEGG and GO analyses also provided perspectives on the mechanism and regulation of PD in these patients. However, there are several limitations in this multicenter study. First, this study involves a limited sample size, which can lead to selection bias and affect the accuracy of our conclusions. Moreover, the follow-up period is limited to 90 days at only four time points. Future studies involving larger cohorts and more extensive and frequent follow-ups are required to validate our findings. Additionally, some tools used to measure PD are self-reported, which may cause response bias. Finally, the number of blood samples used for genetic analysis is limited, and only upregulated circRNAs were validated via PCR, thereby leading to insufficiently robust evidence. The study does also not address how radiotherapy dosage or concurrent chemotherapy regimens might interact with psychological distress trajectories. Therefore, future research requires more rational design with larger samples.