Category: 3. Business

Sampling procedure

In the context of the Horizon Europe Project “InBestSoil”, the data collection focused on arable management practices in Switzerland. Specifically, those practices related to soil health and soil conservation undertaken within the 2022/2023 production season. Farm selection for the survey was based on specific criteria to ensure that the data collection accurately represented arable agricultural practices in Switzerland. These criteria were designed to target farms that were significantly involved in arable agriculture, which is crucial for assessing arable soil health management practices. Eligible farms were required to meet the following criteria:

• Grow wheat in the preceding season (2021/2022).

• Farm at least 3 hectares of arable land in the preceding season (2021/2022).

• Arable land must have comprised at least 20% of the total farmed area in the preceding season (2021/2022).

We entered a data sharing agreement with the Federal Office of Agriculture to enable our survey campaign via access to contact information of all farmers who met the above selection criterion (see the supplementary material in the data repository for a copy of this contract)¹. The Federal Office of Agriculture implemented our selection criterion on the agricultural data that they collect on a yearly basis from the direct payment applications of all Swiss farmers. Note, at the time of our application to the Federal Office of Agriculture, data for the production season 2022/2023 was not available. This is why we use data from the preceding production season for specifying the selection criteria, as this was the latest data available at the time, from which the Federal Office of Agriculture could make an assessment of which farm contact details to share with us for the survey.

In August 2023, we received the contact details of 15,023 farmers who qualified for the survey from the Federal Office of Agriculture’s records. The information we received included the email address, farm identification number, language spoken, name and form of address. However, as per our data sharing agreement with the Federal Office of Agriculture, this data was allowed exclusively for our use in this project and cannot be shared with any outside partner not party to the aforementioned data sharing contract. The contact data of farmers that was received from the Federal Office of Agriculture will be kept for the duration of the InBestSoil project and stored securely on private institutional servers in encrypted files. All contact information will be deleted at the conclusion of the project (December 2026) and all data presented herewith is strictly anonymised to protect the data and identities of the farmers who took part in the survey. Moreover, we have taken measures to prevent any farmers from being identified via their answers (for example variables such as manager age, wheat areas grown, location etc. have been classified into more homogenous categorical groups), which means that the data we present here is slightly different to the data that we have available for our own analyses, as agreed under the data sharing agreement with the Federal Office of Agriculture.

Survey design and content

While adoption of agricultural practices certainly varies with farm characteristics such as size, labour availability, or participation in agri-environmental schemes, these factors alone are not sufficient to explain farmer behaviour. There is no single set of drivers that consistently predicts adoption across studies or regions⁴³. Instead, adoption depends strongly on local contexts, and the interplay of economic, social, and psychological factors⁴⁴. To capture the complexity of adoption behaviour, the survey included questions on farmers’ priorities, perceptions, self-assessed competencies, and personal goals, as well as their exposure to peer practices, participation in training and advisory services, and sources of information. These dimensions are important because farmers do not make decisions in isolation; their attitudes towards risk, innovation and environmental values can influence their decisions alongside financial considerations. Such data contribute to a more thorough understanding of the multifaceted factors influencing soil health-related decisions. The inclusion of these variables also offer valuable insights into the barriers and drivers of sustainable soil management, essential for shaping targeted and effective agricultural policies and support programs.

The full survey is available within the data repository in French, German and English¹. The final survey was developed over the course of a year, including revisions resulting from three rounds of consultation with external stakeholders, internal consultation and testing with farmers. All participants in the survey were asked to give their informed consent by ticking a box in the online questionnaire, confirming their agreement to participate in the study. Additionally, participants consented to the linking of secondary geographical data with their responses, which was also confirmed by ticking a separate checkbox in the survey. Once the participants had agreed to these, the survey was administered uniformly following the structure outlined below. All questions appeared in the same order and, only if certain exclusion criteria were met – such as when their previous answer ruled out any further sub-questions – were some sub-questions hidden from the view of participants. Inclusive of all sub-questions, the survey contained 57 questions, and answering the questionnaire took farmers a median time of 23 minutes.

The survey design was based on previously implemented surveys regarding agricultural production practices in Switzerland^{45,46,47,48,49}. Specifically, questions on farm information and participation in soil-related programmes were included to assess farmers’ engagement with policy incentives and voluntary schemes. The inclusion of personal characteristics aimed to understand demographic drivers of management behaviour. The questions on management practices were developed in close collaboration with experts from the soil science and agricultural extension fields, and were cross-checked with relevant literature. Data on milling wheat production and related input use were collected to link agronomic decisions with productivity outcomes. Information on structural farm characteristics, such as farm type, location, and land tenure, provides context for understanding the decision-making environment and potential constraints faced by farmers. Finally, a strong focus was placed on behavioural and attitudinal factors, including information sourcing, perceived risks, and personal goals, to account for the cognitive and motivational dimensions of farmer behaviour. The following section provides an overview of the variables investigated within each of these question groups. The collected data are documented in the accompanying datasets¹. Each question group corresponds to a clearly defined set of columns.

Demographic details (Primary dataset columns B-H)

Age, duration farm responsibility, gender, full time equivalent and whether the farm succession is already secured.

Participation in soil health programmes (Primary dataset columns H-S)

Organic farming support, soil cover scheme, reduced tillage scheme, herbicide-free farming scheme, pesticide-free farming scheme, efficient fertiliser use, wider row planting, beneficial insect strip, precision application, cantonal soil health support, cantonal input reduction support, cantonal investment and equipment support.

Management Practices (Primary dataset columns T-CA)

An overview of all management practices addressed in the survey, including their descriptions and the typical machinery used, is provided in Table 1. Farmers were asked about their knowledge about the practices, the application as well as the frequency of application within the last 10 years and whether they know other farmers that use the practice. The practices covered by our survey were selected based on the input of soil scientists and agricultural extension workers based in Switzerland.

Milling Wheat Production (Wheat dataset columns B-M)

Production standard, hectares of milling wheat grown, yield milling wheat, yield milling wheat over last five seasons, quantity synthetic fertiliser, quantity organic fertiliser, sowing density, number of biostimulant treatments, number of herbicide treatments, number of fungicide treatments, number of insecticide treatments and number of plant growth regulator treatments.

Structural Farm Characteristics (Primary dataset columns CD-CP)

Family members employed, farm focus (arable, livestock, permanent crop, others), full time or part-time farm, percentage of rented land, whether the soil has been assessed and a soil management plan exists.

Training and Advice (Primary dataset columns CQ-CZ)

Advice agricultural adviser, advice agricultural retailer, advice cantonal or national institution, consult other farmers, consult social media channels, consult publications or webpages, participation equipment demonstration, participation farmer discussion or training group, participation farm demonstration, participation course.

Behavioural and Attitudinal Factors (Primary dataset columns DA-EK)

Respondents’ self-assessment of their perceived influence of the weather on crop production and ambitiousness of self-set production goals.

Respondents’ self-assessment of their willingness to take risks in the domains of; agricultural production, investment in agricultural technology and crop protection.

Respondents’ self-assessment of their confidence in being able to; find solutions to arable production challenges and achieve production goals by harvest end.

The respondents self-reported importance of the following aspects in decision making;

Maximising yields, minimising input costs, minimising time or labour requirements, minimising production risks, minimising farm exposure to weeds or pests or diseases, adapting to weather patterns, adapting to farmland conditions, improving soil health or structure or fertility, improving biodiversity, minimising environmental impact, expanding farm land, adapting to crop market developments, adapting to changes in direct payment rates or regulations, seeking professional agronomic, seeking casual advice from friends or colleagues and seeking peer approval.

Ethical approval and pre-registration

The survey campaign and research design were both approved separately by the ETH Zürich Ethics Commission as proposal 2023-N-212 as well as the FiBL Ethics Committee as proposal FSS-2023-006. Copies of the approval letters are included in the supplementary material¹. Before launching our survey, we also submitted two research plans for pre-registration of hypotheses via the online platform AsPredicted operated by the University of Pennsylvania (link: AsPredicted). For further information on these, see AsPredicted #153145 and AsPredicted #153146 that were registered on 29^th November 2023.

Survey implementation

The survey was implemented as an online survey formulated with Lime Survey and distributed via email. All eligible farms received an individualised email addressed personally to the recipient and a survey link, connected with a unique token to enable us to link the farmer responses with secondary data available for each farm. The participants were asked to give their permission for this by approving the terms and conditions we made available to them regarding how their data would be handled. By agreeing to the disclosure agreement, the farmers gave their permission for the anonymised data, that they subsequently provide through the survey, to be used exclusively for science and research purposes. Farmers were also given the option to opt out of the survey at any time, with no explanation needed. To incentivise participation in the survey we offered the opportunity to enrol in a lottery of 100 supermarket vouchers worth CHF 150 each and the option to receive a personalised results report comparing the farmers’ answers to the answers of other similar farms. The individualised reports were administered via a bespoke app created using R-Shiny (see technical validation section below for further details).

Prior to the full survey launch, a pilot survey was conducted on a random sample of 1% of eligible farms (150 farms) to test the survey’s functionality and to refine any issues. The pilot survey launched on 30^th November 2023, and the full survey went live six days later, on 6^th December 2023. The survey was closed on 31^st January 2024, after a response period of nearly two months.

Data cleaning

To minimise errors already at the point of data entry, the survey was designed to allow only predefined values or plausible numeric ranges for most variables. Wherever this was not technically feasible, such as in open-text fields or free numeric input, we conducted systematic data cleaning after data collection. Data cleaning involved addressing inconsistencies and missing values. In cases where values were deemed implausible or outliers, they were either removed or corrected if sufficient data from other columns was available. This cleaning procedure was applied to variables related to plant protection product treatments, yield, sowing density, labour input, and demographic information. We include the following to illustrate the approach we took as an example (note all processing codes are available in the supplementary material which outline these decisions on a line-by-line basis):

If in the labour units column, an entry was listed as 48, which was inconsistent with the farm area, this value was corrected to 4.8 using a related column for recalculation. Similarly, we proceeded for the variable age: if a data entry was obviously wrong, such as a year of birth recorded as 60 instead of the demanded format YYYY (1960), and the farmer had entered the column of farming experience 40 years, the value was corrected to ‘1960’ based on logical inference. If no reliable correction could be made, the value was marked as ‘NA’ (Not Available).

To ensure anonymity, apart from removing precise geographical information we also grouped continuous variables such as age and farming experience into categories (e.g. age_group and years_experience_group). The data was anonymised, and no specific details were included that could link individual responses to specific farms. No randomisation was applied to the data. With regard to the secondary data, we also took measures to prevent identification by rounding the variables to the nearest integer (the codes for the processing of this data are also available in the supplementary material).

Development of an Artificial Intelligence-assisted CRISPR-Cas Scan (AIL-Scan) strategy based on an ESM large language model

We assumed that by embedding the functional feature with protein primary sequences, we could trace the natural evolution rules and identify the CRISPR-Cas proteins in the metagenomics data directly without sequence alignments. To identify the CRISPR-Cas proteins, we developed an Artificial Intelligence-assisted CRISPR-Cas Scan (AIL-Scan) strategy (Fig. 1a). It includes the following steps:

1. 1.

   CRISPR-Cas training data is created by extracting CRISPR-associated (Cas) proteins from the NCBI database, classifying them by genes, and removing redundant sequences.

2. 2.

   Supervised fine-tuning of ESM on the CRISPR-Cas training data based on the biological information to predict the Cas protein.

3. 3.

   Feature analyses of Cas proteins, including cleavage activity, CRISPR-loci type, CRISPR loci-length, direct repeats, spacers, evolutionary analyses, MSA, and structures.

We generated our training data using reviewed NCBI gene data. We annotated the Cas1, Cas2, Cas3, Cas4, Cas5, Cas6, Cas7, Cas8, Cas9, Cas10, Cas12, and Cas13. Non-Cas proteins were extracted according to the following rules, without the annotation of Cas, and removing the proteins with sequence similarity over 40%. The Cas protein database was separated into a training or validation database using CD-HIT-2D with a 40% identity threshold to remove the redundant sequences and avoid overfitting. We collected 76567 non-redundant positive sequences and 13047 non-Cas proteins, which were deposited in NCBI before July 5, 2023 (Supplementary Fig. 1). The maximal protein length is less than 1764 amino acids. To obtain the best classification, we introduced the “focal loss” in the classification to solve the unbalance of the input data. We obtained the best model during the 13^th Epoch of model training and obtained 97.75% accuracy for the ESM 2 model with 650 million (650 M) parameters (Supplementary Fig. 2). Using the 15 billion (15B) parameters model, we achieved the best performance in the 9^th Epoch with 98.22% accuracy (Supplementary Fig. 2). This model maintained consistent performance, achieving an accuracy 97.68% on the independent dataset, i.e. TestSet2024, which contains sequences deposited in NCBI from July 6, 2023, to Oct 28, 2024 (Supplementary Tables 1–3). These results indicate a robust generalization of this model. The accuracy and prediction speed of AIL-Scan is comparable to the CRISPRcasIdentifier, which integrates HMMs and machine learning (Table 1 and Supplementary Fig. 3). CASPredict performed with the highest speed among the four software, although its accuracy is lower than the machine learning based software, i.e., AIL-Scan and CRISPRcasIdentifier. However, the NCBI data has been partially annotated by the HMM model, so we turned to validate AIL-Scan’s capability in recognizing “unseen proteins”. We utilized a recent dataset of 3601 Cas12 family protein sequences²⁰, in which 3521 sequences (97.8%) had less than 90% similarity with the training set, meanwhile 3351 sequences (93.1%) had less than 40% similarity with the training set. This test set is named TestSet2025 and is significantly distinct from the training set in sequence space, making it suitable for evaluating generalization ability. AIL-Scan successfully identified 3182 Cas12 proteins, in contrast, the HMM model identified 1240 sequences, demonstrating the strong generalization capabilities of AIL-Scan. Considering the resource consumption, the 650M model is sufficient for the Cas prediction. We used ESM embeddings to reduce dimensionality with t-SNE for 77684 sequences and discovered that ESM can distinguish the differences in various Cas classifications. The ROC curves and AUC indicate the probability that the positive sample’s decision value is greater than the negative sample’s decision value for all the Cas and non-Cas proteins (Fig. 1b). The test loss and test accuracy also indicate that the model generalizes correctly and performs well on unseen data (Fig. 1c). We evaluated the model robustness using the 5-fold cross-validation. The average accuracy is 0.9786 and the standard deviation is 0.0013 (Supplementary Table 4).

We use the Global Microbial Gene Catalog (GMGC) metagenomic database for the Cas protein discovery²¹. We selected 50,000 bins with high quality from GMGC and extracted 20,000 MAGs, including CRISPR-loci, to test the performance of AIL-Scan. The protein sequences were predicted by Prodigal software²². We collected ca. 20,000,000 protein sequences shorter than 1500 amino acids for prediction. In comparison with the established methods, the AIL-Scan predicts 1379 Cas12a sequences.

Development of a trans-cleavage activity prediction model

The trans-cleavage activity of Cas12a has been used in various applications. Although many CRISPR-Cas12a proteins have been identified, few of them have been tested in the trans-cleavage experiments. Therefore, the main challenge encountered during this study lies in dealing with a small sample size coupled with high-dimensional embeddings, which often leads to convergence issues when employing most models. A total of 69 labeled Cas12a proteins (including three known Cas12a) were included in our analysis (Supplementary Data 1). Their trans-cleavage activities were assessed by the fluorophore-quencher (FQ) reporter assay. The trans-cleavage activity was defined as proteins displaying fluorescence intensity twice that of the negative control. Thirty-three proteins were classified as active in trans-cleavage activity, and the remaining 36 proteins were categorized as inactive. To evaluate the performance of our predictive model, a test set comprising 13 randomly selected proteins (approximately 20% of the sample) was used, while the remaining 56 proteins were employed for training purposes. Initially, we recorded the last embedding layers based on our fine-tuned ESM model for all labeled Cas12a protein sequences. These embeddings (1280 dimensions) were utilized as covariates to predict trans-cleavage activity.

Different forms of decision tree models are evaluated in this task. The results of our study demonstrate that Light Gradient Boosting Machine (LightGBM) achieves the highest accuracy among mainstream machine learning models, with an accuracy rate of 69.2% on the test set trained on embeddings. To address dimensionality-related challenges, principal component analysis (PCA) was employed to extract essential embeddings, with prediction performance evaluated across 2–15 principal components. Alongside PCA, we compared 31 alternative methods, including t-SNE, UMAP, and raw data. Detailed comparisons, training procedures, and results are provided in Table 2, Supplementary Table 5, and the supplementary notes. LightGBM, CatBoost, and RandomForest achieve the accuracy of 92.3% in the test set (12 out of 13 proteins are correctly labeled) with 4, 6, and 8 principal components, respectively. We can see that compared to training models directly with embeddings, extracting essential dimensions with PCA provides higher accuracies in predicting trans-cleavage activity (Supplementary Table 5). However, this model is still limited by the small dataset, more experimental data would improve its prediction accuracy. Additionally, we tested our prediction model on two unreported Cas12a proteins, i.e., the trans-cleavage activity of two Cas12a candidates: ArCas12a_1 (derived from Agathobacter rectale) and LeCas12a_3 (derived from Lachnospira eligens_B). Our model predicted that ArCas12a_1 has trans-cleavage activity but not LeCas12a_3. In the experiment, ArCas12a_1 demonstrated significantly stronger trans-cleavage activity than the negative control, while LeCas12a_3 did not (Supplementary Fig. 4). These experimental outcomes were consistent with our model’s predictions, supporting the generalizability and robustness of the prediction model.

CRISPR-Cas12a loci predicted from the metagenomics

We did further feature analyses of Cas12a candidate proteins. Phylogenetic analysis of Cas12 proteins suggests that the identified Cas12a proteins fall into the Cas12a clade (Fig. 2a). The classical CRISPR-loci, comprising essential elements such as Cas1, Cas2, and Cas4, play a pivotal role in type classification. To delve into these features, we employed AIL-Scan to predict Cas1, Cas2, and Cas4 proteins within the same CRISPR loci adjacent to the Cas12a sequence. Subsequently, we meticulously verified 300 predicted CRISPR loci to gain deeper insights manually. Normally, Cas12a is considered to have a unique CRISPR locus, comprising Cas1, Cas2, and Cas4. Intriguingly, the observed count of Cas1, Cas2, and Cas4 proteins was notably lower than that of Cas12a, suggesting the absence of these small Cas proteins in some Cas12a loci (Fig. 2b, c). Further stratification based on the number of integrase proteins led to the classification of CRISPR loci into eight distinct subtypes. The distribution of integrase proteins across these subtypes exhibited a sparse pattern (Fig. 2d). Notably, subtype VIII lacked any integrase proteins, subtype I encompassed Cas1, Cas2, and Cas4, while subtype VI exclusively featured Cas2. This nuanced classification sheds light on the diversity within CRISPR loci and underscores the intricate variations in the composition of integrase proteins among different subtypes. Our observations may provide unreported perspectives on correlations among different CRISPR-Cas systems and integrase proteins. Remarkably, the analyses using the 1000 predicted CRISPR Cas12a loci without manual verification show a strikingly similar distribution pattern as the result from the 300 manually confirmed ones, indicating this distribution is a universal phenomenon (Supplementary Fig. 5). To provide further insights, we measured the length of CRISPR loci, beginning from the start of the Cas12 protein and concluding at the first spacer. Subtype VIII emerged as the shortest, spanning mere 4200 bp, while subtype I is the longest, extending over 6100 bp. Particularly noteworthy were certain subtype I CRISPR loci exhibiting extraordinary lengths of up to 6700 bp, raising the possibility of harboring enigmatic protein elements (Fig. 2e). Aligned with the integrase variation, the numbers of spacers notably decreased in subtypes IV, VI, and VIII, underscoring the pivotal roles of integrases in spacer capture (Fig. 2f). Despite the divergence in spacer numbers, the stem-loop region corresponding to direct repeat sequences remained conserved (Fig. 2g). This consistent conservation hints at a shared structural element, emphasizing the importance of the stem-loop region in CRISPR loci across different subtypes.

To explore the distribution of the discovered proteins in the organisms, we constructed a phylogenetic tree using 300 candidate Cas12a proteins, which were manually verified, along with three known Cas12a (LbCas12a, FnCas12a, and AsCas12a). 232 Cas12a proteins from the Lachnospiraceae family cluster into one clade. Within this clade, subclade 1 consisted of 62 subtype I Cas12a proteins, 81 subtype VII Cas12a proteins, and a modest representation of other subtypes. Notably, subtype I and subtype IV emerge as the principal constituents within Subclade 2. Furthermore, Subclade 3 is marked by the exclusive presence of 28 subtype VIII Cas12a proteins originating from the Acutalibacteraceae family. It is worth noting, 94.6% of the identified Cas12a proteins originate from enteric microorganisms (Fig. 2h), which may be due to the ease of recovering high-quality genomes from enteric microorganisms. Additionally, the thermostable YmeCas12a (subtype I) is adjacent to subtype I Cas12a proteins (Supplementary Fig. 6).

Cas integrases in CRISPR loci

New insights highlight the structural diversity and functional roles of Cas integrases in CRISPR loci^{23,24,25,26,27}. Cas1, Cas2, and Cas4 are essential for integrating foreign DNA into bacterial CRISPR systems, which generates bacterial immunity²⁶. AlphaFold2²⁸ was applied to predict all protein structures in the eight distinct subtypes, providing insights into their variation, respectively (Fig. 3 and Supplementary Fig. 7). Cas1 proteins, encompassing 92–331 amino acids, are classified into eight types based on structure and sequence (Fig. 3a, b and Supplementary Fig. 7b). Type 8 is the most prevalent Cas1 protein, resembling AfCas1 (PDB: 4N06)²⁹ and its N-terminal and C-terminal domains (NTD, CTD) contain with key catalytic sites in specific helices and loops (Supplementary Fig. 7c). Structural differences across types were analyzed via the Dali server³⁰. The variation in CTD elements does not necessarily hinder foreign DNA acquisition³¹, emphasizing their structural flexibility. Cas2 proteins, containing 70–146 amino acids, also fall into eight subtypes, with type 8 showing notable structural similarities to E. coli Cas2 (PDB: 5DQT)³² but with unique N-terminal helices (Fig. 3c, d and Supplementary Fig. 7d–f). Other subtypes exhibit varied structural deficiencies, such as missing β-sheets or helices, affecting dimer interfaces and potentially altering DNA binding. This diversity underlines Cas2’s adaptability within Cas1–Cas2 complexes (Supplementary Fig. 7f)³³. Cas4 proteins, comprising 79–206 amino acids, exhibit eight types (Fig. 3e, f and Supplementary Fig. 7g, h), with type 8 resembling I-C Cas4 (PDB: 8D3Q)²⁴ but lacking specific helices critical for protospacer cleavage. Structural differences across subtypes, such as missing helices or β-sheets, impact spacer insertion and integration within CRISPR systems (Supplementary Fig. 7i). These findings broaden our understanding of Cas4 structural variations and their functional implications in bacterial immunity. The detailed structural features of integrases are analyzed in the Supplementary Note.

Cas12a proteins in the subtypes

The differences in the Cas12a structures are key features of the Cas12a subtypes. We analyzed the motifs of the Cas12a sequences and discovered conserved and distinct motifs in the different subtypes, which are key for the Cas12a functions (Supplementary Fig. 8). The analysis revealed that the catalytic residues within the RuvC and Nuc domains are highly conserved among all subtypes, reflecting their critical roles in enzymatic function. Specifically, the first catalytic aspartate in the triad resides within the conserved motif IGIFRGEERN. The second catalytic glutamate displays subtype-specific distributions, appearing as MED in subtypes I, IV, V, and VI, as M/LEN/D in subtype II, and as MEK/D in subtype VIII. The third catalytic aspartate is consistently located in the motif DADANG, specifically at the second “D”. Additionally, a highly conserved TSKIDP motif was identified across all subtypes, indicating a shared functional mechanism. Other conserved motifs showed variability among subtypes, suggesting distinct sequence characteristics while maintaining overall catalytic and structural integrity. We also built the structure models of 300 Cas12a proteins using AlphaFold2, except for the failed construction, and calculated the root mean square fluctuation (RMSF) for all candidate Cas12a proteins within one subtype (Supplementary Fig. 9). The detailed analyses are appended in the Supplementary Notes. The RMSF reflects the residue-wise structural difference within one subtype. The results suggested that, despite an overall conserved structural architecture, specific regions within the proteins exhibit variability that may reflect structural adaptations specific to each subtype.

Cas12a proteins have distinct cis– and trans-cleavage activities

Cas12a processes the pre-crRNA transcripts into mature crRNA by its endoribonuclease activity. Then the Cas12a–crRNA complex efficiently cis-cleaves a double-stranded DNA (dsDNA), which is initiated by a PAM motif recognition. The cleaved DNA segment that remains bound then induces non-specific degradation of single-strand DNA (ssDNA) (Fig. 4a).

Therefore, we evaluated the RNA binding efficiency, DNA binding efficiency, cis– and trans-acting DNase activities of sixteen Cas12a proteins from eight subtypes derive from Anaeroglobus micronuciformis (AmCas12a), Eubacterium_G ventriosum (EvCas12a_1 and EvCas12a_2), Erysipelatoclostridium sp. (EspCas12a), Ruminococcus_E sp. (RspCas12a_1 and RspCas12a_2), Agathobacter rectale (ArCas12a), Lachnospira eligens (LeCas12a_1 and LeCas12a_2), UBA3388 sp. (UBACas12a), RC9 sp. (RCCas12a), CAG-127 sp. (CAGCas12a), Ruminococcus_E bromii_B (RbrCas12a_1, RbrCas12a_2, RbrCas12a_3 and RbrCas12a_4) (Fig. 4, Supplementary Fig. 10 and Supplementary Table 6). Remarkably, the direct repeat sequence of these candidate Cas12a proteins is conserved alongside their celebrated counterparts, i.e., LbCas12a (Fig. 2g and Supplementary Fig. 11). Therefore, we chose LbCas12a as the positive control in the following assays, as well as its crRNA scaffold in the screening step. All the Cas12a proteins show RNA and DNA binding ability as expected (Fig. 4b, c, Supplementary Fig. 10c, d, and Supplementary Table 7). However, the DNA binding ability of subtype I and subtype VIII are higher than other Cas12a proteins. According to the inherent trans-DNase activity of Cas12a, as well as the 4 bp PAM length, we developed a simple and efficient PAM detection method. We constructed 6 short dsDNA target arrays by annealing 256 kinds of PAM sequence primer pairs in each well, which target EMX1 site1, DNMT1 site1, FANCF site1, MERS site1, eGFP site1, and eGFP site 3 (Supplementary Table 8). Each dsDNA target was incubated with candidate Cas12a proteins, crRNA and FAM-BHQ reporter to detect fluorescence of each reaction system (Fig. 4d). Using this assay, we determined the PAM preference of EvCas12a_2, AmCas12a, RspCas12a_2, CAGCas12a and RbrCAS12a_1, EcCas12_2, RspCas12a_2, and CAGCas12a recognize T rich PAM, but AmCas12a prefer G-start PAM, RbrCas12a_1 recognize 5-GTV-3 PAM (Fig. 4e–i and Supplementary Figs. 11, 12).

To corroborate the cis-acting DNase activity of candidate Cas12a proteins, we incubated Cas12a proteins with a crRNA and a linearized plasmid dsDNA. All linearized dsDNA were degraded by candidate Cas12a proteins with comparable efficiency to LbCas12a at 37 °C, with the exception of RCCas12a (Fig. 5a and Supplementary Fig 13a). Sanger sequencing of the cleaved DNA ends revealed that AmCas12a introduced INDELs at 18 in NTS and 23 in TS, consistent with other Cas12a orthologs (Supplementary Fig. 13e, f). However, most Cas12a variants exhibited diminished DNase activity, resulting in the production of uncleaved DNA at room temperature (RT), except for subtype VIII Cas12a proteins, which lack integrases. (Fig. 5b and Supplementary Fig. 13b). Subtype II Cas12a variants are slightly less active than LbCas12a in single-strand (ssDNA) degradation, while EspCas12a, EvCas12a_1, EvCas12a_2, and ArCas12a exhibited moderate activity. In contrast, the other Cas12a variants displayed notably lower activity (Fig. 5c and Supplementary Fig. 13c). Most of these Cas12a proteins represent considerable cis cleavage activity but are a bit different in trans-cleavage activity compared to LbCas12a. The ion preference assay reveals that these Cas12a proteins can be activated by Mn²⁺, similar to the LbCas12a³⁴. Divalent Mg ions prove ineffective in activating the trans ssDNA cleavage activity of low-activity Cas12a variants, and Mn²⁺ cation emerges as the catalyst for their trans DNase activity. (Fig. 5d and Supplementary Figs. 13d and 14) To investigate the genome-editing ability of candidate Cas12a in eukaryotic cells, we selected 6 target sites with canonical PAM, which can be recognized by all the tested Cas12a (Fig. 5e and Supplementary Table 9). AmCas12a exhibits an average editing efficiency of 49.6% across six sites, with remarkable peaks at sites 3 (85.4%) and 6 (84.9%). In contrast, EvCas12a_2 displays an average editing efficiency of 20.3%, with its highest performance observed at site 1 (25.8%). RspCas12a_2 and RbrCas12a_2, which lack integrase in the loci, yield modest average editing efficiencies of 14.3% and 17.8%, respectively, with notable peaks at site 3 (26.3% and 37.3%, respectively). ArCas12a shows comparable average editing efficiencies with AmCas12a (45.4%), which gets notable peaks at site 3 (75.8%). LeCas12a_1 shows an average editing efficiency of 6.2% and a maximum efficiency of 25.7% at site 2. UBACas12a exhibits nearly negligible editing efficiency, with the highest activity reaching 2.1%. At site 4, CAGCas12a and LeCas12a_2 demonstrate peak genome-editing efficacy, at 81.7% and 73.8%, respectively, with mean editing efficiencies of 28.8% and 26%. AsCpf1 attains an impressive average editing efficiency of 65.5%, with its maximum at site 6 (84.7%). Finally, LbCas12a shows an average editing efficiency of 25.6% and a maximum efficacy of 53.5% at site 6.

The AmCas12a–crRNA binary complex

The protein sequence identity of 16 candidate Cas12a proteins to AsCas12a, FnCas12a, and LbCas12a are low, ranging from 30%-46% (Fig.6a and Supplementary Fig. 15). In the three-dimensional structural landscape, Cas12a proteins within the same subclade exhibit a high degree of structural similarity. However, AmCas12a presents a subtle deviation, distinguishing itself somewhat from its subclade I Cas12a counterparts (Fig. 6d, f and Supplementary Fig. 15).

To understand the molecular details underlying the RNA binding behavior of AmCas12a, we achieved the cryo-EM map of the crRNA binding complex, which consists of AmCas12a and a 44-nt crRNA, at 2.9 Å resolution (Fig. 6b, c, Supplementary Figs. 16 and 17, and Supplementary Table 10). The AmCas12a–crRNA structure maintains a bilobed architecture (Fig. 6c), similar to other Cas12a structures^35,36. Nonetheless, it is noteworthy that the AmCas12a–crRNA complex exhibits a distinct conformation when juxtaposed with its counterparts. Specifically, an observable rotational variance is discernible within the REC domain of AmCas12a when compared to the LbCas12a–crRNA and FnCas12a–crRNA complexes. Relative to LbCas12a and FnCas12a, the REC1 domain of AmCas12a presents a deviation of 7.3° and 9.4°, respectively. Simultaneously, the REC2 domain of AmCas12a manifests a rotational disparity of 4.8° and 6.2°, respectively (Supplementary Fig. 17d, e).

As observed in the LbCas12a and FnCas12a crRNA binary structures, the repeat-derived pseudoknot in the 5’ handle of the crRNA is ordered. However, the crRNA conformation is markedly different from that of the crRNA bound by LbCas12a or FnCas12a. Due to the flexibility of the spacer-derived part of crRNA, it’s almost unclear in the Cas12a–crRNA binary complex^35,36. Notably, an extra RNA stem formed by A(1)–A(5) and U(18)–U(22) within the crRNA spacer region makes a part of spacer region including seed sequence well-defined in the central cavity of AmCas12a and adopt an A-form-like helical conformation, but A(−10)–G(−6) and G(6)–A(15) nucleotides of crRNA are unclear (Fig. 6b and Supplementary Fig. 18). To accommodate the double RNA stem substrate, the REC lobe of AmCas12a rotates away from the NUC lobe. Unsurprisingly, the docking of crRNA to Alphafold-generated AmCas12a causes a severe clash in the REC domain (Supplementary Fig. 15c). The attainment of conformational integrity within the extra RNA stem is orchestrated by intricate interplays involving the ribose and phosphate moieties of the crRNA backbone, engaging in multiple interactions with specific residues within the WED, REC1, and RuvC domains of AmCas12a (Fig. 6e). These include residues T19, H751, K522, and H861 from the WED domain, Y50 and R168 from the REC1 domain and Q1003 from the RuvC domain, all of which are conserved with Cas12a orthologs, except Q1003 which form a hydrogen bond with the phosphate of U(18) (Supplementary Fig. 18). Distinct from the FnCas12–crRNA complex, the spacer segment of crRNA major interacts with the WED domain of AmCas12a.

Compared to the LbCas12a–crRNA complex and FnCas12a–crRNA complex, the divalent Mg ions are in the same location (Supplementary Fig. 17a–c). Consistent with a seed sequence-dependent mechanism of DNA targeting and in broad agreement with previous analyses of AsCas12a, LbCas12a activities in vivo, and FnCas12a activities in vitro^35,37,38, cleavage of DNA substrates with single-nucleotide mismatches in the seed segment was almost completely impaired, while mismatches in the PAM-distal region of the DNA target were mostly tolerated (Fig. 6g).

Specific detection of single-nucleotide mutation by AmCas12a

Cas12a is a promising tool in the next-generation molecule diagnosis, however, it suffers from the PAM limitation³⁹. The oncogene SNP only has a small sequence window to probe, the traditional PAM, TTTV, could not cover all the SNPs. Therefore, we tested whether the AmCas12a can distinguish the SNPs without a traditional PAM. (Fig. 7a) The oncogene mutants, KRAS c.34 G > T (G12C), did not contain the available TTTV in the adjacent sequences (Fig. 7b). Among the Cas12a proteins that have undergone PAM preference testing, AmCas12a, EvCas12a_2, CAGGCas12a, and RbrCas12a_1 showed potential for recognizing the G12C mutation. The results revealed that AmCas12a exhibited the best performance (Supplementary Fig. 20). We designed the crRNA targeting the SNP (Fig. 7b). According to the fluorescence intensity, we selected the crRNAs inducing the strongest signals, i.e., crRNA 1 for the KRAS mutant (Fig. 7c). The AmCas12a can detect ten copies of the KRAS mutant (Fig. 7d). Furthermore, we diluted the target mutant and evaluated the sensitivity of detection. The AmCas12a can even distinguish 0.1% KRAS mutant in the wild-type gene background, which is more sensitive than the Sanger sequencing (Fig. 7e, f).

In this section, we first present the experimental settings of our study. Next, we evaluate our method by comparing it with state-of-the-art methods and report both quantitative and qualitative results. Finally, we explore the impact of basic blocks and various modes on SR performance. Additionally, we assess the effectiveness of our method on edge devices.

Experimental setup

Datasets and Metrics: The dataset contains 25,892 valid infrared image samples, all acquired using five self-developed infrared imaging systems equipped with a 640(times)512 uncooled IRFPA. These 640(times)512 images serves as Ground Truth HR references during the training phase. The corresponding LR images are generated using non-overlapping average pooling operations (e.g., a 2(times)2 kernel for (times)2 downsampling) instead of bicubic interpolation, motivated by our proposed readout circuit structure prior that models the physical infrared imaging process characterized by row-wise scanning and column-wise readout. Average pooling is more consistent with this mechanism and better preserves the spatio-temporal correlations in infrared images, effectively avoiding the artifacts and distortions commonly introduced by interpolation-based downsampling. Furthermore, due to the high cost and limited availability of megapixel-level infrared imaging systems, it is not feasible to obtain real infrared images at a resolution of 1280(times)1024 that are perfectly aligned with the corresponding low-resolution counterparts. To evaluate SR performance at this scale, we generate pseudo HR references using the Upscayl, an image upscaling tool based on an open-source large-scale AI model. Although originally designed for natural image enhancement, Upscayl can reconstruct plausible high-frequency textures that serve as reasonable references for evaluating the quality of our reconstructions images. This approach facilitates the assessment of the performance of our method in the absence of true HR infrared image. The core infrared detectors of all imaging devices have the following key performance parameters: a pixel size of 17 µm, a 640(times)512 focal plane array, a noise equivalent temperature difference (NETD) of 25 mK, a time constant of 8 ms, a frame rate of 50 Hz, and a response wavelength range of 8-14 µm. We utilize 2,500 images to evaluate the performance of different approaches. The supplementary file presents the infrared image datasets utilized in this work, including images acquired from a commercial cooled infrared imaging system, synthetically generated high-resolution infrared images, and validation data obtained from a self-developed uncooled infrared detector. Average peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM) are employed as evaluation metrics.

Implementation Details: Different configurations of our proposed EIRSR are presented in Table 1. Data augmentation includes horizontal/vertical flips and random rotations of 90(^circ), 180(^circ), and 270(^circ). The kernel sizes of all convolutions are limited to 3 and 1. The batch size is set to 32, and the size of each ({{I}_{LR}}) is set to (text {48}times text {48}) during the training phase. We employ Adam³⁸ optimizer to train the model with ({{beta }_{text {1}}}=text {0}text {.9}), ({{beta }_{text {2}}}=text {0}text {.999}). The initial learning rate is set to (text {5}times text {1}{{text {0}}^{text {-4}}}) and decays following a cosine learning rate. ({{{L}}_{SR}}) is used to optimize the model over (text {5}times text {1}{{text {0}}^{text {5}}}) iterations, with the initial relaxation factor (alpha) set to 0.2 and is halved every (text {1}times text {1}{{text {0}}^{text {5}}}) iterations. Our method is implemented in PyTorch, and all experiments are conducted on a single GeForce RTX 4090 GPU.

Comparison with state-of-the-art methods

To evaluate the effectiveness of EIRSR, we compare it with several advanced efficient SR methods, including RFDN², LatticeNet²⁸, SwinIR¹⁷, ELAN¹⁸, BSRN³⁰, ESRT²⁹, HAT¹⁹, LKDN³¹, RGT²⁰, PLKSR³², and SeemoRe¹⁴. Table 2 presents the quantitative comparison results of PSNR and SSIM for (times text {2/}times text {3/}times text {4}) upscale factors, along with the number of parameters and Multi-Adds. At scale (times text {2}), the number of parameters and Multi-Adds for our method are 2.54(%) and 2.9(%) of those for the second-ranked method, and 2.77(%) and 2.56(%) of those for the third-ranked method. At scale (times text {3}), the number of parameters and Multi-Adds for our method are 2.74(%) and 3.84(%) of those for the second-ranked method, and 2.51(%) and 3.72(%) of those for the third-ranked method. At scale (times text {4}), the number of parameters and Multi-Adds for our method are 2.54(%) and 3.15(%) of those for the second-ranked method, and 2.77(%) and 2.91(%) of those for the third-ranked method. The results demonstrate that EIRSR, which utilizes a CNN-Transformer structure, surpasses previous leading models in PSNR and SSIM metrics. Comparison results reveal that, for infrared images, Transformer architecture like as HAT¹⁹, RGT²⁰, SwinIR¹⁷, and ESRT²⁹ outperform the convolutional architecture. This advantage may be attributed to the relationship between the infrared imaging process and the feature split of the Transformer.

In Fig. 5, we present a visual comparison of EIRSR and other efficient methods on (times text {2/}times text {3/}times text {4}). It is evident that the HR images reconstructed by EIRSR exhibit more accurate texture details, particularly along the edges. Compared to Ground Truth (GT) images at scale (times text {2/}times text {3}), our reconstructed images show overall smoothness with no obvious artifacts on the edges, which demonstrates superior visual quality and can be attributed to the incorporation of an image enhancement regularization control term in the loss function. The difference maps at scales (times text {2}) and (times text {3}) show that our method has the smallest discrepancy with the GT images in terms of both overall structure and fine details. The difference map at scale (times text {4}) shows that our method does not generate obvious artifacts. It is worth noting that at a scale of (times text {4}), most methods do not perform well due to the inherent lack of details in infrared images, as the algorithm cannot generate non-existent details. All comparative experiments demonstrate the effectiveness of our method. Furthermore, it is important to emphasize that Transformer-based methods outperform CNNs in infrared image SR, as demonstrated by HAT¹⁹ and RGT²⁰. Although the ViT block in both methods employs window-based MHSA, it is essential to recognize that ViT is fundamentally linked to the mechanism of infrared imaging.

We present a comparison of the typical SR methods IRSRMamba¹² and PSRGAN⁴ on infrared images, as illustrated in Fig. 6. Our method demonstrates superior performance compared to the current SR methods for infrared images, particularly in terms of image details and visual perception. Furthermore, the differential analysis through error mapping demonstrates that our reconstructed images maintain the closest structural fidelity to the ground truth references in terms of global feature consistency.

Ablation study

In this section, we perform ablation studies on important design elements in the proposed EIRSR to explore the impact of different blocks on infrared image reconstruction performance. Table 3 shows the results.

Effectiveness of CCB and its internal component USCAB. We conduct an ablation study to evaluate the effectiveness of CCB and its internal component USCAB, as shown in Table 3. By masking the CCB, EIRSR reduces to using only the RCTB component. In parallel model, PSNR decreases by 1.42(%) and SSIM decreases by 1.44(%) (see (#)1 and (#)5), while in serial mode, PSNR decreases by 1.65(%) and SSIM decreases by 1.09(%) (see (#)7 and (#)11). This finding confirms that, in the hybrid architecture, the local features extracted by CCB are crucial for establishing effective long-range dependencies between pixels using the RCTB, ultimately improving the network’s performance. Furthermore, we mask USCAB to assess its impact on performance. In parallel model, PSNR degrades by 0.73(%) and SSIM degrades by 1.06% (see (#)1 and (#)2), whereas in the serial mode, PSNR degrades by 0.89(%) and SSIM degrades by 0.25(%) (see (#)7 and (#)8). The USCAB module operates interactively in both channel and spatial dimensions, facilitating the extraction of potential correlations between pixel locality and feature channels. This dual interaction significantly contributes to the improvement of performance in the SR task, especially when integrated with RCTB for global context modeling. These results validate that the CCB, especially when integrated with the USCAB, substantially improves the model’s capacity to extract local features and reinforce local contextual representations. Such localized enhancements are essential for facilitating the global dependency modeling in RCTB, thereby improving both reconstruction quality and overall SR performance in infrared imaging.

Effectiveness of RCTB. We compare the impact of RCTB on SR performance in three cases. When RCTB ((times) rows) operates only in columns, in parallel model, the PSNR degrades by 0.78(%) and the SSIM degrades by 1.1(%) (see (#)1 and (#)3). In serial mode, PSNR degrades by 1.08(%) and SSIM degrades by 0.44(%) (see (#)7 and (#)9). Conversely, when RCTB ((times) cols) operates solely in rows, the parallel model shows a PSNR degradation of 1.32(%) and an SSIM degradation of 1.3(%) (see (#)1 and (#)4), while in serial mode, the PSNR degrades by 1.56(%) and the SSIM degrades by 1.09(%) (see (#)7 and (#)10). In the absence of RCTB, the parallel model exhibits a PSNR degradation of 4.7(%) and an SSIM degradation of 4.12(%) (see (#)1 and (#)6), and in serial mode, the PSNR degrades by 5.0(%) and the SSIM degrades by 3.24(%) (see (#)7 and (#)12). These comparative results demonstrate that applying the transformer to either rows or columns alone is less effective than applying it to both.

The effectiveness of the RCTB is crucial for enhancing the performance of infrared image SR. RCTB is specifically designed to capture long-range dependencies by modeling correlations across rows and columns, a feature particularly important for infrared images where spatio-temporal relationships are critical for accurate reconstruction. Ablation experiments demonstrate that applying RCTB to both rows and columns yields superior performance compared to applying it to rows or columns individually. This indicates that fully applying RCTB enables it to capture interdependencies between pixels across both dimensions, thereby providing comprehensive image features and improving SR quality. Notably, masking RCTB results in significant performance degradation. For instance, in the parallel mode, the PSNR degrades by 4.7(%) and SSIM by 4.12(%), and in the serial mode, PSNR degrades by 5.0(%) and SSIM by 3.24(%). This sharp performance drop underscores the importance of RCTB in capturing global context, solidifying its role as an essential component of our framework. Furthermore, the combination of RCTB and CCB enables a comprehensive feature extraction approach. CCB is responsible for efficiently extracting high-frequency local details, while RCTB handles the long-range global dependencies. By integrating these two modules, our model harnesses their complementary strengths: the CCB enhances local feature representation, whereas the RCTB captures global spatio-temporal correlations grounded in the infrared imaging process. These results underscore the critical role of the RCTB design, which is inspired by the readout characteristics of IRFPA detectors, as essential for performance improvement and not substitutable by conventional designs relying solely on image content or network architecture. The synergy between CCB and RCTB allows the model to capture both fine textures and global coherence, which is critical for infrared image reconstruction. Another critical aspect of RCTB’s design is that it is based on the IRFPA readout circuit. The IRFPA circuit operates by scanning infrared images row-by-row and column-by-column, and RCTB is designed based on this prior knowledge. By incorporating this prior, RCTB effectively models spatio-temporal correlations between pixels across rows and columns, aligning with the structure of the IRFPA readout circuit. This design allows the network to capture semantically richer features from infrared images, significantly improving performance in infrared SR tasks.

In summary, the integration of RCTB with CCB significantly enhances the model’s ability to capture both local and global features. By leveraging the IRFPA readout circuit’s characteristics, RCTB further improves the model’s ability to handle complex dependencies in infrared images, establishing it as a crucial component for high-performance infrared super-resolution.

As illustrated in Fig. 7, we analyze the cosine correlation between rows and columns in features based on the EIRSR-parallel model in Table 3, revealing several noteworthy findings. First, a comparison between Fig. 7a and Fig. 7b demonstrates that integrating CCB enhances RCTB’s ability to capture high-level semantic information by preserving both row and column correlations. This enhancement is attributed to CCB’s efficient extraction of local features, which supports the long-range dependencies modeled by RCTB across rows and columns. In contrast, when CCB is used without USCAB, as shown in the feature correlation maps of Fig. 7b, the model’s ability to maintain row and column correlations within RCTB is significantly reduced, underscoring the critical role of USCAB in preserving these dependencies. Furthermore, as shown in Fig. 7a and Fig. 7c, the absence of row splitting in RCTB leads to a reduction in both row and column correlations, with row correlations being more significantly weakened. This highlights the importance of the row-splitting mechanism in RCTB, which is essential for preserving strong dependencies between pixels across different rows. Similarly, the comparison between Fig. 7a and Fig. 7d demonstrates that removing column splitting results in a reduction in both row and column correlations, with column correlations experiencing a more pronounced decline. This underscores the importance of column splitting in RCTB for capturing global dependencies between columns, a critical factor for accurately modeling pixel relationships in infrared images. Additionally, the comparison of the red and green boxes in Fig. 7e reveals that RCTB demonstrates a superior ability to model row and column correlations compared to CCB alone. This further emphasizes RCTB’s unique capacity to capture long-range dependencies and spatial-temporal correlations across rows and columns, a key factor in enhancing SR performance. Overall, the results shown in Fig. 7 indicate that column correlations are more prominent than row correlations. This phenomenon can be attributed to the architecture of the IRFPA readout circuit, where each column of pixels shares a single readout channel, resulting in stronger column-wise correlations. By leveraging this inherent structure, RCTB enhances the modeling of interdependencies across rows and columns, leading to significant improvements in infrared image SR performance.

Loss function. We introduce a regularization control term into the loss function and dynamically adjust this loss function using a relaxation factor (alpha) during training, which yields interesting results, as illustrated in Fig. 8. Sub-image (a) shows that without the introduction of control terms, our method produces more noise than GT images, resulting in unsatisfactory outcomes. In sub-image (b), we observe that with the introduction of the control item, ({{tau }_{prep}}left( I_{HR}^{i} right))represents (I_{HR}^{i}) processed by guided filtering and demonstrates superior performance compared to the GT image processed directly by guided filtering. Furthermore, in sub-image (c), where ({{tau }_{prep}}left( I_{HR}^{i} right)) refers to the application of guided filtering and detail enhancement (Laplacian sharpening) on (I_{HR}^{i}), our method outperforms the direct application of guided filtering and detail enhancement on GT images in the terms of image detail. The computational profiling conduct on the RK3588 Core Board reveals clear temporal characteristics: standalone guided filtering and Laplacian sharpening operations require 18.79 ms and 6.701 ms respectively under single-threaded mode. Our integrated architecture, which combines these preprocessing operators, demonstrated 37.815 ms processing latency. Compared to the conventional sequential approach (18.79 ms + 6.701 ms + 37.815 ms), the proposed end-to-end implementation achieves a 40.27(%) reduction in total execution time. The experimental analysis of the loss function encourages us to investigate the integration of the infrared image preprocessing algorithm into the network in future research by incorporating the control term into the loss function, with the aim of reducing the computational cost associated with infrared imaging system preprocessing and minimizing overall processing latency.

SR comparison under different readout modes

To validate the effectiveness of the proposed spatio-temporal readout prior, we conduct a comparative experiment using infrared images acquired by a self-developed infrared imaging system and a commercial infrared imaging system, operating in rolling shutter and global shutter readout modes, respectively. The IRFPA in the self-developed system operates in a rolling shutter readout mode, which performs row-wise scanning and column-wise readout, in contrast to the commercial system that adopts a global shutter readout mode. The commercial system features a pixel size of 15 µm, a 640(times)512 focal plane array, and a noise equivalent temperature difference (NETD) (le) 17 mK. A total of 2,500 images are used for validation in both the rolling shutter and global shutter imaging systems. The average results of the quantitative comparison are summarized in Table 4. As shown in Table 4, the spatio-temporal readout prior-based method achieves PSNR improvements of 6.94(%) and 9.65(%), and SSIM improvements of 2.65(%) and 6.68(%) on the rolling shutter imaging system, compared to its performance on the global shutter system, under (times)2 and (times)4 upscaling factors, respectively.

Representative comparison images are selected from different imaging systems under (times)2 and (times)4 upscaling factors, and their corresponding PSNR, SSIM, and difference maps are computed, as illustrated in Fig. 9. As shown in Fig. 9, the proposed method achieves better performance on the self-developed imaging system that incorporates spatio-temporal readout priors, whereas its effectiveness is less pronounced on the global shutter imaging system, which lacks such priors. Specifically, under the (times)4 SR scenario, it fails to reconstruct the vertical structural components of the glass curtain wall on the global shutter imaging system, leading to a reconstruction that retains only the horizontal stripe patterns, with the vertical features entirely absent. The comparative results across different imaging modes support the effectiveness of the proposed method that incorporates spatio-temporal readout priors. Moreover, these findings imply that accounting for hardware-level imaging characteristics can be beneficial to the performance of SR tasks. This further underscores the design specificity of our network design for row-wise scanning and column-wise readout IRFPAs, in which spatio-temporal readout priors play a crucial role in guiding the reconstruction process. While this specificity contributes to significant performance improvements on row-column scanned systems, it also underscores the need to adapt our framework for other imaging sensor architectures–such as global shutter or event-based imaging systems–where differing physical imaging process and spatio-temporal dynamics may necessitate alternative modeling approaches.

Edge device deployment

To validate the effectiveness of our model on edge devices, we optimized EIRSR-T as EIRSR-T-opt and evaluated it on an edge inference device: RK3588 Core Board, an embedded system-on-module (SoM) from Rockchip, which features three integrated NPU cores. We assessed the performance of models at a scale of (times text {2}) in single-process mode, utilizing 16-bit floating point precision during inference. For each input image size, we executed the models for two hours to avoid the warm-up effect, the results are presented in Table 5. As the size of the input image increases, there is a corresponding rise in power consumption, memory usage, memory read/write operations, and runtime for the model. In single-threaded mode, models with low power consumption can be deployed to edge devices. However, our optimized model cannot achieve real-time processing speeds in the single-threaded mode with an input size of 1280(times)1024. In this case, real-time SR for large images can be achieved through multi-core and multi-threaded processing, but this approach significantly increases the resources consumption of edge devices. Table 5 illustrates that memory usage and memory read/write operations are significant bottlenecks that limit the model’s performance. To facilitate deployment on edge devices, we have summarized several guidelines for model optimization. Specifically, the following strategies are recommended: consider operator fusion whenever possible, implement weight sharing during model quantization, ensure that the number of feature channels is a multiple of four, adopt a general 3(times)3 convolution kernel, reduce the number of heads in MHSA, maximize the split of row and column features, and utilize operators that are optimized for the specific hardware platform.

Motivation and applicability of the hardware prior

Our method is inspired by the row-wise scanning and column-wise readout mechanism of our self-developed uncooled IRFPA detectors. This readout mechanism introduces inherent temporal correlations among row pixels and spatial correlations among column pixels during the image formation process, both of which are explicitly exploited in our model design. To capture these correlations, we propose the RCTB, which applies self-attention separately along the row and column dimensions. This design aligns closely with the physical imaging mechanism and enables the network to effectively capture pixel-level dependencies introduced by the readout circuitry, thereby yielding improvements in SR performance, as demonstrated in our ablation studies. While our method is tailored to IRFPAs exhibiting such readout characteristics, this class of imaging sensors is widely deployed in low-power, cost-sensitive, and edge-oriented infrared imaging systems. Therefore, the proposed method has considerable potential for practical deployment.

We also acknowledge that the method is not directly applicable to global shutter mode imaging sensors, which lack the row-column spatio-temporal dependencies leveraged in our design. As demonstrated in the “SR Comparison under Different Readout Modes,” the method exhibits reduced effectiveness. Nonetheless, the principle of incorporating hardware-level priors into network design can be extended to other imaging architectures through appropriate modifications, which we aim to explore in future work. Compared with previous infrared SR methods that focus solely on images or networks, our approach is the first to integrate the imaging circuitry structure priors into the network architecture. This allows for more efficient modeling of spatio-temporal correlations that are consistent with the hardware, leading to enhanced reconstruction fidelity and robustness in the infrared SR task.

Introduction

As the core technical force of China’s healthcare service system, the optimal allocation of pharmacists’ human resources is not only crucial for ensuring medication safety and improving healthcare quality but also serves as a key measure in achieving the goal of health equity. The World Health Organization emphasized in its Global Strategy for Human Resources for Health¹ that the global density of healthcare personnel, including pharmacists, is generally insufficient—particularly in low-income countries and rural areas—leading to significant bottlenecks in medication accessibility and safety. In China, adverse drug reactions (ADRs) result in approximately 2.5 million hospitalizations annually, with 500,000 classified as serious ADRs and around 190,000 deaths.² On the other hand, studies have demonstrated that optimization of patients’ medication regimens can significantly reduce healthcare costs.³ Pharmacists thus play an essential role in ensuring rational and safe medication use while simultaneously reducing healthcare expenses.

Since the initiation of the 2009 healthcare reform in China, the human resources for pharmacists have grown considerably, with their numbers increasing from 342,000 in 2009 to 531,000 in 2022—a 55.26% increase. However, the growth rate of pharmacists’ human resources and their proportion within healthcare professionals remain notably lower than other healthcare categories. Additionally, a serious shortage persists in healthcare institutions, where only 45% of the nation’s pharmacists are responsible for approximately 73% of pharmacy service workloads.⁴ Furthermore, due to regional disparities in development, China faces challenges in achieving balanced allocation of pharmacists’ human resources across different regions and healthcare tiers. At the same time, there have been significant changes in the factors affecting the distribution of human resources for health, changes that have, to some extent, exacerbated regional differences in the allocation of health resources. There has been a gradual shift from a single administratively-dominated to a diversified pattern, with the early distribution being heavily influenced by administrative factors, and with the development of the market economy, economic and geographic factors have gradually become dominant, with developed regions attracting large numbers of medical personnel and less developed regions experiencing a relative lack of resources. In addition, the aging population and rising burden of chronic diseases further exacerbate these challenges, creating a widening gap in the demand for pharmacy services. Addressing the uneven regional and hierarchical distribution of pharmacists’ human resources has thus become a critical research issue for both the Chinese government and academic scholars.

Since 2009, the Chinese government has undertaken a new round of healthcare system reform, with primary objectives centered around reducing urban-rural and regional disparities in health resources and promoting equalization of basic public health services. In 2015, efforts were intensified to strengthen the primary healthcare system and implement a tiered diagnosis and treatment system that aims to distribute healthcare resources more evenly across various healthcare institution levels.⁵ The emphasis on equalization of basic public health services was reaffirmed in the Healthy China 2030 strategic outline. By March 2023, the government released the Opinions on Further Improving the Medical and Health Care Service System, which detailed staffing standards for specialized public health institutions, improved pharmacy service quality, and enhanced drug supply security measures.

Considerable scholarly attention has been directed toward equitable allocation of healthcare resources. Studies have analyzed fairness at national and regional levels using approaches such as Lorenz curves, Gini coefficients, and Theil indices, revealing disparities in healthcare systems both nationally and within specific provinces.^6–11 Chen et al¹² posited that these disparities in China primarily stem from imbalances in economic development and population distribution. Other researchers have investigated the efficiency of health resource allocation through methods like Data Envelopment Analysis (DEA) and Malmquist models, uncovering inefficiencies that hinder the service capacity of healthcare human resources.^13,14

In terms of pharmacists’ human resource allocation, Ni et al¹⁵ applied metrics such as Gini coefficients and Theil indices to explore manpower distribution, while Wei et al¹⁶ employed the Vector Autoregression (VAR) model to identify factors affecting pharmacists’ growth, finding that hospital visit numbers, per capita healthcare expenditure, and per capita GDP significantly influence the development of pharmacists’ resources. Wu et al¹⁷ constructed two models—the pharmacy service demand method and the economic indicator method—to estimate pharmacist demand in healthcare institutions. Zhou et al¹⁸ utilized agglomeration models to assess the equilibrium of pharmacists’ human resource allocation, while Qiao et al¹⁹ and Zhu et al²⁰ employed spatial analysis methods to examine the geographic distribution characteristics of pharmacists. He et al²¹ further applied VAR models and Gini coefficients to analyze disparities in pharmacists’ human resources in Tianjin, identifying significant spatial autocorrelation in pharmacist distribution.

Although these studies have contributed valuable insights, they primarily focus on the fairness and spatial effects of pharmacists’ human resource allocation, with limited exploration of regional and tiered differences. Moreover, existing studies largely treat pharmacist allocation on a single-tier basis, neglecting disparities across different healthcare system levels (eg, tertiary hospitals versus primary healthcare institutions). Failure to consider these tiered differences may lead to oversimplified conclusions and hinder accurate assessments of pharmacist shortages across healthcare tiers.

In summary, while previous literature has addressed fairness and spatial distribution issues in pharmacists’ human resource allocation in China, there remains a gap in research concerning regional and tiered disparities, as well as in-depth analyses of spatial relationships. To address these shortcomings, this study employs the Dagum Gini coefficient, kernel density analysis, and spatial β-convergence models to comprehensively examine regional differences, dynamic evolution, and convergence trends in pharmacists’ human resource allocation across two tiers: hospitals and primary healthcare institutions (PHC). This study provides valuable insights into the distribution and allocation of pharmacist resources in China, underscoring the critical need for targeted strategies to address disparities and support primary healthcare institutions in achieving equitable access to pharmacist services.

Materials and Methods

Data Sources

In this paper, pharmacists per 1,000 population are used as the base metric for all subsequent calculations. Specifically, the total number of pharmacists is derived from the China Health and Health Statistics Yearbook (2018–2023), the China Health and Family Planning Statistics Yearbook (2016–2017), and the China Health Statistics Yearbook (2013–2015). Population data, as well as data for control variables, are obtained from the China Statistical Yearbook (2013–2023).

The study period is defined as 2012–2022, and the spatial scope focuses on 31 provinces, autonomous regions, and municipalities directly under the central government in mainland China (excluding Hong Kong, Macao, and Taiwan). For regional analysis, the paper employs the commonly accepted division of China into three major economic zones: eastern, central, and western regions. Specifically: Eastern Region: Beijing, Tianjin, Hebei, Liaoning, Shanghai, Jiangsu, Zhejiang, Fujian, Shandong, Guangdong, and Hainan; Central Region: Shanxi, Jilin, Heilongjiang, Anhui, Jiangxi, Henan, Hubei, and Hunan; Western Region: Inner Mongolia, Guangxi, Chongqing, Sichuan, Guizhou, Yunnan, Shaanxi, Gansu, Qinghai, Ningxia, Xinjiang, and Tibet.

Methods

Dagum Gini Coefficient

The traditional Gini coefficient and Dagum Gini coefficient both measure inequality, but differ in analytical capacity. While the traditional approach effectively quantifies income inequality within homogeneous populations, it cannot identify disparity sources in complex cross-regional contexts. In contrast, the Dagum Gini coefficient provides methodological improvements that enable precise decomposition of regional imbalances. This advanced approach proves particularly valuable when examining inequalities across geographical regions, demographic subgroups, or national boundaries, as it facilitates nuanced analysis of contributing factors.²² In this paper, Dagum’s Gini coefficient and its decomposition method are utilized to analyze the disparities in the allocation of pharmacists’ human resources across various levels in China. This method decomposes the overall disparity into three components: intra-regional disparity, inter-regional disparity, and hypervariable density (hypervariable density quantifies the contribution of resource distribution overlaps between differentially ranked regions to the aggregate imbalance, capturing how disparities between high- and low-resource provinces influence overall inequality when their endowment distributions intersect).²³ The specific formula is shown below:

(1)

(2)

In Equation (1), G represents the overall Gini coefficient, where a larger value indicates greater disparities in the allocation of pharmacists’ human resources. The parameters are defined as follows: mm denotes the number of regional divisions, is the number of provinces in the ath subregion, is the number of provinces in the bth subregion, represents the pharmacists’ human resource allocation in the ith province within the aa-th subregion, refers to the allocation in the rth province within the bth subregion, n is the total number of provinces, and is the mean value of pharmacists’ human resource allocation across all provinces; In Equation (2), , , and , represent the contributions of intra-region variation, inter-region variation, and contribution of hypervariable density, respectively. These components are used to quantify the contributions of intra-regional disparities, inter-regional disparities, and cross-regional interactions to the overall variation in the allocation of pharmacists’ human resources.²⁴ For details on the specific decomposition of the formula, refer to the related study by Dagum.²⁵

Kernel Density Estimation

The kernel density estimation method was employed to analyze the distribution dynamics and trends in the allocation of pharmacists’ human resources in China. The formula is:

(3)

Equation (3) defines the parameters as follows: N represents the number of samples, h denotes the bandwidth, x is the mean value, and refers to the kernel function. The kernel density curve provides critical insights into the research subject during the study period, including its distribution position, trend, and extensibility. Specifically, the distribution position reflects the growth of pharmacists’ human resource allocation, where higher values indicate a higher allocation level. The distribution trend reveals disparities in pharmacists’ human resource allocation, as indicated by the height of the wave peaks, while the number of peaks identifies the degree of polarization in the allocation. Finally, the distribution extensibility captures the gap between the highest and lowest allocation levels among provinces, offering insights into regional differences in pharmacists’ resource distribution.²⁶

Global Moran’s Index

The global Moran index can judge the spatial distribution characteristics between variables from the significance level, in which the positive significance sign implies that there is a centralized distribution characteristics between variables, and the negative significance characteristics represent the dispersed characteristics. The size of the global Moran index can also judge the size of the spatial link between variables, the closer the positive value is to 1, the stronger the concentration, the closer the negative value is to −1, the stronger the dispersion.²⁷ The formula is:

(4)

In Equation (4), is the observed value of region i, is the observed value of region j, and is the spatial weight matrix with spatial adjacency of 1. >0 is a positive spatial correlation, which indicates that there is a significant positive spatial correlation of pharmacists’ human resource allocation.²⁸

Spatial β Convergence Model

The concept of traditional convergence originates from the neoclassical school of economic growth theory and encompasses three primary types: absolute convergence, conditional convergence, and club convergence. Absolute convergence assumes that all regions possess identical initial endowments, whereby less developed regions exhibit higher growth rates, allowing them to catch up with more developed regions over time and eventually achieve a relatively stable state of development. Conditional convergence expands on absolute convergence by incorporating variations in initial endowments, acknowledging that differences in these factors can significantly affect the rate and likelihood of convergence.²⁹ However, in the context of the spatial distribution of pharmacists’ human resources across regions in China, the traditional β-convergence model fails to account for spatial correlation effects that may exist between these regions. To address this limitation, this study integrates spatial correlation into the traditional β-convergence model, drawing on the methodology proposed by Castellanos-Sosa et al.³⁰ Specifically, the Spatial Autoregressive (SAR) model, Spatial Error Model (SEM), and Spatial Durbin Model (SDM) are employed to conduct a spatial convergence analysis. The formula is:

(5)

(6)

(7)

In equations (5–7), is the observation of the ith region in period t, represents an independent and identically distributed residual term, and is the error term with spatial autocorrelation. is the spatial autoregression coefficient, is spatial autocorrelation coefficient of error item, and is the spatial lag coefficient of . In addition, is the element in spatial weight matrix, denotes the spatial interaction of , and represent province and time-fixed effects. is the coefficient of control variable. is the spatial lag coefficient of control variable.³¹ is the control variable, The main components are the level of economic development (GDP), the level of investment in health (IVH), the urbanization rate (URB), the population density (POP), and the level of government intervention (GOV),^32–34 The definitions and data sources for these variables are shown in Table 1. if β<0, it means that there isβconvergence.

Table 1 Definition and Descriptive Statistics of Control Variables

Results

Temporal Dynamics of Pharmacists’ Human Resource Allocation

As shown in Figure 1, the human resources of pharmacists in China experienced substantial growth between 2012 and 2022, with the total number increasing from 358,511 to 504,401, reflecting an average annual growth rate of 3.4%. At different institutional levels, pharmacist resources in Chinese hospitals increased from 231,249 to 332,849 during the same period, with an average annual growth rate of 3.7%. Meanwhile, pharmacist resources in PHC rose from 127,262 to 171,552, yielding an average annual growth rate of 3%. Furthermore, the disparity between hospitals and PHCs exhibited a fluctuating upward trend, rising from 1.82 in 2012 to a peak of 2.03 in 2018, before slightly declining to 1.94 between 2019 and 2022. This shift can be attributed to two post-pandemic policy developments in China. First, governments at all levels substantially increased investments in health resources, with particular emphasis on primary healthcare institutions. Second, nationwide implementation of the medical consortium model enhanced the integration of grassroots facilities, promoting the redistribution of pharmacist human resources to primary care settings and ultimately reducing inter-tier disparities across healthcare organizations.

Figure 1 Regional Differences in Pharmacists’ Human Resources in China.

Overall, despite the relatively rapid development of pharmacist human resources in China, the growth rate and total quantity remain significantly higher at the hospital level compared to PHCs. The disparity between pharmacist resources across different healthcare tiers has widened. However, the outbreak of the COVID-19 pandemic has notably altered China’s health resource allocation patterns and affected the trajectory of this discrepancy.

From a regional perspective, the total pharmacist resources across different tiers were substantially higher in the eastern region, followed by the central and western regions. Specifically, pharmacist resources in the eastern region increased from 167,161 to 236,291, reflecting an average annual growth rate of 3.5%. The central region saw an increase from 106,059 to 132,027, with an average annual growth rate of 2.2%, while the western region demonstrated the most rapid growth, expanding from 85,291 to 136,083, with an average annual growth rate of 4.8%. Notably, the western region surpassed the central region in resource quantity in 2021. Regarding tier disparities, all three major regions displayed an initial rise followed by a decline, albeit with differing inflection points: the eastern region peaked in 2016, the western region in 2017, and the central region in 2019. On average, the eastern and western regions maintained a value of 1.97 from 2012 to 2022, whereas the central region recorded a higher average value of 2.12.

Regional Differences and Sources of Differences in Human Resource Allocation for Pharmacists in China

Regional Differences in Pharmacist

As shown in Figure 2A, the Gini coefficient for hospital pharmacist human resource allocation in China decreased from 0.147 in 2012 to 0.103 in 2022, demonstrating a reduction in regional disparities over the study period. Regionally, the eastern, central, and western regions all experienced declines in their hospital pharmacist Gini coefficients: 0.149, 0.083, and 0.14 in 2012 compared to 0.13, 0.053, and 0.081 in 2022, with average annual reduction rates of −1.3%, −4.3%, and −5.2%, respectively. Notably, the eastern region consistently had a higher Gini coefficient than the national average. The central region consistently had the smallest Gini coefficient, whereas the western region exhibited the fastest decline.

Figure 2 Changes in the Gini Coefficient of Pharmacists’ Human Resource Allocation in China. (A) Changes in the Gini coefficient of human resource allocation among hospital pharmacists in China, (B) Changes in the Gini coefficient of human resource allocation among PHC pharmacists in China.

As shown in Figure 2B, the Gini coefficient for pharmacist human resource allocation in China’s PHCs displayed a fluctuating downward trend, dropping from 0.207 in 2012 to 0.192 in 2021, before slightly rebounding to 0.195 in 2022. This indicates that, while regional disparities in PHC pharmacist allocation have decreased overall, the magnitude of the decline was smaller compared to hospitals. Regionally, significant heterogeneity was observed in PHC pharmacist human resource allocation. The eastern region exhibited a unique upward trend, with its Gini coefficient rising from 0.171 in 2012 to 0.189 in 2022, reflecting an annual growth rate of 1.0%. In contrast, the central and western regions demonstrated declining trends, with Gini coefficients decreasing from 0.141 and 0.246 in 2012 to 0.095 and 0.182 in 2022, respectively, corresponding to average annual reduction rates of −3.9% and −3.1%.

Overall, regional disparities in pharmacist human resource allocation in both hospitals and PHCs have generally decreased. However, disparities in PHCs remain larger than those at the hospital level, particularly in the eastern and western regions. Economic development disparities between regions are a major contributing factor. For example, in 2012, the western region’s weaker economic base made it challenging to allocate sufficient pharmacist resources to PHCs, compounded by lower salaries leading to workforce migration to wealthier areas within the region. Conversely, the eastern region’s economic growth (eg, in Zhejiang and Jiangsu provinces in the Yangtze River Delta) provided greater financial capacity for PHC expansion, increasing disparities between provinces within the region.

Interregional Differences in Pharmacist

Figure 3A illustrates the trend of inter-regional differences in human resource allocation for hospital pharmacists in China from 2012 to 2022. Over the study period, all three major inter-regional gaps (Eastern & Western, Eastern & Central, and Central & Western) exhibited a decreasing trend. Between 2012 and 2021, the largest inter-regional gap was observed between the Eastern and Western regions, while in 2022, this shifted to the disparity between the Eastern and Central regions.

Figure 3 Interregional Differences in Pharmacists’ Human Resource Allocation in China. (A) Interregional differences in human resource allocation among hospital pharmacists in China, (B) Interregional differences in human resource allocation among PHC pharmacists in China.

Figure 3B 2highlights the trend of inter-regional differences in the human resource allocation of pharmacists in PHCs in China over the same period. Both the Eastern & Western and Central & Western disparities displayed a downward trend, whereas the Eastern & Central disparity showed a significant upward trajectory. During 2012–2013, the Eastern & Central disparity was the smallest among the three pairs of regions. From 2014 to 2017, this disparity became the second largest, while in 2018–2022, it emerged as the most pronounced inter-regional gap.

Sources of Differential Contributions by Pharmacists

Figure 4A hows the sources of contribution to the human resource allocation gap for hospital pharmacists in China from 2012 to 2022. During the study period, disparities in hospital pharmacist allocation were primarily driven by inter-regional disparity, which remained the largest and most stable source of variation, followed by intra-regional disparity and hypervariance density, which was consistently the smallest contributor. Specifically, inter-regional disparities exhibited a fluctuating downward trend from 42.52% in 2012 to 34.76% in 2020, followed by an upward trend from 2021 to 2022, reaching 41.10%. This pattern is consistent with the sources of differences observed in 2012.

Figure 4 Sources of Gini coefficient variance contributions. (A) Sources of differences in hospital Gini coefficients, (B) Sources of differences in PHC Gini coefficients.

Figure 4B resents the sources of contribution to the variance in human resource allocation for PHC pharmacists. During the study period, the disparities in PHC pharmacist allocation showed substantial changes. In general, intra-regional gaps remained relatively stable, while inter-regional disparities and hypervariance density fluctuated significantly between 2012 and 2016, alternately becoming the largest source of variation. From 2016 to 2022, inter-regional disparities exhibited a consistent upward trend, increasing from 33.01% to 43.44%, while hypervariance density declined from 34.62% to 26%.

Overall, the disparities in pharmacist human resource allocation for both hospitals and PHCs followed a broadly similar trend: a gradual increase in inter-regional disparity, a decrease in hypervariance density, and relatively stable intra-regional disparity.

Kernel Density Curve Estimates

Kernel Density Profile Analysis of Hospital Pharmacists

Figure 5 illustrates the dynamic evolution of hospital pharmacists’ human resource allocation distribution across China and its three major regions during the study period. First, the center point of the overall curves for the country and the three regions gradually shifted to the right, indicating a steady increase in the allocation of pharmacists’ human resources.

Figure 5 Kernel density distribution of human resource allocation for hospital pharmacists. (A) national nuclear density, (B) Nuclear density in the eastern region, (C) Nuclear density in the central region, (D) Nuclear density in the western region.

Second, the main peak height of the distribution curve for the country and the three regions consistently increased, while the curve width narrowed, suggesting a trend of reduced disparity. However, from 2019 to 2022, the peak height of the central and western regions decreased, accompanied by a widening curve width, indicating increased disparities in these regions.

Third, the distribution curves for the national, eastern, and western regions exhibited “right-dragging” characteristics, expanding further to the right. This suggests that provinces with higher levels of human resource allocation strengthened their leading advantage.

Finally, the eastern region transitioned from a single-peak distribution to a distinct bimodal pattern, reflecting a pronounced gradient effect within the region. Conversely, the central and western regions transitioned from multi-peak distributions to single-peak patterns, indicating a weakening of multi-polarization within these regions.

Kernel Density Profile Analysis of PHC Pharmacists

Figure 6 illustrates the dynamic evolution of human resource allocation distribution for PHC pharmacists across China and its three major regions during the study period. First, the center points of the overall curves for China and the three regions gradually moved to the right, indicating a steady increase in the allocation of PHC pharmacists. Second, The height of the main peak of the distribution curve for China and the eastern region declined, while the curve width increased, suggesting a widening disparity. In contrast, the central and western regions showed an increase in peak height and a narrowing of width, indicating a reduction in disparities. However, between 2019 and 2022, disparities in the central and western regions widened slightly, likely influenced by the COVID-19 pandemic. Third, the distribution curves for China, the central, and the western regions exhibited “right-dragging” indicating that provinces with higher levels of human resource allocation experienced better development. The eastern region, however, displayed a coexistence of “right-dragging” and “left-dragging” phenomena, suggesting that the absolute disparity between advantaged and disadvantaged provinces within the region has been increasing. Finally, both the national and eastern regions transitioned from a single-peak to a bi-peak distribution, indicating an increase in multi-polarization. Meanwhile, the central and western regions did not show significant changes in their multi-peak patterns.

Figure 6 Kernel density distribution of human resource allocation for Primary Health Care pharmacists. (A) national nuclear density, (B) Nuclear density in the eastern region, (C) Nuclear density in the central region, (D) Nuclear density in the western region.

Spatial Convergence Analysis of Human Resources of Pharmacists in China

Spatial Correlation Test

Given the theoretical spatial correlation between pharmacist human resource allocation across provinces, spatial econometric analysis was conducted. Using Stata 17 software, the global Moran index for pharmacist human resource allocation in hospitals and PHCs in China was calculated from 2012 to 2022. As shown in Table 2, the global Moran index is positive in all years and statistically significant at the 5% level, indicating that pharmacist human resource allocation in China does not follow a completely random distribution. Instead, it exhibits a positive spatial autocorrelation, meaning that the development of pharmacist human resources in a given province is influenced not only by its own development level but also by neighboring provinces.

Table 2 The Spatial Correlation Test Results

From a temporal perspective, the Moran index for both hospitals and PHCs shows a downward trend, suggesting that the spatial clustering of pharmacist human resource allocation is gradually weakening. This decline may be attributed to the uneven regional development in China. As indicated by the kernel density measurements, pharmacist human resource allocation currently displays gradient disparities and tail-dragging phenomena, with provinces developing at varying speeds. These differences increase overall variability, thereby contributing to the weakening of spatial correlation over time.

Spatial Convergence Model Test

Convergence analysis examines whether there is a catch-up effect in pharmacist human resource allocation across provinces in China. If the allocation of pharmacist human resources across regions demonstrates convergence, it indicates that the gap between advantaged and disadvantaged provinces will continue to narrow; otherwise, the gap will widen. To explore this, the study employs a spatial β-convergence model to analyze the spatial convergence of pharmacist human resources across China and its three major regions from 2012 to 2022.

Before performing spatial β-convergence regression, the spatial fitness of the model must be tested. Following the research framework of J. P. Elhorst,^35, the paper uses three key tests: (1) the LM test is conducted based on OLS regression to evaluate spatial error effects and spatial lag effects within the model; (2) the Hausman test determines whether the spatial regression model fits fixed effects or random effects; and (3) the LR test and Wald test assess whether the model is applicable to the spatial Durbin model. As more models are involved in this paper, presenting the test results is more complicated, so only the test method is presented here, if there is a need for the test results can be asked to the corresponding author.

Absolute β Convergence Analysis

Considering the spatial correlation of pharmacist human resource allocation across the nation and within the three major regions—East, Central, and West—this study employs a spatial model to test for absolute β-convergence. As shown in Table 3, the β coefficients for pharmacist human resource allocation in hospitals and PHCs across the country and the three regions are significantly negative. This indicates a tendency toward absolute convergence, meaning that, after excluding external factors such as economic, social, and geographic influences, regions with lower levels of pharmacist human resource allocation exhibit higher growth efficiencies, allowing them to catch up with high-level regions and ultimately converge to a steady state.

Table 3 Absolute β Convergence Test of Pharmacist Human Resources Allocation

In terms of convergence speed, the national, eastern, and western regions demonstrate faster convergence in pharmacist human resource allocation in hospitals compared to PHCs. Conversely, the central region shows a significantly higher convergence speed in PHCs than in hospitals.

Condition β Convergence Analysis

As shown in Table 4, the β coefficients for the national level and all three regions are negative and statistically significant at the 10% level, suggesting the presence of conditional β convergence in pharmacists’ human resource allocation across hospitals and PHCs. This indicates that convergence trends persist even after controlling for economic, institutional, and demographic factors. Specifically, under absolute convergence, the convergence rates of pharmacists’ human resource allocation at both national and regional levels were significantly higher in hospitals than in PHCs. Regional rankings of convergence rates differed between healthcare tiers: at the hospital level, the northeast region exhibited the highest rate, followed by central, national, and western regions; whereas at the PHC level, the central region led, followed by national, eastern, and western regions. Notably, the western region consistently demonstrated the lowest convergence rates across both tiers.

Table 4 Conditional β Convergence Test of Pharmacist Human Resources Allocation

From a control variable perspective, GDP exhibited a significantly negative association only at the national hospital level and in the eastern PHC level, suggesting that economic development may paradoxically inhibit convergence. This phenomenon could be attributable to pharmacists’ preference for economically advantaged provinces, potentially widening interregional disparities through resource concentration in developed areas.³⁶ Health investment intensity showed a statistically significant positive association only at the national hospital level, indicating its effectiveness in promoting pharmacist resource convergence. Conversely, at the PHC level, a significantly negative association emerged in the central region, likely due to governmental prioritization of hospital-level infrastructure over PHC development.³⁷ Urbanization rate positively correlated with convergence at both national and eastern hospital levels, yet showed negative associations at the PHC level (except in the central region). This dichotomy may reflect urban healthcare systems’ capacity to attract PHC resources through superior facilities, particularly problematic in highly urbanized areas. The central region’s contrasting pattern could stem from its moderate urbanization allowing balanced resource distribution.³⁸ Government intervention demonstrated significantly negative effects on hospital-level convergence in both national and western regional analyses, potentially linked to policy emphasis on service scale expansion over workforce retention, driving pharmacist migration to private sectors.³⁹ Population density showed significantly negative associations with hospital-level convergence only at the national level and in the eastern region, likely exacerbated by rigid staffing quotas under China’s health ministry regulations (eg, mandatory 8% pharmacist-to-staff ratio in general hospitals). These institutional constraints, combined with lagged workforce adjustments amidst growing demand, may dilute per capita resources and accelerate staff turnover.⁴⁰

Discussion

Disparities Among Different Tiers of the Healthcare System Represent a Critical Determinant in the Allocation of Health Resources

Currently, there are more studies on the equity of China’s health system, such as Chai et al,⁴¹ Ao et al⁴² and Guo et al¹⁶ who made studies on the equity of China’s overall health human resources, rural health resources and pharmacist manpower, respectively, and concluded that the equity problems of China’s health resource allocation are widespread, and that the problems facing them are consistent across different healthcare systems. However, this study found that in the same health resources, there are large differences in the allocation of health resources at different levels, firstly, different levels of healthcare systems show different characteristics in regional differences, and secondly, the same factor does not have a consistent impact on different levels of healthcare systems. This is similar to Bin et al,²⁰ but they only explored the equity of pharmacist resource allocation in urban and rural area distribution from the perspective of spatial correlation, and concluded that the differences in pharmacist resource allocation mainly stemmed from the Chinese government’s negligence of geographic distribution equity. However, this study used the spatial β-convergence model to further investigate that different tiers of healthcare systems have significant differences in the speed of convergence, and that differences from different levels of healthcare systems are also important factors affecting pharmacist resource allocation in China. Kai et al⁴³ studied a single city, which also verified this phenomenon, and found that the human resources of pharmacists at the hospital level were significantly better than those at the primary care level. This suggests that differences in healthcare systems at different levels are widespread in China, and that distinguishing between the problems faced in different levels of healthcare systems is more conducive to narrowing the problem of regional differences in health resource allocation in China.

Marked Differences in the Allocation of Health Resources Between Medical Tiers and Regions

Despite substantial growth in China’s pharmacist workforce over the past decade, persistent structural disparities across healthcare tiers and regions remain evident. Pharmacist distribution remains heavily skewed toward hospitals, exhibiting an average disparity of approximately twofold compared to PHCs. This situation runs counter to the pyramid structure advocated for the rational allocation of healthcare resources, where primary healthcare institutions should hold a dominant position in resource allocation. Instead, the significant dominance of hospitals creates an imbalance in resource distribution, reducing the efficiency of healthcare services and making it difficult to meet the population’s demand for basic pharmacy services.^44, Regional heterogeneity in these disparities is pronounced: the economically advanced eastern region demonstrates relatively smaller hospital-PHC gaps, whereas the western region’s lower economic development correlates with wider disparities. The central region confronts compounded challenges including elevated population density coupled with moderate economic development, exacerbating its hierarchical resource disparities. China’s 13th Five-Year Plan (2016–2020) contributed to mitigating hierarchical disparities through its Basic Public Healthcare Services Equalization Initiative, which prioritized resource redistribution to PHCs.⁴⁰ While the COVID-19 pandemic accelerated short-term resource reallocation across tiers. This change is driven by the government’s substantial investments and temporary recruitment practices implemented in response to the COVID-19 pandemic. However, with the decline in local revenues and the conclusion of the pandemic, economically disadvantaged regions may find it challenging to sustain increased investments or may even face reductions in healthcare funding over the long term. This could further widen the gap between economically advantaged and disadvantaged localities. Therefore, efforts to increase healthcare staffing in the primary healthcare system should not focus solely on short-term investments but instead prioritize establishing a stable and sustainable growth mechanism.

Regional Differences with Significant Hierarchical Heterogeneity

Nationwide regional disparities in pharmacist distribution demonstrated a declining trajectory during the observation period (2012–2022), yet significant stratification-linked heterogeneity persists. Hospital-level disparities exhibited consistent contraction across all regions, aligning with Ni et al’s findings.¹⁵ Conversely, PHC-level disparities displayed divergent patterns: eastern PHC-level disparities displayed an upward trajectory, contrasting with narrowing gaps in central and western regions—a finding discordant with Bin et al’s western-focused equity analysis.²⁰ Kernel density estimation revealed strengthened primacy effects in eastern hospital clusters, while PHC-level analysis exposed widening absolute gaps between advantaged (eg, Shanghai, Zhejiang) and disadvantaged provinces. This siphoning phenomenon intensifies intra-regional pharmacist resource polarization, particularly affecting western rural areas.⁴⁵ Notably, the COVID-19 pandemic disproportionately impacted PHC resources, accelerating regional disparity expansion in central-western regions during 2019–2022.

Pharmacist Human Resources in PHC are Facing Double Siphoning

From the previous analysis, it is evident that pharmacist human resources in this PHC are simultaneously affected by the siphoning effect of dominant provinces and high-level medical institutions. This leads to a loss of PHC pharmacist human resources and exacerbates regional imbalances, a phenomenon amplified during the COVID-19 pandemic and other challenges. Despite increased investment and infrastructure development by the Chinese government in primary healthcare institutions, governmental interventions, based on control variables, exhibit an inhibiting effect on PHCs’ convergence. This is primarily because local government interventions in PHCs are often capital-focused, and economically developed regions possess stronger capabilities and efficiency in utilizing governmental funds, thereby amplifying their advantages in resource utilization. The efficiency advantage in economically developed regions further exacerbates the siphoning effect of advantaged provinces.⁴⁶ Moreover, the Chinese government’s focus on hospital-level healthcare infrastructure intensifies the siphoning effect on PHC pharmacist human resources, exacerbating PHC instability.⁴⁷ At the same time, the Chinese government has made various efforts to address this issue. For example, Chongqing Municipality has explicitly included a pharmacy workforce growth plan in its government work objectives. It has allocated special funds and established separate recruitment channels for pharmacists. Additionally, it has implemented a unified personnel management system, integrating pharmacists from both hospitals and PHCs while appointing them based on specific needs. These measures have proven effective in improving local pharmacist human resources and enhancing the pharmacy service capacity of PHCs.

Suggestions

Overall, China’s human resource allocation for pharmacists faces challenges such as an insufficient total number of pharmacists, increasing pressure on pharmacist services, and significant regional and hierarchical disparities. To address these issues, the Chinese government should formulate a medium- and long-term plan for the development of pharmacist human resources with reasonable growth targets. This plan should include encouraging the establishment of new pharmacy colleges and universities in the central and western regions and expanding enrollment in pharmaceutical programs to increase the overall supply of pharmacist human resources18. To address regional disparities in the distribution of pharmacists, the government should establish cross-regional mobility mechanisms. These mechanisms may include special subsidies and optimized talent introduction policies to facilitate the redistribution of surplus pharmacist resources from the eastern region to the central and western regions. Additionally, increased health funding and policy support should be provided to the central and western regions to create opportunities for new pharmacy graduates and attract talent.⁴⁸ Finally, to address hierarchical disparities, efforts should focus on supporting the development of primary health care institutions. This includes optimizing resource allocation systems, prioritizing the basic needs of primary health care institutions in terms of the number of pharmacists, hardware facilities, pharmacy service capacity, and providing special salary subsidies to pharmacists working at the grassroots level.

Limitations

Due to data limitations, this study focuses on the period from 2012 to 2022. While this timeframe effectively captures trends in the allocation of pharmacist human resources during this period, it may not fully reflect the current situation. Rapid changes in China’s healthcare policies, economic development, and demographic structure could have significantly altered the resource allocation landscape. For instance, certain regions may now be experiencing either a surplus or shortage of pharmacists. In particular, the eastern region may face an oversupply of pharmacists—a development not accounted for in this study.

Conclusion

This study found that the overall number of pharmacists in China continues to grow, but there are still significant differences in resource allocation. First, there are large differences in the allocation of different healthcare organizations, with the number of pharmacists in hospitals almost twice as many as in primary healthcare organizations, and this gap is particularly prominent in the central region. Secondly, there are significant regional differences, with the eastern region having a significant leading edge in resources and the central and western regions in a state of catching up, but the number of pharmacists in the west is still insufficient. Further analysis reveals that regional differences show different trends, with the regional gap in the distribution of hospital pharmacists gradually narrowing, but the gap in primary healthcare institutions intensifying, the internal imbalance in the east worsening, and the gap in the central and western regions widening due to the impact of epidemics. Spatial β-convergence indicated that there was a trend of convergence in pharmacist allocation across regions in China, but the rate of convergence was significantly higher in hospital tiers than in PHCs. These results suggest that pharmacist resource allocation is still facing both structural and regional challenges. This study goes beyond previous studies by revealing the dynamic changes in resource distribution and the impact of epidemics on regional disparities in terms of hierarchical and dynamic dimensions, emphasizing the lag and structural imbalance in resource adjustment in primary institutions. The findings have important implications for China’s health policy making, especially in terms of improving the equity of healthcare resources, optimizing the allocation of primary pharmacists and reducing regional disparities. Future workforce planning needs to further study the long-term trends in resource distribution and the effects of policy interventions on different tiers of the healthcare system, as well as strengthen support for central, western, and grassroots institutions, in order to promote a more balanced allocation of resources and to provide a guarantee for the achievement of the goal of universal healthcare.

Abbreviation

PHC, Primary healthcare institutions.

Data Sharing Statement

If someone wishes to access the original data, they should contact the corresponding author.

Funding

This work is supported by Hefei Maternal and Child Health Hospital hospital-level scientific research projects(yb2023_2_9), The second batch of key research projects on medical and health system reform by the Chinese Society of Health Economics, “Research on the Current Status of Coordination of the ‘Three Medical’ Policies and Guarantee Mechanisms” (CHEATGZZ20250203), Hefei Social Science Planning Project (HFSKYY202559), Research Project for Universities in Anhui Province (2023AH050717) Quality Project for Graduate Students in Anhui Province (RC2400000774) and Quality Project for Cultivating People in the New Era of the Education Department of Anhui Province (2023cxcys116). The funders had no involvement in study design, data collection, statistical analysis and manuscript writing.

Disclosure

The authors report no conflicts of interest in this work.

References

1. Organization WH. Global Strategy on Human Resources for Health: Workforce 2030. Geneva: World Health Organization; 2016.

2. Ling H, Huiqing M, Xiaojie S, et al. Construction of a data mining model for signaling serious adverse drug reactions based on real-world data–an empirical example of bilewood preparations. China Food Drug Admin Magazine. 2024;(10):102–111. doi:10.3969/j.issn.1673-5390.2024.10.012

3. Aitken M, Gorokhovich L. Advancing the responsible use of medicines: applying levers for change. SSRN Electron J. 2013.

4. Chen W, Yuqi F, Lulu S, Xi C, Hongmei L, Yingbo Z. Reflections on pharmacy professional and technical staffing in medical institutions in China. Chin Health Quality Manage. 2024;31(04):7–10+20. doi:10.13912/j.cnki.chqm.2024.31.4.02

5. Jiaqi W, Qiuchan Z, Yang ZG, Liu Y, Lili S. Analysis of the service capacity of primary healthcare institutions in Guangzhou. China Rural Health. 2023;15(10):16–19. doi:10.20126/j.cnki.1674-361X.2306-016

6. Park HR, Park SQ, Kim JH, Hwang JC, Lee GS, Chang JC. Geographic analysis of neurosurgery workforce in Korea. J Korean Neurosurg Soc. 2018;61(1):105–113. doi:10.3340/jkns.2017.0303.006

7. Silva M, de Almeida ATC, de Arauujo IT. Equity analysis of resource distribution for the popular pharmacy program. Revista de Saude Publica. 2019;5350. doi:10.11606/S1518-8787.2019053000731

8. Mollahaliloglu S, Yardim M, Telatar TG, Uner S. Change in the geographic distribution of human resources for health in Turkey, 2002–2016. Rural Remote Health. 2021;21(2):6478. doi:10.22605/RRH6478

9. Ismail M. Regional disparities in the distribution of Sudan’s health resources. East Mediterr Health J. 2020;26(9):1105–1114. doi:10.26719/emhj.20.056

10. Wang Z, He HY, Liu X, Wei HK, Feng QM, Wei B. Health resource allocation in Western China from 2014 to 2018. Arch Public Health. 2023;81(1):30. doi:10.1186/s13690-023-01046-x

11. Zhou MH. Equity and efficiency of health resource allocation in Sichuan Province, China. BMC Health Serv Res. 2024;24(1):1439. doi:10.1186/s12913-024-11946-5

12. Chen BQ, Jin FL. Spatial distribution, regional differences, and dynamic evolution of the medical and health services supply in China. Front Public Health. 2022;101020402. doi:10.3389/fpubh.2022.1020402

13. Liu T, Li JX, Chen J, Yang SL. Regional differences and influencing factors of allocation efficiency of rural public health resources in China. Healthcare. 2020;8(3):270. doi:10.3390/healthcare8030270

14. Gong YH, Ma D, Feng W. Study on the allocation efficiency of medical and health resources in Hainan Province: based on the super-efficiency SBM-Malmquist model. PLoS One. 2024;19(2):e0294774. doi:10.1371/journal.pone.0294774

15. Ni HF, Yuan M, Wan P, et al. The allocation and fairness of hospital pharmacist human resources in China: a time-series study. BMC Health Serv Res. 2024;24(1):1370. doi:10.1186/s12913-024-11882-4

16. Guo W, Zhang H, Lin YQ, et al. Analysis of relevant factors influencing the development of pharmacist human resources in China: based on the time series vector auto-regression model. Risk Manage Healthcare Policy. 2025;18:505–515. doi:10.2147/RMHP.S484621

17. Ning W, Guangpeng Z, Tiemin Z, Yuan D. Construction and optimization of a pharmacist demand measurement model for healthcare organizations. Chinese Health Econ. 2019;38(05):51–55. doi:10.7664/CHE20190513

18. Mingyue Z, Jiaming Z, Xuan L, Shijin Z. Research on the current situation and equity of pharmacist allocation in China’s healthcare institutions based on agglomeration degree. Chin Hospital Manage. 2024;44(06):56–60.

19. Ruochen Q, Qianggang B. Spatial distribution and dynamic changes of pharmacist resources in China. China Pharm. 2022;33(15):1898–1900+1910. doi:10.6039/j.issn.1001-0408.2022.15.20

20. Bin Z, Ying M, Ning Z, et al. Research on the allocation and optimization strategy of pharmacists and technicians in China based on spatial analysis. Chinese Health Service Manage. 2022;39(07):520–526.

21. YongZhi H. Study on the allocation status and fairness of hospital pharmacists in Tianjin after the new medical reform. China Pharm. 2014;25(13):1159–1161. doi:10.6039/j.issn.1001-0408.2014.13.03

22. Jie Z, Xiaona L. Measurement and decomposition of household income inequality of herders in Inner Mongolia-based on Dagum’s gini coefficient approach. J Inner Mongolia Agri Univ. 2024;26(06):1–10. doi:10.16853/j.issn.1009-4458.2024.06.001

23. Shen YQ, Sun ZS. Estimating the spatial correlation and convergence of China’s healthcare resources allocation: evidence from the Yangtze River Delta Region. Arch Public Health. 2022;80(1):207. doi:10.1186/s13690-022-00958-4

24. Zha QF, Liu Z, Song ZH, Wang J. A study on dynamic evolution, regional differences and convergence of high-quality economic development in urban agglomerations: a case study of three major urban agglomerations in the Yangtze river economic belt. Front Environ Sci. 2022;101012304. doi:10.3389/fenvs.2022.1012304

25. Dagum C. A new approach to the decomposition of the Gini income inequality ratio. Empirical Econ. 1997;22(4):515–531. doi:10.1007/BF01205777

26. Li S, Wenxiao F. Spatial differences and distributional dynamic evolution of science and technology financial efficiency in China. Manage Rev. 2021;33(01):44–53+67. doi:10.14120/j.cnki.cn11-5057/f.2021.01.003

27. de Campos AC, Lopes L, Carreira C. Spatial autocorrelation of exports and R&D expenditures in Portugal. J Knowledge Econ. 2024;15(2):8632–8653. doi:10.1007/s13132-023-01425-3

28. Sun J, Fan P, Wang K, Yu Z. Research on the impact of the industrial cluster effect on the profits of new energy enterprises in China: based on the Moran’s I index and the fixed-effect panel Stochastic frontier model. Sustainability. 2022;14(21). doi:10.3390/su142114499

29. Chen MH, Zhang TW, Chu QR, et al. Convergence analysis of inclusive green growth in China based on the spatial correlation network. Sustainability. 2023;15(16):12344. doi:10.3390/su151612344

30. Castellanos-Sosa FA, Cabral R, André Varella M. Energy reform and energy consumption convergence in Mexico: A spatial approach. Structural Change and Economic Dynamics. 2022;61:336–350. Available from: https://doi.org/10.1016/j.strueco.2022.03.007

31. Cavalieri M, Ferrante L. Convergence, decentralization and spatial effects: an analysis of Italian regional health outcomes. Health Policy. 2020;124(2):164–173. doi:10.1016/j.healthpol.2019.12.001

32. Tang H, Ao RJ, Li YL. The spatio-temporal pattern and its influencing factors of production efficiency of health resources in China. Front Public Health. 2024;121376518. doi:10.3389/fpubh.2024.1376518

33. Zhou Y, Zhao KX, Han JL, Zhao SD, Cao JY. Geographical pattern evolution of health resources in China: spatio-temporal dynamics and spatial mismatch. Trop Med Infect Dis. 2022;7(10):292. doi:10.3390/tropicalmed7100292

34. Yu YY, Lu JC, Dou XF, Yi YH, Zhou L. Spatial distribution and influencing factors of CDC health resources in China: a study based on panel data from 2016–2021. Front Public Health. 2024;121331522. doi:10.3389/fpubh.2024.1331522

35. Elhorst JP. Specification and estimation of spatial panel data models. Int Regional Sci Rev. 2003;26(3):244–268. doi:10.1177/0160017603253791

36. Xinghuan W, Lili C, Xinglong Z, Jun F. Evaluation of the fairness of health resource allocation in Guangxi’s specialized public health institutions based on HRAD-GIS method and its spatial difference analysis. Chin J Health Statis. 2024;41(05):742–746. doi:10.11783/j.issn.1002-3674.2024.05.023

37. Shizhong T, Mengmeng T, Xiaoyue L. Fiscal expenditure, government competition and health care service levels – an empirical analysis based on provincial panel entropy Tobit models. J Yunnan Univ Finance Econom. 2022;38(08):19–36. doi:10.16537/j.cnki.jynufe.000807

38. Yueran Y, Fei L. Research on the coupled and coordinated relationship between the degree of health technician agglomeration and the level of urbanization in China. Med Soc. 2024;37(06):59–65. doi:10.13723/j.yxysh.2024.06.009

39. Ruizhi P, Qiannan L. Expensive medical treatment and government expenditure in the context of healthy China: an analysis based on the perspective of structural imbalance of resource allocation. J Xi’an Jiaotong Univ. 2024;44(02):104–116. doi:10.15896/j.xjtuskxb.202402010

40. Yongyi G, Jing Z, Yunhan S, Yaru L, Xinran W, Xinyu L. Research on the equity of human resource allocation and its driving path in Chinese traditional medicine in China–an analysis based on fsQCA methodology. Chin J Health Policy. 2024;17(10):46–51. doi:10.3969/j.issn.1674.2982.2024.10.007

41. Chai YL, Xian GW, Kou RX, et al. Equity and trends in the allocation of health human resources in China from 2012 to 2021. Arch Public Health. 2024;82(1):175. doi:10.1186/s13690-024-01407-0

42. Ao YB, Feng QQ, Zhou ZL, Chen YF, Wang T. Resource allocation equity in the China’s rural three-tier healthcare system. Int J Environ Res Public Health. 2022;19(11):6589. doi:10.3390/ijerph19116589

43. Kai H, Weifang Y. Survey and analysis of human resources of Chinese pharmacists in Ningbo medical institutions. J Trad Chin Med Manag. 2025;33(1). doi:10.16690/j.cnki.1007-9203.2025.01.016

44. Zhiyong D, Chenxiao Z. Ten years of “New Medical Reform”: achievements, dilemmas and path choices of China’s medical and health care development. Reform. 2020;(09):149–159.

45. Di L, Yumei H. A study on the allocation and equity of health manpower resources in Beijing municipality. Hospital Manage Forum. 2024;41(11):15–18.

46. Weiwei W. New quality productivity, urban-rural income gap and shared prosperity. J Tech Econom Manage. 2025;(02):52–57.

47. Lina Y, Wei Z, Yaoyao J, Baiyue A, Yajuan L. Analysis of the current situation of talent aggregation and development strategies in public hospitals. Chin Hospital Manage. 2025;45(01):29–32.

48. Yu S, Shuang W, Li Q. A study of regional differences and dynamic evolution of the equity of general practitioner allocation in China. Health Econom Res. 2024;41(09):16–20. doi:10.14055/j.cnki.33-1056/f.2024.09.003

Introduction

Vitamin D₃ (VD₃), particularly its hormonally active form, calcitriol (1,25-dihydroxycholecalciferol), has immunomodulatory properties.¹ Its direct effects on nearly all immune cells, including lymphocytes, monocytes, and macrophages, as well as other cells within the tumor microenvironment (eg, fibroblasts and vascular endothelial cells), have been demonstrated.^1,2 Consequently, the impact of VD₃ on cancer development depends not only on its direct effects on cancer cells, extensively described in the scientific literature,^3–6 but also on its broader effects on the immune system and the tumor microenvironment.^7–10

While the immunosuppressive effects of VD₃ may offer potential benefits in cancer treatment;¹¹ some researchers suggest these same properties could also be detrimental.^12,13 Lymphocytes express the vitamin D receptor (VDR) upon activation, while dendritic cells and macrophages express it constitutively, making VD₃ a key modulator of immune and inflammatory responses.^1,12,14 Moreover, Th17 lymphocytes express VDR, and the proinflammatory cytokine IL-17A is modulated by VD₃ in both mouse and human T lymphocytes. Most studies suggest that calcitriol reduces Th17 cell recruitment and IL-17 secretion via the VDR-mediated pathway.^15–20 However, in our studies in young mice, calcitriol and its analogs (eg, tacalcitol) enhanced the lung metastatic potential of 4T1 mouse mammary gland cells,²¹ while screening revealed enhanced expression of genes associated with Th17 lymphocytes in the spleen: IL-17A (Il17a), RAR-related orphan receptor α (Rora), RAR-related orphan receptor γ (Rorc), IL-21 (Il21), IL-17 receptor E (Il17re), and IL-1 receptor type I (Il1r1).²² In contrast, applying the same treatment to aged, ovariectomized (OVX) mice bearing 4T1 tumors led to a temporary reduction in lung metastases²³ and decreased Rorc expression.²⁴

CD4⁺ splenocytes coming from young tacalcitol-treated mice and stimulated ex vivo to induce Th17 cells produced higher levels of IL-17A than those from untreated control mice, whereas the opposite effect was observed in aged OVX mice.²⁴ The effect of IL-17 varies depending on the stage of tumor development. In the context of chronic cancer and inflammation, IL-17’s tumor-promoting activity—mainly through the enhancement of angiogenesis—often surpasses its anticancer functions, including the stimulation of cytotoxic T lymphocytes and other immune cells that target tumors.^25,26

Current evidence indicate that osteopontin (OPN) is essential for dendritic cells to facilitate Th17 cell differentiation²⁷ and IL-17 production.²⁸ OPN regulates IL-17 expression through its receptors,²⁹ and the Spp1 (OPN) gene promoter region contains a vitamin D response element (VDRE).³⁰ The effects of VD₃ metabolites and analogs are mediated through VDR, which, upon ligand binding, dimerizes with the retinoid X receptor (RXR) and interacts with VDREs in the target genes, thereby influencing their transcription.³¹

We hypothesized that VD₃ regulates Th17 cell differentiation through VDR’s influence on OPN, potentially affecting tumor progression. Therefore, we aimed to analyze the role of OPN receptors in the differentiation of Th17 cells harvested from murine mammary gland tumors (4T1 and 67NR) treated with calcitriol and tacalcitol.

Materials and Methods

Cells

Mouse mammary gland cancer cell lines (4T1 and 67NR) were cultured under standard conditions to establish orthotopic tumor models. 4T1 and 67NR cells were procured from the American Type Culture Collection (ATCC, Rockville, MD, USA) and the Barbara Ann Karmanos Cancer Institute (Detroit, MI, USA), respectively. 4T1 cells were maintained in RPMI 1640 medium (Gibco™, Thermo Fisher Scientific, Waltham, MA, USA) with 10% (v/v) fetal bovine serum (FBS; HyClone, GE Healthcare, Chicago, IL, USA), 1 mM sodium pyruvate, and 3.5 g/L glucose (both Sigma-Aldrich, St. Louis, MO, USA). The 67NR cell line was cultured in Dulbecco’s modified Eagle medium (DMEM; Gibco, Scotland, UK) containing 10% (v/v) calf bovine serum (CBS; ATCC, Rockville, MD, USA), 1% amino acids, and 2 mM L-glutamine (both Sigma-Aldrich Chemie GmbH, Steinheim, Germany). Culture media were further supplemented with 100 µg/mL streptomycin and 100 U/mL penicillin (Sigma-Aldrich Chemie GmbH, Steinheim, Germany, and Polfa Tarchomin S.A., Warsaw, Poland, respectively). All cell lines were incubated at 37 °C in a humidified atmosphere of 5% (v/v) CO₂.

Mice

BALB/c/Foxp3GFP mice were generated for this study by systematically interbreeding BALB/c mice (Animalab, Poznań, Poland; Hirszfeld Institute of Immunology and Experimental Therapy, Polish Academy of Sciences [HIIET PAS], Wrocław, Poland) with C57BL/6/Foxp3GFP (C57BL/6-Tg(Foxp3-GFP)90Pkraj) mice, which express green fluorescent protein (GFP) under the Foxp3 gene promoter.³² Progeny were screened for Foxp3GFP reporter expression,³³ and all animals were housed at HIIET PAS, Wrocław, Poland.

Animal experiments were performed with approval from the Local Ethics Committee for Animal Experiments, Wrocław, Poland (permissions No. 50/2020 for experimental procedures and No. 44/2019, covering transgenic animal testing during breeding). All procedures adhered to the 3R principles, Directive 2010/63/EU, and national regulations. Mice were maintained in a specific pathogen-free facility with a 12/12-h light/dark cycle and provided SAFE 132 fodder (SAFE, Rosenberg, Germany) ad libitum.

This investigation utilized 6- to 8-week-old and 36- to 40-week-old mice. To establish the postmenopausal model, aged mice (30–35 weeks old) underwent ovariectomy (OVX) 5–10 weeks before the start of experiments. Ovariectomy was performed under general anesthesia using 3% (v/v) isoflurane inhalation (Aerane isofluranum, Baxter, Canada) in synthetic air (200 mL/min) and a buprenorphine injection (0.2 mg/kg; Orion Pharma Poland, Warsaw, Poland), as previously described.²³ As a control for the procedure, mice underwent sham surgery. The uterine weight confirmed successful ovariectomy (Supplementary Figure S1A). Post-surgical care included buprenorphine (0.1 mg/kg) for 24 h and as needed subsequently, with Dermatol powder (Galenic Laboratory, Olsztyn, Poland) applied to promote wound healing.²³

Scheme of Animal Studies

Orthotopic implantation of 4T1 (1×10⁴) or 67NR (2×10⁵) cells into the second right mammary fat pad was performed on day 0. Beginning 7 days post-tumor cell inoculation, mice received calcitriol (0.5 µg/kg) or tacalcitol (1 µg/kg) via oral gavage three times weekly. Tumor volume and body weight were monitored tri-weekly throughout the study. Tumor volume (TV) was determined using the following formula: , where a and b represent the shorter and longer tumor diameters, respectively, measured with a caliper. On day 14 (4T1) or day 18 (67NR), blood flow analysis was performed using contrast-enhanced ultrasonography (CEUS). Mice in the 4T1 cancer model were sacrificed on days 22–24 after cancer cell inoculation, while 67NR tumor-bearing mice were sacrificed on days 27–28. Tumors, spleens, lungs, livers, lymph nodes, brains, bone marrow, and blood were collected from all animals for further analysis (Figure 1).

Figure 1 Schematic overview of the experimental design and key methodological steps. Created in BioRender. Filip-Psurska, (B) (2025) https://BioRender.com/49ip4xx. Abbreviations: CEUS, contrast-enhanced ultrasonography; OPN, Spp1, osteopontin; VDR, vitamin D receptor; T-bet, Tbx21, T-Box Transcription Factor 21; RORa, RAR-related orphan receptor A; RORg t, RAR-related orphan receptor C (Rorc/ROR-γt); FoxP3, forkhead box P3; Gata3, GATA Binding Protein 3; Stat3, Stat5, signal transducer and activator of transcription 3; 5a, IL-17, Il17, interleukin 17a. Note: *Indicates aged mice.

For each tumor model, four independent experiments were conducted, with five mice per treatment group in each experiment. The exact number of animals or replicates used in each assay is stated in the figure legends.

Contrast-Enhanced Ultrasonography (Tumor Blood Perfusion)

Contrast-enhanced ultrasonography (CEUS) was used to assess tumor vascularization and blood flow dynamics. The MicroMarker™ contrast agent (VisualSonics, Ontario, Canada), approximately 10⁹ particles per mouse, was reconstituted in 1 mL of sterile 0.9% saline solution. Animals were anesthetized via continuous inhalation of 2–3% isoflurane (Baxter, Deerfield, Germany) in synthetic air (200 mL/min) and secured on the treatment platform. Following application of air bubble–free gel, the tumor’s central cross-section was visualized in the transverse plane with an MS-250S scanhead (Vevo 2100 ultrasound imaging system; VisualSonics). A 100 µL bolus of the contrast agent was administered intravenously, and the initial imaging sequence (bolus phase) was recorded at approximately 15 frames per second. Once the contrast signal within the tumor stabilized (ca. 50s), microbubbles in the field of view were disrupted using burst mode, and a subsequent imaging sequence (replenishment phase) was acquired. Post-imaging, mice were maintained in a warm environment until complete recovery. Data were analyzed using Vevo LAB 1.7.1 software with the VevoCQ modality (VisualSonics).

CD31 Immunohistochemistry (Blood Vessel Analysis)

Immunohistochemistry was performed to quantify tumor-associated blood vessels via CD31 staining. Detection of CD31 protein was performed using a primary anti-CD31 antibody (CD31 Polyclonal Antibody, series ZA4177363, ZC4226264, Thermo Fisher Scientific, Waltham, MA, USA) and a peroxidase-labeled secondary antibody (Goat anti-Rabbit IgG (H+L) Secondary Antibody, HRP, series YH375237, Thermo Fisher Scientific, Waltham, MA, USA). Microscopic evaluation of specimens was conducted to determine the number of blood vessels visible in the field of view. In manual assessment, a fully formed vessel or a group of positive cells was counted as one vessel, while single cells were excluded from the count. These analyses were performed by Sorbolab sp. z o.o., Poznań, Poland.

Blood Mononuclear Cell Isolation and Blood Morphological and Biochemical Parameters

Blood samples from experimental mice were collected into VACUTTE tubes containing lithium heparin (Greiner Bio-One, 454089; A-Biotech, Wrocław, Poland). Plasma was obtained by centrifuging the samples at 2000×g for 15 min at 4°C and subsequently stored at −80°C for further investigations.

For mononuclear cell isolation, blood cells were suspended in Hank’s Balanced Salt Solution (HIIET PAS, Wrocław, Poland) and subjected to gradient density centrifugation (Ficoll Paque Premium 1.084; Sigma-Aldrich, St. Louis, MO, USA) at 400×g for 40 min at room temperature. The collected mononuclear cells were washed with phosphate-buffered saline (PBS; HIIET PAS, Wrocław, Poland), centrifuged, and promptly utilized for cytometric staining.

Whole blood parameters were analyzed using a Mythic 18 hematology analyzer (Cormay, Warsaw, Poland). Plasma concentrations of calcium (Ca²⁺), creatinine (CRE2), aspartate aminotransferase (AST), and alanine aminotransferase (ALT) were determined using a Cobas c111 ISE biochemistry analyzer (Roche Diagnostics, Warsaw, Poland).

Cell Isolation for Further Analysis

Cells were isolated from various tissues by mechanical and enzymatic dissociation for downstream analyses.

Lymph nodes and spleens were mechanically dissociated by passing them through a 40 µm cell strainer into PBS (HIIET PAS, Wrocław, Poland) with 2% FBS (Sigma-Aldrich, St. Louis, MO, USA). The resultant cell suspension was centrifuged, and the supernatant discarded. For spleen samples, erythrocytes were lysed using Red Blood Cell Lysis Buffer (Sigma-Aldrich, St. Louis, MO, USA) at a 1:1 ratio for 1 min, followed by centrifugation and resuspension in PBS. The pellet containing separated cells was resuspended in PBS + 2% FBS and counted. The entire pellet of lymph node–isolated cells and a portion of spleen-derived cells were resuspended in 90% FBS + 10% DMSO; Tocris Bioscience, Bristol, UK) and frozen at −80°C. All steps were performed on ice using chilled reagents.

Tumor, brain, lung, and liver tissues were dissected using a scalpel and suspended in IMDM; Thermo Fisher Scientific, Waltham, MA, USA). DNase (Roche, Basel, Switzerland) and collagenase IA (Clostridium histolyticum collagenase; Sigma-Aldrich, St. Louis, MO, USA), both at a concentration of 1 mg/mL were added to the suspension and the mixture was incubated with shaking for 1 h at 37°C. Following enzymatic digestion, the suspension was filtered through a 40 µm cell strainer and centrifuged; the supernatant was subsequently removed. The cell pellet was resuspended in PBS containing 2% FBS and enumerated.

Remaining undigested tissue samples were placed in tubes with a ceramic homogenization bead (MP Biomedicals, Santa Ana, CA, USA) along with RIPA buffer containing protease and phosphatase inhibitor cocktails 2 and 3 (Sigma-Aldrich, St. Louis, MO, USA). Mechanical homogenization was performed using a FastPrep-24 Instrument homogenizer (5 m/s for 30s; MP Biomedicals, Santa Ana, CA, USA). Samples were incubated on ice for 30 minutes, flash-frozen in liquid nitrogen for 1 minute, and stored at −80°C. Except for the 37°C incubation, all procedures were conducted on ice using pre-chilled reagents.

Bone marrow: Femurs and tibiae were harvested from mice on ice in PBS containing antibiotics (100 U/mL penicillin, 100 µg/mL streptomycin; Polfa Tarchomin S.A., Warsaw, Poland, and Sigma-Aldrich Chemie GmbH, Steinheim, Germany, respectively). Bones were rinsed in PBS and cleaned of muscle tissue using a scalpel. After removing the epiphyses, the bone marrow was flushed out with PBS + 2% FBS. The collected suspension was centrifuged, and the supernatant was removed, followed by washing with PBS, centrifugation, and supernatant removal. The pellet containing separated cells was either resuspended in PBS + 2% FBS for immediate use or in 90% FBS + 10% DMSO and frozen in liquid nitrogen for storage.

Clonogenic Assay for Metastasis Identification

This assay was used to quantify metastatic cells within the lungs, liver, bone marrow, brain, spleen, and lymph nodes. Cells procured from these organs were propagated in medium appropriate for the specific tumor cell line initially implanted in the mice. Aliquots were seeded onto Ø 100 mm dishes with 10 mL of medium at the following densities: bone marrow (10×10⁶ cells); lymph nodes (0.5×10⁶ cells); lungs (0.25×10⁶ cells); liver (4.5×10⁶ cells); and spleen (5×10⁶ cells). For cultures derived from 4T1 tumor-bearing animals, 6-thioguanine (Sigma-Aldrich, St. Louis, MO, USA) was introduced 24 hours post-plating to a final concentration of 1 µM/mL. Media were refreshed one to two times weekly, contingent upon cellular proliferation rates. Incubation periods were tissue-specific: 3 weeks for preparations from lungs, spleen, liver, and lymph nodes, and 2 weeks for those from spleen and bone marrow. Material from identical tissue types, obtained from three murine groups matched for age and inoculated tumor cell line, was processed for equivalent durations.

After incubation, plates were rinsed with PBS, and colonies were stained using a 1% crystal violet solution in 80% methanol (Avantor, Gliwice, Poland) for 30 minutes at room temperature. Subsequently, plates were washed with distilled water and air-dried at room temperature for 24 hours. Images of the stained colonies were captured with the ChemiDoc Imaging System (BioRad, Hercules, CA, USA). The percentage of dish area covered by colonies was quantified using ImageJ 1.53q software (Wayne Rasband and contributors, National Institutes of Health, USA).

Magnetic Separation of CD4⁺ Splenocytes

CD4⁺ lymphocytes were isolated from mouse spleens using a magnetic separation kit (Miltenyi Biotec, Auburn, CA, USA). Splenocytes were centrifuged at 300×g for 7 min at 4°C and counted. Subsequently, the cells were resuspended in separation buffer (PBS with 0.5% bovine serum albumin (BSA) and 2 mM EDTA; pH 7.2; HIIET PAS, Wrocław, Poland) and incubated with an anti-CD16/CD32 blocking antibody (TruStain FcX, BioLegend, San Diego, CA, USA) for 5 min at 4°C to prevent non-specific Fc receptor binding. This step was omitted for cells designated for Th17 differentiation.

Magnetic separation proceeded according to the manufacturer’s instructions. Briefly, anti-CD4 magnetic beads (L3T4) MicroBeads (Miltenyi Biotec, Auburn, CA, USA) were added to the cell pellet, followed by a 30-minute incubation at 4°C. Post-incubation, cells were washed with separation buffer and centrifuged at 300×g for 7 minutes at 4°C. After aspirating the supernatant, the cell pellet was resuspended in fresh separation buffer. MS columns (Miltenyi Biotec, Auburn, CA, USA) were placed in a magnetic stand and equilibrated with separation buffer. The cell suspension was then applied to the column. After the flow-through was collected, the column was washed twice with separation buffer. To elute the magnetically retained CD4⁺ T cells, the column was removed from the magnetic field, and 1 mL of separation buffer was applied. The eluted cells were then counted and utilized for subsequent analyses.

Ex vivo Th17 Differentiation Assay

Th17 cells were generated ex vivo by differentiating CD4⁺ splenocytes isolated from young mice bearing 4T1 tumors. Initially, 24-well plates were coated with anti-CD3 and anti-CD28 antibodies (both 5 µg/mL; BioLegend, San Diego, CA, USA) in sterile PBS (500 µL/well) and overnight incubated at 4°C. Next, after washing plates twice with PBS, 0.5×10⁶ CD4⁺ splenocytes were seeded per well. The culture medium consisted of IMDM supplemented with GlutaMax and β-mercaptoethanol (5×10⁻⁵ M) (both Thermo Fisher Scientific, Waltham, MA, USA), 10% FBS (HyClone, GE Healthcare, Chicago, IL, USA), 100 µg/mL streptomycin, 100 U/mL penicillin, IL-6 (50 ng/mL), and TGF-β (1 ng/mL) (both BioLegend, San Diego, CA, USA). These cells were cultured for 4 days under conditions of 5% CO₂, 37°C, and 95% humidity.

After the 4-day differentiation, Th17 lymphocytes were collected for further analysis, and the corresponding supernatants were stored at −80°C for subsequent investigations. For flow cytometry, a subset of differentiated Th17 cells was additionally incubated with PMA (50 ng/mL), ionomycin (500 ng/mL), and brefeldin A (5 µg/mL) for 4 hours under 5% CO₂, 37°C, and 95% humidity.

Osteopontin Receptor Blocking During Th17 Differentiation

Blocking antibodies were used to study the role of OPN receptors (CD29, CD51, CD44) in Th17 cell differentiation. After CD4⁺ splenocyte separation from young mice bearing 4T1 tumors, cells were incubated with blocking antibodies targeting specific proteins (Table 1) in IMDM medium containing GlutaMax (Thermo Fisher Scientific, Waltham, MA, USA), 10% FBS; HyClone, GE Healthcare, Chicago, IL, USA), 100 U/mL penicillin, 100 µg/mL streptomycin, and β-mercaptoethanol (5 × 10⁻⁵ M, Thermo Fisher Scientific, Waltham, MA, USA). Incubation was performed for 30 min in uncoated 24-well plates under conditions of 5% CO₂, 37°C, and 95% humidity. After incubation, the cells and medium were collected, supplemented with differentiation factors (IL-6 at 50 ng/mL and TGF-β at 1 ng/mL; BioLegend, San Diego, CA, USA), and transferred to plates precoated with anti-CD3 and anti-CD28 antibodies (both 5 µg/mL, BioLegend, San Diego, CA, USA). The cells were then cultured for 4 days. At the final step of the experiment, PMA (50 ng/mL), ionomycin (500 ng/mL), and brefeldin A (5 µg/mL; all from Sigma-Aldrich, St. Louis, MO, USA) were introduced to each well for a final 4-hour stimulation. After completing the differentiation process, induced Th17 lymphocytes were analyzed using flow cytometry to assess intracellular IL-17 and IFNγ levels.

Table 1 Antibodies are Used to Block Osteopontin Receptors

Flow Cytometry (Extracellular and Intracellular Staining)

Flow cytometry was employed to quantify immune cell subsets and intracellular cytokine expression in isolated cells. Extracellular markers (CD3, CD4, CD25, CD29, CD44, CD51, CD61) and isotype controls were stained in specimens from each mouse using specific antibodies (Table 2). Cell viability was assessed with Fixable Viability Dye eFluor™ 780 (Invitrogen, Waltham, MA, USA), and non-specific antibody binding was prevented using TruStain FcX™ (anti-mouse CD16/32) (BioLegend, San Diego, CA, USA) to block Fc receptors.

Table 2 List of Antibodies and Fluorochromes Used in Flow Cytometry Analyses

Staining for extracellular markers: 1×10⁵ cells, initially suspended in PBS, were pelleted by centrifugation (350×g, 7 min, 4°C). The cell pellet was then suspended in 100 µL of Fixable Viability Dye eFluor™ 780 solutions in PBS + 2% FBS and incubated for 30 min in darkness at 4°C. After incubation, 700 µL PBS was added to each tube and centrifuged for 7 min, 350×g, at 4°C. Then, the cells were incubated in the dark for 5 min at 4°C with 50 µL of TruStain FcX™ blocking antibody. Subsequently, without an intervening centrifugation step, 50 µL of the appropriate extracellular staining antibody solution (in PBS + 2% FBS) was added, and samples were incubated for 30 minutes in darkness at 4°C. Following a final wash with 700 µL PBS and centrifugation, the cell pellet was resuspended in 200 µL of PBS + 2% FBS. Samples were analyzed using an LSR Fortessa flow cytometer (BD Biosciences, San Jose, CA, USA), with data compensation and analysis performed using FACSDiva software.

Intracellular staining: Cells derived from the same tissues were stimulated to produce cytokines. A total of 5×10⁵ cells were suspended in a stimulation medium (IMDM + GlutaMax + 10% FBS + brefeldin A (5 µg/mL) + PMA (50 ng/mL) + ionomycin (500 ng/mL)) and incubated for 4 h under 5% CO₂, 37°C, and 95% humidity. For staining, 1 × 10⁵ cells were resuspended in 100 µL of Fixable Viability Dye eFluor™ 780 solution prepared in PBS with 2% FBS and incubated in the dark for 30 min at 4°C. Cells were subsequently fixed by adding 500 µL of Fixation Buffer (BioLegend, San Diego, CA, USA) and incubating for 20 minutes in darkness at room temperature. Permeabilization was achieved by washing the cells three times with 500 µL of 1× Intracellular Staining Perm Wash Buffer (BioLegend, San Diego, CA, USA), centrifuging at 350×g for 7 minutes at room temperature after each wash. Next, 100 µL of TruStain FcX™ blocking antibody in 1× Intracellular Staining Perm Wash Buffer was applied for a 5-minute incubation in darkness at 4°C. The permeabilized cells were stained with antibodies against CD3, CD4, IL-17A, and IFNγ (Table 2), diluted in 100 µL of permeabilization buffer, and incubated at room temperature for 30 minutes in darkness. Post-incubation, cells were centrifuged, resuspended in PBS with 2% FBS, and analyzed on an LSR Fortessa flow cytometer (BD Biosciences, San Jose, CA, USA). FACSDiva software was used for data processing.

Quantitative PCR (qPCR) for Gene Expression

Gene expression in CD4⁺ lymphocytes was analyzed using quantitative real-time PCR (qPCR). RNA was isolated from CD4⁺ splenocyte lysates—previously stimulated with PMA and ionomycin and suspended in TRI-Reagent—using the Direct-zol™ RNA Miniprep Kit (ZYMO RESEARCH, Tustin, CA, USA), in accordance with the manufacturer’s instructions. DNase digestion was performed directly on the isolation columns during the procedure. RNA concentration was measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Reverse transcription to cDNA was performed with the SuperScript IV VILO Master Mix Kit (Thermo Fisher Scientific, Waltham, MA, USA) in a Veritii 9902 thermal cycler (Life Technologies, Carlsbad, CA, USA), employing the following thermal profile: 10 min at 25°C, 10 min at 50°C, and 5 min at 85°C. Subsequent gene expression analysis utilized customized TaqMan Array 96-Well FAST Plate Mouse FoxP3 plates, pre-coated with probes for specific target genes (Table 3).

Table 3 List of Probes Used in the Real-Time PCR Reaction

For gene expression analysis, 18.5 ng of cDNA (4T1 model) or 10 ng of cDNA (67NR model) per reaction was used. Each reaction mixture contained the specified amount of cDNA, 2× TaqMan™ Gene Expression Master Mix (Thermo Fisher Scientific, Waltham, MA, USA), and RNase-free DEPC-treated water. The reactions were carried out using the ViiA™ 7 device (Thermo Fisher Scientific, Waltham, MA, USA) with the cycling conditions: 40 cycles of 95°C for 15s and 60°C for 60s. Expression data were normalized using endogenous controls (B2m, Actb, and Sdha for 4T1; B2m and Actb for 67NR), and relative quantification (RQ) was determined by the ΔΔCt method, with samples from the control group serving as calibrators. Data analysis was performed using ExpressionSuite Software v1.3 (Thermo Fisher Scientific, Waltham, MA, USA).

Capillary Western Blot (Jess Simple Western)

Capillary Western blotting (Jess system) was used to quantify protein expression in CD4⁺ T cells. Protein concentration in isolated CD4⁺ splenocytes (after incubation with PMA and ionomycin) cell lysates was measured using the Bio-Rad DC Protein Assay Kit (Bio-Rad, Hercules, CA, USA). Samples were initially centrifuged (10,000×g, 10 min, 4°C), and the supernatant was transferred to a fresh tube. Concentration of lysates was performed using Amicon Ultra-0.5 mL Centrifugal Filter Unit columns (Sigma-Aldrich, St. Louis, MO, USA). Following a subsequent centrifugation (14,000×g, 15 min, 4°C), the resulting supernatant was collected. Protein levels were quantified against a standard curve generated with bovine serum albumin (BSA; 2–0.125 mg/mL; Bio-Rad, Hercules, CA, USA) using the Lowry method. Absorbance was measured on a Synergy H4 Hybrid Multi-Mode Microplate Reader (BioTek Instruments, Inc., Winooski, VT, USA).

Samples were prepared in 0.1× Sample Buffer (ProteinSimple, Bio-Techne, Minneapolis, MN, USA), and 4× Master Mix (EZ Standard Pack 1 or 3; ProteinSimple, Bio-Techne, Minneapolis, MN, USA) was added at a 4:1 ratio. The total protein amounts used for each assay were as follows: ERK/p-ERK (1.12 µg), p38/p-p38 (1.25 µg), and VDR (3.0 µg). Prepared samples were incubated at 95°C for 5 min using a heating block and then cooled on ice. Primary antibodies (Table 4) were prepared using a milk-free antibody diluent (ProteinSimple, Bio-Techne, Minneapolis, MN, USA) at the specified dilutions. For multiplex assays testing two proteins in the near-infrared (NIR) fluorescence and chemiluminescence (CHEMI) detection channels simultaneously, secondary antibodies were diluted 1:1000 (eg, AntiMouse NIR Detection Antibody/AntiRabbit NIR Detection Antibody and AntiMouse HRP Secondary Antibody/AntiRabbit HRP Secondary Antibody; ProteinSimple, Bio-Techne, Minneapolis, MN, USA).

Table 4 Antibodies Used for Immunodetection Using Jess Simple Western Method

To normalize the obtained results to the total protein level, the RePlex solution and Total Protein Detection solution (ProteinSimple, Bio-Techne, Minneapolis, MN, USA) were used following the manufacturer’s protocol. Additionally, luminol-S and peroxidase solutions (ProteinSimple, Bio-Techne, Minneapolis, MN, USA) were mixed at a 1:1 ratio. All solutions were applied to a 25-capillary plate (12–230 or 66–440 kDa separation modules; ProteinSimple, Bio-Techne, Minneapolis, MN, USA) designed for Jess Simple Western analysis. The analysis was performed using the Jess Simple Western System (ProteinSimple, Bio-Techne, Minneapolis, MN, USA). Data processing was conducted using Compass software for Simple Western analysis.

ELISA (Proteins and 25(OH)D₃ Quantification)

The levels of selected molecules were analyzed using the ELISA method in various samples, including plasma from 4T1- or 67NR-tumor-bearing mice to determine VEGF and 25(OH)D₃ levels, as well as supernatants from stimulated CD4⁺ splenocyte cultures and cultures of differentiated Th17 lymphocytes derived from the spleens of young, 4T1 tumor-bearing mice to determine OPN levels. Each ELISA assay was performed according to the respective manufacturer’s protocols (Table 5). Analyte concentrations in the samples were calculated from standard curves generated using CurveExpert 1.4 software.

Table 5 List of ELISA Kits Used in the Study

Statistical Analysis

All statistical analyses were conducted using GraphPad Prism version 7.1. Data distribution normality was evaluated using the Shapiro–Wilk test. Specific tests used: One-Way ANOVA with post hoc correction for normally distributed data, and Kruskal–Wallis test with correction for non-normal distributions. The specific statistical post hoc tests employed for individual data analyses are detailed in the corresponding figure legends. Inter-group differences were considered statistically significant when the p-value was less than 0.05.

Results

Metastasis of Mouse BC Cells to Lung, Liver, Spleen, Brain and Bone Marrow

The results of clonogenic assay have shown that tacalcitol increased metastases to the lung and liver in young 4T1 tumor-bearing mice. In aged mice, calcitriol reduced lung metastases, while both calcitriol and tacalcitol decreased liver metastases (Figure 2A–D). 4T1 cells were also detected in the brains of both young and aged mice, in the bone marrow of young mice, and in the spleens of aged mice, without significant effects of treatment (Supplementary Figure S2A–C). The growth of 67NR cells was observed in bone marrow cultures from both young and aged mice and in lymph node cultures from young mice. Treatment did not affect the growth of these cells (Supplementary Figure S2D and E).

Figure 2 Metastasis of 4T1 mouse mammary gland tumors to lungs and liver. Surface area of 4T1 colonies cultured from (A) lungs and (B) livers harvested from 4T1 tumor-bearing young and aged OVX mice. Representative images of stained colonies: (C) lung and (D) liver. Cells isolated from organs were cultured in a medium appropriate to the tumor cell line inoculated into the mice. The percentage of the dish area covered by colonies was analyzed using the ImageJ 1.53q program. N = 7–10. Statistical analysis: Dunn’s test except (B) aged OVX: Holm-Sidak’s multiple comparison test; *p < 0.05. Abbreviations: Ctrl, control; Cal, calcitriol; Tacal, tacalcitol; OVX, ovariectomized.

Characteristics of 4T1 and 67NR Tumor Progression Upon Calcitriol and Tacalcitol Treatment

To investigate the effects of calcitriol and tacalcitol on breast cancer progression, we assessed tumor growth, vascularization, and potential toxicity of treatment using caliper measurements, contrast-enhanced ultrasonography, CD31 immunohistochemistry, blood morphology, plasma biochemistry, and ELISA assay.

Treatment with calcitriol and tacalcitol did not affect 4T1 or 67NR tumor growth in either young or aged OVX mice. The body weight of mice decreases significantly only in 4T1 and transiently in 67NR tumor-bearing young mice treated with calcitriol (Supplementary Figure S1C–F).

To further assess the toxicity of VD₃ compounds, blood morphological and biochemical parameters were analyzed. 4T1 tumor growth induces leukocytosis,³⁴ a phenomenon confirmed in this study in both young and aged mice, with increased leukocytes, including lymphocytes, monocytes, and granulocytes, compared to healthy young or sham-operated aged mice. In young mice, tacalcitol further increased the number of these cells. 67NR tumors caused a significant but less pronounced increase in white blood cell count. The number of erythrocytes and platelets, as well as hemoglobin levels, remained unchanged (Supplementary Table S1). Calcitriol, but not tacalcitol, increased Ca²⁺ plasma level. The level of 25(OH)D₃ did not change significantly upon tumor growth and treatment. Calcitriol increased creatinine plasma levels in 4T1-bearing aged mice and 67NR-bearing young mice. Calcitriol and tacalcitol did not influence the plasma level of alanine transaminase (ALT) nor aspartate aminotransferase (AST). The level of ALT was decreased in tumor-bearing mice (Supplementary Table S2).

Time-intensity curve (TIC) parameters reflecting blood flow dynamics in tumor tissue were analyzed. These included: peak enhancement (PE), indicating the maximum intensity on the TIC and representing blood volume; time to peak (TTP), defined as the interval from baseline to maximum intensity; mean transit time (mTT), corresponding to the center of gravity of the best-fit function of echo power (or fitted signal); wash-in area under the TIC curve (WiAUC); wash-in rate (WiR), maximum slope, which reflects the rate of blood inflow; and the wash-in perfusion index (WiPI), calculated as WiAUC divided by rise time (RT, the time from the onset of enhancement to PE), serving as a measure of overall blood flow.

Treatment of young mice bearing 4T1 or 67NR tumors led to a similar pattern in TIC parameters. An increase in WiR and PE indicated enhanced blood inflow into tumor tissue (Figure 3A, B and Supplementary Figure S3). A lower TTP compared to controls suggested a longer time required to reach maximal inflow intensity (Figure 3A and B). Conversely, mTT was reduced in both treatment groups in 67NR tumor-bearing mice, whereas tacalcitol increased this parameter in 4T1 tumors, suggesting either a slower or faster reperfusion process (blood inflow from adjacent tissue), respectively (Supplementary Figure S3A and C). In aged OVX mice bearing 67NR tumors, most TIC parameters showed opposite effects compared to young mice (Figure 3C and Supplementary Figure S3D). In aged mice bearing 4T1 tumors, treatments did not significantly affect blood flow parameters (Supplementary Figure S3B).

Figure 3 Tumor angiogenesis parameters measured in young and aged OVX mice bearing 4T1 or 67NR mouse mammary gland cancer. (A–C) contrast-enhanced ultrasonography (CEUS) analysis of blood flow after contrast agent injection in (A) 4T1 bearing young mice, (B) 67NR bearing young mice, and (C) 67NR bearing aged OVX mice. N = 4–5 mice. Parameters describing blood flow in tumor tissue: PE—peak enhancement; WiR—wash-in-rate; TTP—time to peak, WiAUC—wash-in area under the curve. (D) Immunohistochemical examination of tumor CD31 expression in young mice bearing 4T1 and 67NR tumors. N = 4–5 mice. (E) Vascular endothelial growth factor (VEGF) plasma level in young mice bearing 4T1 and 67NR tumors; ELISA assay. N = 9–10 mice. (F) Representative pictures of WIR parameter measured in young and aged ovariectomized mice bearing 67NR tumors. (G) Representative images of CD31 tumor tissue staining. Scale bars = 50 µm. Statistical analysis: (A–C and E) Sidak’s, (D) Dunnett’s multiple comparison tests; *p < 0.05; **p < 0.01. Abbreviations: Ctrl, control; Cal, calcitriol; Tacal, tacalcitol; VEGF, vascular endothelial growth factor; OVX, ovariectomized.

Histopathological analysis of blood vessels (CD31 staining) showed an increased number of blood vessels in young mice bearing 67NR tumors (Figure 3D–G). In 4T1 tumors, both in young and aged mice, as well as in 67NR tumors growing in aged mice, treatments did not affect CD31 expression (Figure 3D–G and Supplementary Figure S3E). Plasma VEGF concentration was elevated in calcitriol- and tacalcitol-treated young mice bearing 67NR tumors (Figure 3E) but remained unchanged in aged mice and in both young and aged 4T1 tumor-bearing mice (Figure 3E and Supplementary Figure S3F).

IL-17 Expression in CD3⁺CD4⁺ Lymphocytes

CD4⁺ T cell subsets—including IL-17⁺ (Th17)—were characterized by flow cytometry in various tissues.

To characterize CD3⁺CD4⁺ lymphocytes, these cells were isolated from the lungs, tumor, and spleen of all mice, as well as from the blood and lymph nodes of 4T1 tumor-bearing mice. Supplementary Figure S4E–H summarizes the number of CD3⁺CD4⁺ lymphocytes in all treated mice. Significant differences in the percentage of these cells were observed in tumor tissue and lungs of aged OVX 4T1 tumor-bearing mice: tacalcitol increased CD3⁺CD4⁺ lymphocytes in both tissues, whereas calcitriol increased them only in the lungs (Supplementary Figure S4G).

Among CD3⁺CD4⁺ lymphocytes, we analyzed the population of IL-17-positive cells (Th17 lymphocytes). Tacalcitol increased Th17 lymphocytes in the lung tissue of young 4T1 tumor-bearing mice (Figure 4A). In aged OVX mice, the opposite trend was observed, with a significantly lower percentage of Th17 cells in tacalcitol-treated mice compared to those treated with calcitriol. Conversely, calcitriol increased the percentage of these cells in the tumor tissue of aged OVX 4T1 tumor-bearing mice (Figure 4B).

Figure 4 IL-17 expression (Th17 lymphocytes) in CD3+CD4+ lymphocytes isolated from 4T1 tumor-bearing mice treated with calcitriol and tacalcitol. Flow cytometry analysis of cells from (A) young and (B) aged OVX mice bearing 4T1 mouse mammary gland tumors. (C) Scheme of gating. N = 4–5 except spleens from aged OVX mice N = 10. Statistical analysis: (A) Dunn’s and (B) Dunnett’s multiple comparison tests; *p < 0.05, **p < 0.01. Abbreviations: Ctrl, control; Cal, calcitriol; Tacal, tacalcitol; OVX, ovariectomized.

Analysis of other tissues in 4T1 and 67NR tumor-bearing mice showed no differences among treatment groups (Figure 4 and Supplementary Figure S4A–D). The only exception was the blood cell analysis of young 4T1 tumor-bearing mice, where tacalcitol increased the percentage of IL-17-positive cells compared to calcitriol (Supplementary Figure S4A).

FoxP3 Expression in CD3⁺CD4⁺ Lymphocytes

CD3⁺CD4⁺ cells from tumors, blood, and selected organs were also analyzed for Foxp3 expression (CD3⁺CD4⁺Foxp3⁺) (Supplementary Figure S5). In the tumor tissue of aged OVX 4T1 tumor-bearing mice, tacalcitol treatment increased the percentage of FoxP3⁺ cells (Supplementary Figure S5A). In lung tissue, calcitriol (p = 0.0535) and tacalcitol (p = 0.0688) showed a trend toward increasing the Treg population in young 4T1 tumor-bearing mice (Supplementary Figure S5B). The percentage of FoxP3⁺ splenocytes decreased in aged 4T1 tumor-bearing mice treated with tacalcitol (Supplementary Figure S5C). In addition to tumors, regional lymph nodes and blood from 4T1 tumor-bearing mice were analyzed (Supplementary Figure S5D and E). The only significant change was observed in the blood of aged mice, where calcitriol decreased the percentage of FoxP3⁺ cells (Supplementary Figure S5E).

Characteristics of CD3⁺CD4⁺ Lymphocytes in Terms of Expression of OPN Receptors

CD3⁺CD4⁺ lymphocytes were analyzed for the expression of selected molecules described as OPN receptors, including three integrins (CD29, CD51, and CD61) and the CD44 receptor. In young mice bearing either 4T1 or 67NR tumors, the expression of these molecules was not significantly affected in various organs (Supplementary Figures S6A and S7A). The only observed change was an increase in CD44⁺ splenocytes in calcitriol-treated young 4T1 tumor-bearing mice (Supplementary Figure S6A). In aged 4T1 tumor-bearing mice, tacalcitol significantly decreased the percentage of CD29⁺ cells in tumor tissue, while calcitriol increased CD29⁺ cells in the blood (Supplementary Figure S6B). Tacalcitol increased the percentage of CD51⁺ cells in tumor tissue but decreased it in the spleen of aged 4T1 tumor-bearing mice. In the blood of these mice, calcitriol increased the percentage of CD44⁺ cells (Supplementary Figure S6B). In aged 67NR tumor-bearing mice, the only observed change was an increase in CD44⁺ cells in the spleen (Supplementary Figure S7B). Supplementary Figure S8 illustrates the gating strategy used for the OPN receptor expression analysis by flow cytometry.

Selected Genes and Protein Expression in CD3⁺CD4⁺ Splenocytes

Further characterization of CD3⁺CD4⁺ lymphocytes isolated from the spleen included the analysis of selected genes related to T lymphocyte subpopulation differentiation and the mechanisms of vitamin D₃ action (Figures 5, 6, Supplementary Figures S9 and S10). In 4T1 tumor-bearing mice, calcitriol and tacalcitol significantly increased Rorc expression, while tacalcitol also increased Tbx21 expression (Figure 5A and B). In aged mice, calcitriol increased Stat5a and Vdr expression (Figure 5C and E). Tacalcitol increased Spp1 mRNA expression in young 4T1 tumor-bearing mice (Figure 5F).

Figure 5 Selected genes and proteins expression in CD4+ T cells isolated from spleens of 4T1 tumor bearing young and aged OVX mice. (A) Rorc, (B) Tbx21, (C) Stat5a, (D) Il17a, (E) Vdr, (F) Spp1, (G) p-ERK/ERK ratio, (H) p-p38/p-38 ratio. (A–F) Analysis of gene expression was performed using specially designed TaqMan Array 96-Well FAST Plate Mouse FoxP3 plates coated with probes for the analysis of specific genes. Endogenous controls (B2m, Actb, and Sdha) were selected, and a comparative ΔΔCt analysis was performed. Samples from the group of control mice were selected as calibrators to obtain the RQ parameter. Analysis was performed using ExpressionSuite Software v1.3. (G and H) The Western blot analysis was performed in a Jess Simple Western System. The samples were normalized to total protein, and then the p-ERK/ERK or p-p38/p38 ratio was calculated. N = 4–6. Statistical analysis: (A–F) Sidak’s multiple comparisons test. (G and H) Tukey’s multiple comparisons test; *p < 0.05. Abbreviations: Il17a, interleukin 17A; Spp1, osteopontin; Rorc, RAR related orphan receptor C; Tbx21, T-box transcription factor 21; Vdr, vitamin D receptor; Stat5a, signal transducer and activator of transcription 5A; Ctrl, control; Cal, calcitriol; Tacal, tacalcitol; OVX, ovariectomized.

In aged OVX 67NR tumor-bearing mice, tacalcitol decreased Il17a expression (Figure 6A), while calcitriol reduced Foxp3 expression (Figure 6C). Increased Rora and Gata3 expression was observed in young 67NR tumor-bearing mice (Figure 6B and D). Spp1 expression was decreased by tacalcitol in young 67NR mice (Figure 6E). Analysis of phosphorylated ERK and p38 in CD3⁺CD4⁺ lymphocytes revealed an increase in ERK phosphorylation following calcitriol treatment in aged 4T1 tumor-bearing mice (Figure 5G–H). Additionally, tacalcitol increased VDR expression in young 67NR (Figure 6F) and aged 4T1 tumor-bearing mice (Supplementary Figure S11F). The expression levels of ERK, p-ERK, p38, and p-p38 are presented in Supplementary Figures S11 and S12.

Figure 6 Selected genes and proteins expression in CD4+ T cells isolated from spleens of 67NR tumor bearing young and aged OVX mice. (A) IL17a, (B) Rora, (C) Foxp3, (D) Gata3, (E) Spp1, (F) VDR. (A–E) Analysis of gene expression was performed using specially designed TaqMan Array 96-Well FAST Plate Mouse FoxP3 plates coated with probes for the analysis of specific genes. Endogenous controls (B2m and Actb) were selected and a comparative ΔΔCt analysis was performed. Samples from the group of control mice were selected as calibrators to obtain the RQ parameter. Analysis was performed using ExpressionSuite Software v1.3. (F) The Western blot analysis was performed in a Jess Simple Western System. The samples were normalized to total protein. N = 3–5. Statistical analysis: Sidak’s multiple comparisons test; *p < 0.05. Abbreviations: Il17a, interleukin 17A; Spp1, osteopontin; Rora, RAR related orphan receptor A; Gata3, GATA binding protein 3; Foxp3, forkhead box P3; VDR, vitamin D receptor; Ctrl, control; Cal, calcitriol; Tacal, tacalcitol.

Differentiation of Th17 Lymphocytes in the Presence of CD29, CD51 and CD44 Blocking Antibodies

To evaluate functional responses, CD4⁺ splenocytes were differentiated into Th17 cells ex vivo under polarizing conditions, with additional blockade of OPN receptors (CD29, CD51, CD44).

CD3⁺CD4⁺ splenocytes isolated from young 4T1 tumor-bearing mice were differentiated into Th17 lymphocytes. In tacalcitol-treated young mice, the percentage of lymphocytes expressing IL-17 was significantly higher (Figure 7A, control group); however, this effect was not observed when IFNγ⁺ cells were excluded from the analysis (Figure 7B). In mice treated with calcitriol and tacalcitol, an increased percentage of IFNγ⁺ cells among CD3⁺CD4⁺ splenocytes was noted (Figure 7C). When IL-17⁺ cells were excluded, a significant increase in IFNγ⁺ cells was observed only after calcitriol treatment (Figure 7D). Additionally, the percentage of IL-17⁺IFNγ⁺ double-positive cells increased in splenocytes from calcitriol-treated young mice during Th17 differentiation (Figure 7E).

Figure 7 Characteristic of Th17 lymphocytes ex vivo differentiated with CD29, CD51, and CD44 blockade. Percentage of (A) all IL-17+, (B) IL-17+IFNγ−, (C) all IFNγ+, (D) IL-17−IFNγ+, (E) double positive IL-17+IFNγ+ cells among CD3+CD4+ splenocytes differentiated toward Th17 in control conditions (left part of each graph) and with blocking antibodies (right part of each graph). Left column of graphs: CD29 blocking, middle column: CD51 blocking, and right column, CD44 blocking. CD3+CD4+ splenocytes isolated from 4T1 young mice treated with calcitriol and tacalcitol (blocked or not with antiCD29, antiCD51, and antiCD44 antibodies before seeding) were differentiated ex vivo on plates coated with antiCD3 and antiCD28 antibodies with the presence of IL-6 and TGF-β. On the last day, cells were incubated with PMA (50 ng/mL) + ionomycin (500 ng/mL) + brefeldin A (5 µg/mL) and incubated for 4 h. The cells were analyzed by flow cytometry for intracellular expression of IL-17 and IFNγ. N = 4–5. Statistical analysis: Sidak’s multiple comparisons test; *p<0.05. Abbreviations: IL-17, interleukin 17; IFNγ, interferon γ.

During culture in the presence of Th17 differentiation factors, blocking antibodies against CD29, CD51, and CD44 were used. Blocking CD29 did not affect IL-17 expression in control mice but inhibited tacalcitol-induced stimulation of IL-17⁺ cells (Figure 7A, left graph). Furthermore, when the IL-17⁺IFNγ⁻ cell population was analyzed under CD29 blockade in cells from both calcitriol- and tacalcitol-treated mice, a significant decrease in the percentage of these cells was observed (Figure 7B, left graph). Blocking CD29 also led to an increased percentage of IFNγ⁺ cells across all groups of mice (Figures 7C and 8C). However, the increased percentage of IFNγ⁺ lymphocytes as an effect of calcitriol and tacalcitol treatment was preserved despite CD29 blockade. IL-17⁻IFNγ⁺ cells with blocked CD29 generally behaved similarly to the entire IFNγ⁺ population, however, the effect of stimulating IFNγ expression by calcitriol and tacalcitol was more pronounced (Figure 7D). The pattern of IL-17⁺IFNγ⁺ double-positive cells after CD29 blockade was similar to that observed in other IFNγ-expressing cell subpopulations (Figure 7E).

Figure 8 Osteopontin level in culture media of CD4+ splenocytes and induces Th17 cells from young mice bearing 4T1 mouse breast cancer. (A) CD4+ splenocytes; (B) CD4+ splenocytes differentiated ex vivo in Th17 differentiation condition. N = 5. Statistical analysis: Sidak’s multiple comparison test; *p<0.05. Abbreviations: Ctrl, control; Cal, calcitriol; Tacal, tacalcitol.

Blocking CD51 did not significantly affect Th17 differentiation in young mice but led to an increase in IL-17⁺ splenocytes from calcitriol-treated mice compared to controls (Figure 7A, middle graph). Analysis of the IL-17⁺IFNγ⁻ cell subpopulation showed a decreased level of these cells in the CD51-blocked control group compared to non-blocked control cells. However, IL-17⁺IFNγ⁻ cells were stimulated following CD51 blocking when derived from both calcitriol- and tacalcitol-treated mice (Figure 7B, middle graph). Conversely, after CD51 blocking, IFNγ stimulation was inhibited in cells from calcitriol- and tacalcitol-treated mice (Figure 7C). This inhibition was confirmed for calcitriol treatment in IL-17⁻IFNγ⁺ cells. Under CD51-blocking conditions, tacalcitol stimulated IL-17⁺IFNγ⁻ lymphocytes (Figure 7D). The pattern of IL-17⁺IFNγ⁺ double-positive cells after CD51 blocking was similar to that of the IL-17⁻IFNγ⁺ subpopulation (Figure 7E).

Blocking CD44 did not significantly affect the percentage of IL-17⁺ control cells (Figure 7A, right graph). However, IL-17⁺IFNγ⁻ cells decreased in the CD44-blocked control group compared to non-blocked controls, while all IFNγ⁺ populations increased (Figure 7B). The stimulatory effect of calcitriol and tacalcitol under CD44-blocking conditions was observed for the IL-17⁺IFNγ⁻ cell population, similar to the IL-17⁺ cells (Figure 7A and B). In contrast to the overall IFNγ⁺ population, calcitriol decreased IL-17⁻IFNγ⁺ cells under CD44-blocking conditions (Figure 7D). The pattern of IL-17⁺IFNγ⁺ double-positive cells after CD44 blocking was similar to that of the IL-17⁻IFNγ⁺ subpopulation (Figure 7E). Supplementary Figure S13 presents the gating scheme for IL-17 and IFNγ intracellular analysis by flow cytometry.

OPN Level in Culture Supernatants

Supernatants from CD4⁺ splenocyte cultures and induced Th17 lymphocytes from young 4T1 tumor-bearing mice were analyzed for OPN secretion (Figure 8A and B). An increased concentration of OPN was observed in the culture of induced Th17 cells derived from tacalcitol-treated mice (Figure 8B).

Discussion

Our previous studies have shown that calcitriol and its analogs (tacalcitol, PRI-2205) exert differential effects on the metastatic progression of 4T1 breast cancer depending on the age of the mice.^21,23 We have attributed this to the immunomodulatory activity of vitamin D₃ derivatives.^22,35 It is well established that aging alters immune system function;³⁶ for example, dendritic cells—which are essential for Th17 cell maturation—may lose or diminish their activity with age.³⁷ Moreover, in our earlier work, we have demonstrated that vitamin D₃ compounds can modulate immune cell populations differently in young versus aged mice, both in healthy individuals³⁸ and in those bearing tumors.^7,24,35 The age-related hormonal changes (eg, ovariectomy in aged mice) also impact the outcome of VD₃-modulated immune responses.²⁴ These findings suggest that age-related changes in the immune microenvironment significantly influence the response to vitamin D₃-based treatments.

Previous studies demonstrated that calcitriol can stimulate blood flow in 4T1 mouse mammary gland tumors at a late stage of tumor development (days 21–24 after cell implantation).²¹ However, no changes in blood flow were observed in young mice bearing 67NR tumors at that time following calcitriol treatment.¹⁰ In the present study, we assessed blood flow in tumors at an early stage of development, revealing increased tumor blood flow and blood flow abnormalities in young mice bearing both metastatic 4T1 and nonmetastatic 67NR tumors. In addition, in young mice bearing 67NR tumors, these blood flow abnormalities were accompanied by a calcitriol-induced increase in vascular density, as well as an increase in pro-angiogenic VEGF plasma levels following both calcitriol and tacalcitol treatment. Interestingly, in aged mice, the effect on blood flow in 67NR tumors was the opposite, whereas no significant impact of the tested compounds was observed in 4T1 tumors, either at this early stage or in advanced tumors.²³

Previously, we also demonstrated that calcitriol, tacalcitol, or a cholecalciferol-rich diet enhanced 4T1 lung and bone marrow metastasis in young mice^7,21,39 while reducing this process in aged OVX mice.²³ Here, in addition to the lungs, we confirmed these previously reported findings by showing pro-metastatic effects in young mice and antimetastatic effects in aged OVX mice bearing 4T1 tumors in the liver.

We previously linked these changes in blood flow and metastasis to an increased expression of genes encoding transcription factors associated with Th17 differentiation in IL-6 and TGF-β-stimulated splenocytes from 4T1-bearing, tacalcitol-treated young mice. This effect was absent in aged mice.²⁴ Furthermore, calcitriol is known to regulate the expression of the OPN (Spp1) gene.⁴⁰ In our study, we observed a modulatory effect of calcitriol and its analogs on Spp1 expression. Specifically, in lung tissue from aged OVX mice bearing 4T1 tumors, Spp1 expression was decreased by calcitriol analogs on day 28 of the experiment,³⁵ whereas an increase in Spp1 expression was observed in young mice.²²

In the present study, we analyzed IL-17 expression in CD3⁺CD4⁺ T cells across different tissues in mice bearing both metastatic and nonmetastatic tumors. We also examined the expression of various molecules known as OPN receptors in these cells. Interestingly, tacalcitol increased the percentage of CD3⁺CD4⁺IL-17⁺ cells in the lungs of young 4T1 tumor-bearing mice while decreasing their percentage in aged OVX mice, aligning with the observed lung Spp1 expression.^22,35 IL-17, through NF-κB-mediated expression of MMP-2 and MMP-9, is a key driver of breast cancer invasiveness and metastasis.^41,42 Furthermore, OPN is known to promote Th17 lymphocytes and inflammation.^28,43

In experimental autoimmune encephalomyelitis (EAE), inhibition of Th17 cells following calcitriol treatment has been reported in young mice.¹⁹ Additionally, calcitriol has been shown to inhibit Th17 cell generation in vitro.⁴⁴ However, in our study, CD3⁺CD4⁺ splenocytes from young 4T1 tumor-bearing mice treated with tacalcitol (or both compounds) exhibited increased Rorc, Tbx21, and Spp1 mRNA expression. This may lead to the fact that they are more likely to differentiate to Th17 and produce IFNγ under Th17 differentiation conditions (Figure 9).

Figure 9 Calcitriol/tacalcitol stimulated tumor metastasis and blood flow in young mice bearing 4T1 murine breast cancer via impact on Th17 cells. The role of osteopontin receptors in this process is highlighted. Abbreviations: IL-17, interleukin 17; IFNγ, interferon γ; Spp1, osteopontin; Rorc, RAR-related orphan receptor C; Tbx21, T-box transcription factor 21; Vdr, vitamin D receptor; Stat5a, signal transducer and activator of transcription 5A. Note: *Published previously.21,22

The effects of calcitriol on CD3⁺CD4⁺ lymphocytes varied depending on mouse age and tumor type. Specifically, CD3⁺CD4⁺ splenocytes from aged OVX 4T1 tumor-bearing mice treated with calcitriol showed increased Stat5a and Vdr mRNA expression. Additionally, calcitriol stimulated p-ERK expression in CD3⁺CD4⁺ lymphocytes from these mice. Stat5a, a transcription factor essential for Th2 and Treg cell differentiation but not for Th17 or Th1 differentiation,⁴⁵ may explain the lack of an effect of calcitriol or tacalcitol on IL-17 expression in splenocytes. In young 67NR tumor-bearing mice, calcitriol increased Rora and Gata3 expression, while tacalcitol decreased Spp1. Conversely, in aged OVX 67NR tumor-bearing mice, tacalcitol reduced Il17a expression, and calcitriol decreased Foxp3. These opposing effects may contribute to the absence of significant differences in the CD3⁺CD4⁺ cell phenotype (IL-17 expression) in 67NR tumor-bearing young or aged mice.

IFNγ- and IL-17-co-expressing T cells have been documented in various human inflammatory autoimmune diseases.^46–48 A Th17-derived Th1 subset has also been identified within tumors, showing a positive correlation with the survival status of cancer patients.⁴⁹ Additionally, Muranski et al showed that Th17-polarized cells mediated the antitumor effect against advanced B16 melanoma in an IFNγ-dependent manner.⁵⁰ In our study, splenocytes differentiated toward Th17, harvested from young mice treated with calcitriol and/or tacalcitol, exhibited increased expression of both IL-17 and IFNγ and produced higher levels of OPN. Moreover, calcitriol increased also double positive T lymphocytes percentage.

In further studies aiming to explain the mechanism that may be responsible for the increasing effect of calcitriol and tacalcitol in young mice bearing 4T1 tumors on IL-17 and IFNγ expression during ex vivo Th17 differentiation, we used blocking antibodies targeting selected integrins and the CD44 molecule. Previous studies have shown that CD44 blockade, and to a lesser extent blockade of β1 integrin subunit (CD29), significantly reduced Th17 differentiation, whereas blocking CD51 (αv) had no such effect, indicating that OPN exerts its effects on Th17 differentiation primarily through CD44 and CD29.⁵¹ Other studies have highlighted the role of CD61 (β3) as an OPN receptor involved in Th17 differentiation.²⁹

CD29, CD51, and CD44 were selected for these studies because their expression on CD4⁺ cells from various organs was modulated to some extent by calcitriol and/or tacalcitol. Additionally, all three receptors have been implicated in Th17 differentiation.^43,51–53 Among the CD29 subunits, α4β1, and α5β1 are the predominant integrins expressed on T lymphocytes.^54,55 Furthermore, blocking CD29 or CD44 in CD4⁺CD45RA⁻CD25⁺ T cells stimulated with OPN, from acute coronary syndrome patients significantly decreased Th17 differentiation.⁴³

Our results suggest that calcitriol and tacalcitol may regulate Th17 differentiation in 4T1 tumor-bearing mice through CD29. Specifically, CD29 signaling appears to be responsible for stimulating Th17 (but not Th1) differentiation in response to tacalcitol. In contrast, CD51 and CD44 may counteract the pro-differentiation effects of calcitriol on Th17 cells. Additionally, CD51 signaling appears to be important for calcitriol- and/or tacalcitol-induced IFNγ (Th1) stimulation under Th17 differentiation conditions.

CD44 signaling appears to play a distinct role in the effects of calcitriol and tacalcitol on IFNγ⁺ cells. In young mice, blocking CD44 did not alter the stimulatory effect of calcitriol on all IFNγ-expressing cells. However, in IL-17⁻IFNγ⁺ cells, CD44 blockade reversed the effect of calcitriol from stimulation to inhibition. Tacalcitol, which stimulated only the total IFNγ⁺ cell population, lost this ability when CD44 was blocked, suggesting that CD44 plays a role in tacalcitol’s Th1 pro-differentiation activity.

Interestingly, in our previous studies, we have shown increased expression of CD44, Il23, and Irf4 (IL-23, IRF4, Th17 immune response drivers⁵⁶) in tumor-associated macrophages (TAMs) isolated from 4T1 tumors growing in young mice and treated with calcitriol, alongside elevated OPN levels in tumor tissue. Moreover, these TAMs exhibited M2 macrophage characteristics. Additionally, 4T1 cells expressed high levels of COX-2 and released PGE₂ following in vitro calcitriol treatment.⁷ Other studies using the 4T1 cancer model have shown that tumor-released PGE₂ induces IL-23 production in the tumor microenvironment, causing Th17 cell expansion.⁵⁷ Both IL-17 and PGE₂ may, in turn, stimulate M2 macrophage polarization.⁵⁸

These multidirectional effects of VD₃ in young 4T1 tumor-bearing mice make tumor-conducing microenvironment leading to increased metastatic potential.

Conclusion

Our study demonstrates that calcitriol and tacalcitol modulate Th17 cell differentiation and tumor progression in an age-dependent manner in murine models of breast cancer. These effects are linked to the expression and activity of OPN receptors on CD4⁺ T cells, particularly CD29, CD44, and CD51. We found that in young mice, CD29 contributed to the increased IL-17 expression in ex vivo differentiated Th17 cells derived from tacalcitol-treated mice, while CD51 and CD44 appeared to play opposing roles. Furthermore, IFNγ production in ex vivo differentiated Th17 cells from both calcitriol- and tacalcitol-treated young mice was mediated through CD51 integrin.

Taken together, our findings suggest that the balance between IL-17 and IFNγ in Th17 cells, shaped by vitamin D₃ signaling and OPN receptor engagement, may determine whether Th17 cells exert pro- or anti-tumor effects.

Ethics Statement

For the animal experiments, permissions No. 50/2020 and No. 44/2019 for testing transgenic animals during breeding were obtained from the Local Ethics Committee for Animal Experiments, Wroclaw, Poland. Committee Office: Institute of Immunology and Experimental Therapy, Polish Academy of Sciences, Weigla St.12, 53-114 Wrocław, Poland. All experiments were conducted by the 3R principles, Directive 2010/63/EU of the European Parliament and the Council, and national regulations.

Acknowledgments

Dedicated to Professor Hector F. DeLuca of the University of Wisconsin Madison on the occasion of his 95th birthday.

Funding

This work was supported by research grant no 2019/35/B/NZ5/01250 from the National Science Centre, Poland.

Disclosure

The authors report no conflicts of interest in this work.

References

1. Charoenngam N, Holick MF. Immunologic effects of vitamin D on human health and disease. Nutrients. 2020;12(7):1–28. doi:10.3390/NU12072097

2. El-Sharkawy A, Malki A. Vitamin D signaling in inflammation and cancer: molecular mechanisms and therapeutic implications. Molecules. 2020;25(14). doi:10.3390/MOLECULES25143219

3. Swami S, Krishnan A, Wang J, et al. Dietary vitamin D₃ and 1,25-dihydroxyvitamin D3 (calcitriol) exhibit equivalent anticancer activity in mouse xenograft models of breast and prostate cancer. Endocrinology. 2012;153(6):2576–2587. doi:10.1210/en.2011-1600

4. Swami S, Krishnan AV, Wang JY, et al. Inhibitory effects of calcitriol on the growth of MCF-7 breast cancer xenografts in nude mice: selective modulation of aromatase expression in vivo. Horm Cancer. 2011;2(3):190–202. doi:10.1007/s12672-011-0073-7

5. Marcinkowska E, Wallace GR, Brown G. The use of 1α,25-dihydroxyvitamin D3 as an anticancer agent. Int J Mol Sci. 2016;17(5). doi:10.3390/ijms17050729

6. Zhao C-N, Li Y, Meng X, et al. Insight into the roles of vitamins C and D against cancer: myth or truth? Cancer Lett. 2018;431:161–170. doi:10.1016/j.canlet.2018.05.039

7. Stachowicz-Suhs M, Łabędź N, Anisiewicz A, et al. Calcitriol promotes M2 polarization of tumor-associated macrophages in 4T1 mouse mammary gland cancer via the induction of proinflammatory cytokines. Sci Rep. 2024;14(1):3778. doi:10.1038/s41598-024-54433-x

8. Łabędź N, Anisiewicz A, Stachowicz-Suhs M, et al. Dual effect of vitamin D3 on breast cancer-associated fibroblasts. BMC Cancer. 2024;24(1):209. doi:10.1186/s12885-024-11961-z

9. Stachowicz-Suhs M, Łabędź N, Milczarek M, et al. Vitamin D3 reduces the expression of M1 and M2 macrophage markers in breast cancer patients. Sci Rep. 2024;14(1):22126. doi:10.1038/s41598-024-73152-x

10. Łabędź N, Stachowicz-Suhs M, Psurski M, et al. Modulation of fibroblast activity via vitamin D3 is dependent on tumor type—studies on mouse mammary gland cancer. Cancers. 2022;14(19):4585. doi:10.3390/cancers14194585

11. Liu W, Zhang L, Xu H-J, et al. The anti-inflammatory effects of vitamin D in tumorigenesis. Int J Mol Sci. 2018;19(9):2736. doi:10.3390/ijms19092736

12. Adorini L, Daniel K, Penna G. Vitamin D receptor agonists, cancer and the immune system: an intricate relationship. Curr Top Med Chem. 2006;6(12):1297–1301. doi:10.2174/156802606777864890

13. Cao Y, Du Y, Liu F, et al. Vitamin D aggravates breast cancer by inducing immunosuppression in the tumor bearing mouse. Immunotherapy. 2018;10(7):555–566. doi:10.2217/imt-2017-0131

14. Martens PJ, Gysemans C, Verstuyf A, Mathieu C. Vitamin D’s effect on immune function. Nutrients. 2020;12(5). doi:10.3390/NU12051248

15. Hayes CE, Hubler SL, Moore JR, Barta LE, Praska CE, Nashold FE. Vitamin D actions on CD4+ T cells in autoimmune disease. Front Immunol. 2015;6(March):1–22. doi:10.3389/fimmu.2015.00100

16. Huang P, He XY, Xu M. Correlation of vitamin D3 with the expression of RORγt and Foxp3 mRNAs in the peripheral blood of myasthenia gravis patients. Neuroimmunomodulation. 2020;27(2):97–103. doi:10.1159/000510861

17. Slominski AT, Kim TK, Li W, Yi AK, Postlethwaite A, Tuckey RC. The role of CYP11A1 in the production of vitamin D metabolites and their role in the regulation of epidermal functions. J Steroid Biochem Mol Biol. 2014;144:28–39. doi:10.1016/J.JSBMB.2013.10.012

18. Thangamani S, Kim M, Son Y, et al. Cutting edge: progesterone directly upregulates vitamin D receptor gene expression for efficient regulation of T cells by calcitriol. J Immunol. 2015;194(3):883–886. doi:10.4049/jimmunol.1401923

19. Joshi S, Pantalena L-C, Liu XK, et al. 1,25-dihydroxyvitamin D3 ameliorates Th17 autoimmunity via transcriptional modulation of interleukin-17A. Mol Cell Biol. 2011;31(17):3653–3669. doi:10.1128/MCB.05020-11

20. Slominski AT, Kim TK, Hobrath JV, et al. Endogenously produced nonclassical vitamin D hydroxy-metabolites act as “biased” agonists on VDR and inverse agonists on RORα and RORγ. J Steroid Biochem Mol Biol. 2017;173:42–56. doi:10.1016/j.jsbmb.2016.09.024

21. Anisiewicz A, Pawlik A, Filip-Psurska B, et al. Unfavorable effect of calcitriol and its low-calcemic analogs on metastasis of 4T1 mouse mammary gland cancer. Int J Oncol. 2017;52(1):103–126. doi:10.3892/ijo.2017.4185

22. Pawlik A, Anisiewicz A, Filip-Psurska B, et al. Calcitriol and its analogs establish the immunosuppressive microenvironment that drives metastasis in 4T1 mouse mammary gland cancer. Int J Mol Sci. 2018;19(7):2116. doi:10.3390/IJMS19072116

23. Anisiewicz A, Filip-Psurska B, Pawlik A, et al. Calcitriol analogues decrease lung metastasis but impair bone metabolism in aged ovariectomized mice bearing 4T1 mammary gland tumours. Aging Dis. 2019;10(5):977. doi:10.14336/AD.2018.0921

24. Pawlik A, Anisiewicz A, Filip-Psurska B, et al. Divergent effect of tacalcitol (PRI-2191) on Th17 cells in 4T1 tumor bearing young and old ovariectomized mice. Aging Dis. 2020;11(2):241. doi:10.14336/AD.2019.0618

25. Song Y, Yang JM. Role of interleukin (IL)-17 and T-helper (Th)17 cells in cancer. Biochem Biophys Res Commun. 2017;493(1):1–8. doi:10.1016/j.bbrc.2017.08.109

26. Gulubova M, Ananiev J, Ignatova M, Halacheva K. Pro-tumor and anti-tumor functions of IL-17 and of TH17 cells in tumor microenvironment. Acta Medica Bulg. 2016;43(2):68–79. doi:10.1515/amb-2016-0019

27. Shan M, Yuan X, Song LZ, et al. Cigarette smoke induction of osteopontin (SPP1) mediates TH17 inflammation in human and experimental emphysema. Sci Transl Med. 2012;4(117):117ra9. doi:10.1126/scitranslmed.3003041

28. Murugaiyan G, Mittal A, Weiner HL. Increased osteopontin expression in dendritic cells amplifies IL-17 production by CD4+ T cells in experimental autoimmune encephalomyelitis and in multiple sclerosis. J Immunol. 2008;181(11):7480–7488.

29. Zhao Q, Cheng W, Xi Y, et al. IFN-β regulates Th17 differentiation partly through the inhibition of osteopontin in experimental autoimmune encephalomyelitis. Mol Immunol. 2018;93:20–30. doi:10.1016/j.molimm.2017.11.002

30. Lau WL, Leaf EM, Hu MC, et al. Vitamin D receptor agonists increase klotho and osteopontin while decreasing aortic calcification in mice with chronic kidney disease fed a high phosphate diet. Kidney Int. 2012;82(12):1261–1270. doi:10.1038/ki.2012.322

31. Jeon S-M-M, Shin E-A-A. Exploring vitamin D metabolism and function in cancer. Exp Mol Med. 2018;50(4):20. doi:10.1038/s12276-018-0038-9

32. Kuczma M, Podolsky R, Garge N, et al. Foxp3-deficient regulatory T cells do not revert into conventional effector CD4+ T Cells but constitute a unique cell subset. J Immunol. 2009;183(6):3731–3741. doi:10.4049/jimmunol.0800601

33. Scirka B, Szurek E, Pietrzak M, et al. Anti-GITR antibody treatment increases TCR repertoire diversity of regulatory but not effector T cells engaged in the immune response against B16 melanoma. Arch Immunol Ther Exp. 2017;65(6):553–564. doi:10.1007/s00005-017-0479-1

34. DuPre’ SA, Hunter KW. Murine mammary carcinoma 4T1 induces a leukemoid reaction with splenomegaly: association with tumor-derived growth factors. Exp Mol Pathol. 2007;82(1):12–24. doi:10.1016/j.yexmp.2006.06.007

35. Anisiewicz A, Pawlik A, Filip-Psurska B, Wietrzyk J. Differential impact of calcitriol and its analogs on tumor stroma in young and aged ovariectomized mice bearing 4T1 mammary gland cancer. Int J Mol Sci. 2020;21(17):6359. doi:10.3390/ijms21176359

36. Salminen A, Kaarniranta K, Kauppinen A. Immunosenescence: the potential role of myeloid-derived suppressor cells (MDSC) in age-related immune deficiency. Cell Mol Life Sci. 2019;76(10):1901–1918. doi:10.1007/s00018-019-03048-x

37. Agrawal A, Gupta S. Impact of aging on dendritic cell functions in humans. Ageing Res Rev. 2011;10(3):336–345. doi:10.1016/J.ARR.2010.06.004

38. Śnieżewska A, Anisiewicz A, Gdesz-Birula K, Wietrzyk J, Filip-Psurska B. Age-dependent effect of calcitriol on mouse regulatory T and B lymphocytes. Nutr. 2024;16(1):49. doi:10.3390/NU16010049

39. Anisiewicz A, Kowalski K, Banach J, et al. Vitamin d metabolite profile in cholecalciferol-or calcitriol-supplemented healthy and mammary gland tumor-bearing mice. Nutrients. 2020;12(11):1–28. doi:10.3390/nu12113416

40. Shen Q, Christakos S. The vitamin D receptor, Runx2, and the notch signaling pathway cooperate in the transcriptional regulation of osteopontin. J Biol Chem. 2005;280(49):40589–40598. doi:10.1074/jbc.M504166200

41. Shibabaw T, Teferi B, Ayelign B. The role of Th-17 cells and IL-17 in the metastatic spread of breast cancer: as a means of prognosis and therapeutic target. Front Immunol. 2023;14. doi:10.3389/fimmu.2023.1094823

42. Roy L, Sahraei M, Schettini JL, Gruber HE, Besmer DM, Mukherjee P. Systemic neutralization of IL-17A significantly reduces breast cancer associated metastasis in arthritic mice by reducing CXCL12/SDF-1 expression in the metastatic niches. BMC Cancer. 2014;14(1):225. doi:10.1186/1471-2407-14-225

43. Zheng Y, Wang Z, Deng L, et al. Osteopontin promotes inflammation in patients with acute coronary syndrome through its activity on IL ‐17 producing cells. Eur J Immunol. 2012;42(10):2803–2814. doi:10.1002/eji.201242475

44. Ikeda U, Wakita D, Ohkuri T, et al. 1α,25-dihydroxyvitamin D3 and all-trans retinoic acid synergistically inhibit the differentiation and expansion of Th17 cells. Immunol Lett. 2010;134(1):7–16. doi:10.1016/j.imlet.2010.07.002

45. Lee J, Lozano-Ruiz B, Yang FM, Fan DD, Shen L, González-Navajas JM. The multifaceted role of Th1, Th9, and Th17 cells in immune checkpoint inhibition therapy. Front Immunol. 2021;12. doi:10.3389/fimmu.2021.625667

46. Kurschus FC, Croxford AL, P. Heinen A, Wörtge S, Ielo D, Waisman A. Genetic proof for the transient nature of the Th17 phenotype. Eur J Immunol. 2010;40(12):3336–3346. doi:10.1002/eji.201040755

47. Harbour SN, Maynard CL, Zindl CL, Schoeb TR, Weaver CT. Th17 cells give rise to Th1 cells that are required for the pathogenesis of colitis. Proc Natl Acad Sci U S A. 2015;112(22):7061–7066. doi:10.1073/PNAS.1415675112

48. Ronchi F, Basso C, Preite S, et al. Experimental priming of encephalitogenic Th1/Th17 cells requires pertussis toxin-driven IL-1β production by myeloid cells. Nat Commun. 2016;7. doi:10.1038/NCOMMS11541

49. Lei X, Xiao R, Chen Z, et al. Augmenting antitumor efficacy of Th17-derived Th1 cells through IFN-γ-induced type I interferon response network via IRF7. Proc Natl Acad Sci U S A. 2024;121(47):e2412120121. doi:10.1073/PNAS.2412120121/SUPPL_FILE/PNAS.2412120121.SAPP.PDF

50. Muranski P, Boni A, Antony PA, et al. Tumor-specific Th17-polarized cells eradicate large established melanoma. Blood. 2008;112(2):362–373. doi:10.1182/blood-2007-11-120998

51. Chen G, Zhang X, Li R, et al. Role of osteopontin in synovial Th17 differentiation in rheumatoid arthritis. Arthritis Rheum. 2010;62(10):2900–2908. doi:10.1002/art.27603

52. Guan H, Nagarkatti PS, Nagarkatti M. CD44 reciprocally regulates the differentiation of encephalitogenic Th1/Th17 and Th2/regulatory T Cells through epigenetic modulation involving DNA methylation of cytokine gene promoters, thereby controlling the development of experimental autoimmune encephalomyelitis. J Immunol. 2011;186(12):6955–6964. doi:10.4049/jimmunol.1004043

53. Su P, Chen S, Zheng YH, et al. Novel function of extracellular matrix protein 1 in suppressing Th17 cell development in experimental autoimmune encephalomyelitis. J Immunol. 2016;197(4):1054–1064. doi:10.4049/jimmunol.1502457

54. Iwata S, Ohashi Y, Kamiguchi K, Morimoto C. Beta 1-integrin-mediated cell signaling in T lymphocytes. J Dermatol Sci. 2000;23(2):75–86. doi:10.1016/S0923-1811(99)00096-1

55. Puente LG, Ostergaard HL. Beta 1/beta 3 integrin ligation is uncoupled from ERK1/ERK2 activation in cytotoxic T lymphocytes. J Leukoc Biol. 2003;73(3):391–398. doi:10.1189/JLB.0402199

56. Filip-Psurska B, Zachary H, Strzykalska A, Wietrzyk J. Vitamin D, Th17 lymphocytes, and breast cancer. Cancers. 2022;14(15):3649. doi:10.3390/cancers14153649

57. Qian X, Gu L, Ning H, et al. Increased Th17 cells in the tumor microenvironment is mediated by IL-23 via tumor-secreted prostaglandin E2. J Immunol. 2013;190(11):5894–5902. doi:10.4049/jimmunol.1203141

58. Liu L, Ge D, Ma L, et al. Interleukin-17 and prostaglandin E2 are involved in formation of an M2 macrophage-dominant microenvironment in lung cancer. J Thorac Oncol. 2012;7(7):1091–1100. doi:10.1097/JTO.0b013e3182542752

Introduction

Gynecologic malignancies may involve multiple pelvic and abdominal organs, necessitating extensive surgical resection.¹ Open surgery provides superior exposure of the operative field, facilitating comprehensive tumor removal and lymph node dissection to achieve curative intent. Moreover, laparoscopic approaches may carry a risk of tumor cell dissemination.² In open gynecologic oncology surgeries, pain mechanisms can be categorized into visceral and somatic pain. Visceral pain arises from direct tumor invasion, tumor-associated obstruction or perforation, tumor-induced inflammation or infection, and nerve traction or stimulation (activating splanchnic nerves). Somatic pain results from surgical trauma, inflammatory responses, and iatrogenic nerve injury. Collectively, these mechanisms contribute to postoperative pain. Epidural analgesia is widely regarded as a gold standard for postoperative pain management in major abdominal surgeries, including open gynecological malignancies, due to its superior blockade of nociceptive transmission.³ However, its clinical utility is hampered by several persistent limitations. The application of epidural analgesia is limited in duration, and after the removal of the analgesic pump, some patients may still experience severe pain.⁴ To address these challenges, adjuvant additives have been explored to synergize with local anesthetics.^5,6 Among these, N-methyl-D-aspartate (NMDA) receptor antagonists (e.g ketamine, esketamine) hold particular promise. Esketamine, the S(+)-enantiomer of ketamine, exhibits twice the analgesic potency of racemic ketamine and demonstrates prolonged pharmacokinetics (half-life: 15.4 hours), theoretically extending analgesia beyond epidural catheter removal.⁷

Recent advances in opioid-sparing strategies have underscored the potential of esketamine as a promising adjuvant agent. The study by Yan et al⁸ demonstrated that substituting epidural opioids with esketamine in thoracoscopic surgery reduced the incidence of chronic postoperative pain. However, the exclusive use of epidural esketamine, compared to opioids alone, has been associated with an increased incidence of short-term postoperative pain. Therefore, we aim to investigate whether the combination of opioids and esketamine can mitigate acute postoperative pain. Furthermore, although intravenous esketamine has demonstrated efficacy in reducing postoperative moderate-to-severe pain,⁹ its systemic administration may limit targeting of spinal NMDA receptors and increase neuropsychiatric side effects.¹⁰ Xu et al found that intravenous administration of esketamine can lead to neuropsychiatric adverse reactions, such as hallucinations, vivid dreams, and dizziness.¹¹ In contrast, a direct comparison in lower extremity surgery revealed that epidural esketamine (0.2 mg/mL) reduced VAS scores compared with IV administration, without increasing sedation or psychiatric adverse effects.¹² Additionally, epidural esketamine permits direct delivery to spinal dorsal horn NMDA receptors, potentially enhancing analgesia while minimizing systemic exposure.¹³ This hypothesis is supported by pharmacokinetic data demonstrating prolonged esketamine retention in cerebrospinal fluid following epidural administration (half-life: 15.4 hours vs 2.5 hours intravenously) [57]. As far as we know, there have been no clinical trials investigating the use of ketamine as an adjuvant in combination with opioids and local anesthetics for epidural analgesia.

To address the gap, we designed this prospective randomized controlled clinical trial, intending to verify our hypothesis that epidural esketamine could alleviate the rebound pain after removing the analgesic device.

Materials and Methods

Study Design

This is a prospective, randomized, double-blind, parallel-controlled clinical trial conducted at the Fudan University Shanghai Cancer Center in accordance with the CONSORT statement. The study protocol was approved by the Ethics Committee of the Fudan University Cancer Shanghai Center (Approval number: 2407300-10) and registered in the China Clinical Trials Registry (Registry number: ChiCTR2400091875). Written informed consent was obtained from each participant.

Inclusion Criteria

We recruited female participants aged 18–70 years who were scheduled for elective open surgery for gynecological malignancies. The inclusion criteria were as follows: American Society of Anesthesiologists (ASA) physical status I–III and body mass index (BMI) < 35 kg/m².

Exclusion Criteria

The exclusion criteria were as follows: (1) Inability to read or provide written informed consent. (2) Poorly controlled hypertension (ie, systolic blood pressure >180 mmHg or diastolic blood pressure >110 mmHg). (3) Coagulation dysfunction. (4) Contraindications to epidural puncture. (5) Intraocular hypertension or glaucoma. (6) Known or suspected glaucoma. (7) Severe coronary heart disease, history of myocardial infarction within the past 6 months, or chest pain following recent mild activity.

Randomization and Blinding

Enrolled patients were randomly allocated to either the eskeamine group or the control group in a 1:1 ratio. Randomization was performed by an independent researcher using SPSS software. The results of the group allocation were securely sealed in opaque envelopes. An anesthesiologist, not involved the anesthetic management, opened the envelope 1 hour before surgery. The independent anesthesiologist then prepared the analgesic medication and PCEA pump according to the allocated group. The anesthesiologists in charge of perioperative management remained blinded to the group allocation. After surgery, follow-up visits were conducted by a researcher who was also blinded to the group allocation.

Interventions

All of patients were monitored using electrocardiography, non-invasive blood pressure, and pulse oximetry when entering the operating room. Invasive arterial blood pressure and anesthesia depth with bispectral index (BIS) were monitored during the surgical procedures.

After central venous access was established, epidural puncture was performed at the T11-12 interspace, and a 5-cm catheter was inserted into epidural space. A test dose of 1% lidocaine (3 mL) was administered. Subsequently, the block level was assessed after the patient lies supine for 3 minutes. General anesthesia was induced intravenously with 2–3 mg/kg propofol, 0.3 µg/kg sufentanil, and 0.6–0.8 mg/kg rocuronium. And 5 mg dexamethasone and 3 mg granisetron were used before and after surgery, respectively, for postoperative nausea and vomiting (PONV) prevention.

Anesthesia depth was adjusted using target-controlled infusion of propofol, sevoflurane (0.5MAC) inhalation in an oxygen/air mixture, and 0.25% ropivacaine for epidural analgesia to maintain a BIS value between 40 and 60. The epidural administration protocol was as follows: An initial dose of 7–8 mL was administered in divided boluses before skin incision, with intervals of 3–5 minutes. Subsequently, 3–5 mL was administered every 1 hour. At the end of the procedure, 50 mg flurbiprofen was administered for multimodal analgesia. After surgery, sugammadex was administered at a dose of 4 mg/kg to antagonize muscle relaxation.

In postanesthesia care unit (PACU), all patients received a PCEA pump for pain control. For the esketamine group, the analgesic pump contained esketamine 0.4 mg/kg, sufentanil 100 µg, and ropivacaine 300 mg. For the control group, the analgesic regimen consisted of sufentanil 100 µg and ropivacaine 300 mg. The total capacity of the patient-controlled analgesia pump was 200 mL (sufentanil 0.5 µg/mL, ropivacaine 1.5 mg/mL), with a background infusion rate of 4.0 mL/h, a bolus dose of 4.0 mL, and a lockout interval of 15 minutes. Analgesia was maintained via the PCEA pump for 48 hours postoperatively. If a patient pressed the PCEA button during the lockout interval, the request was recorded as an additional press, but no additional analgesics were administered, which was defined as a failed press. Postoperative pain was evaluated using the Numeric Rating Scale (NRS) (NRS; 0 = no pain, 10 = worst imaginable pain).

A PACU nurse, who was not involved in anesthesia or postoperative follow-up, assessed the patient’s NRS, checked for the presence of dizziness, delirium, nausea and vomiting, and recorded the patient’s awakening time and medications administered in the PACU. Patients were discharged from the PACU once their Steward scores were greater than 6.

Outcomes

The primary outcome was the pain score on postoperative day 3 (POD3), assessed using the NRS. This included pain at rest, pain during activity, nocturnal pain, and the most severe pain experienced after removal of the analgesia pump.

Secondary outcomes included pain scores on postoperative days 1 (POD1), 2 (POD2), and 4 (POD4), as well as the Quality of Recovery-15 (QoR-15) score at POD3; the duration of surgery; pain scores and medication use in the PACU; the number of presses and failed presses on the PCEA pump; sleep duration and quality on POD3 and POD4; and complications related to ketamine (dizziness, nausea, vomiting, delirium, sleepiness, dreaming, nightmares). On POD1-4, an anesthesiologist, who was not involved in the management of anesthesia and was unaware of the group assignments, assessed the patient’s postoperative pain. All follow-up evaluations were uniformly scheduled between 10:00 and 11:00 a.m. On postoperative day 2 (POD2), at approximately 11:00 a.m., the epidural analgesia pump was removed by a designated anesthesiologist, who was also responsible for conducting subsequent follow-up assessments. Additionally, the Quality of Recovery-15 (QoR-15) score was assessed on the POD3.

Sample Size Calculation

The sample size was calculated using PASS software (version 2021). Based on previous clinical trials, the mean score of breakthrough pain within 24 hours after removal of the analgesic pump was approximately 5 points. Since pump removal occurred at approximately POD2 11:00 a.m., this 24h period corresponds to POD3 assessments. Our primary outcome was the pain score on Postoperative Day 3 (POD 3). This specifically captured the peak intensity of breakthrough pain occurring between pump discontinuation on POD 2 and POD 3. We assumed that epidural esketamine would reduce the pain score to 4 points, with a standard deviation of 1.5. With a power of 80% and an α level of 0.05, 37 patients per group were required to detect differences. To account for a potential 20% dropout rate, 44 patients per group were finally enrolled in the trial.

Statistical Analysis

Data normality was assessed using the Kolmogorov–Smirnov test. Continuous data were presented as mean ± standard deviation or median (interquartile range), as appropriate. Continuous outcomes were compared between groups using the Mann–Whitney U-test or independent samples t-tests, depending on their distribution. The non-normal distribution data and nocturnal pain scores were analyzed using the Wilcoxon signed-rank test. Categorical data were analyzed using the chi-square test. The significance level was set at 0.05 for all analyses. Statistical analyses were performed using SPSS software (Version 26.0, USA).

Results

A total of 88 patients undergoing open gynecological oncological surgery were assessed for eligibility. One patient, who was transferred to the ICU, was excluded from the analysis. All patients completed the follow-up. The study flowchart is presented in Figure 1.

Figure 1 The flow chart of the study. Abbreviations: ICU, Intensive Care Unit. POD, postoperative day; NRS, Numeric Rating Scale; QoR-15, Quality of Recovery-15.

The patient demographics were comparable between the two groups (Table 1), except that the prevalence of hypertension was higher in the esketamine group than in the control group (22.7% vs 7.0%; P = 0.04). Subgroup analysis by non-hypertension status (Table 2) revealed consistent reductions in POD2 nocturnal pain with esketamine in non-hypertensive patients ((1, 3) vs (1, 4.5); P=0.023).

Table 1 Demographic and Clinical Characteristics at Baseline

Table 2 Principal Findings in Non-Hypertensive Patients

No significant differences in pain intensity scores were found between the two groups at POD3 (Figure 2). However, significant differences in pain intensity at POD2 night were observed between the esketamine group and control group. The pain scores are presented in Table 3. In the esketamine group, less patients required rescue analgesics at POD3 (13.4% vs 32.6%; P = 0.036). Regarding adverse events, the incidence of postoperative dizziness was significantly higher in the esketamine group than in the control group (47.7% vs 23.3%; P = 0.017). However, no significant difference in dizziness was observed between the two groups on POD3. Additionally, no significant differences were observed in bolus analgesia times, incidences of PONV and nightmares. The median QoR-15 scores did not differ between the esketamine group and control group on POD3. Nevertheless, a subgroup analysis on QoR-15 indicated that moderate pain was less duration in the esketamine group than in the control group (median, 10 vs 9; P = 0.048). No significant differences were observed in postoperative severe pain, sleep quality, anxiety, and depression between the two groups.

Table 3 Comparisons of Postoperative Outcomes Between the Two Groups

Figure 2 Violin plots of nocturnal pain intensity in two groups at postoperative 2 and 3 days. Abbreviations: EG, esketamine group; CG, control group. Notes: The nocturnal pain at POD2 was lower in the esketamine group compared with that in the control group (P=0.044). Within the violin plot, the solid line represents the median, and the dotted line represents the quartile.

No significant differences in other pain intensity scores during the postoperative first 4 days, including rest pain, movement pain, and breakthrough pain, which is shown in Figure 3.

Figure 3 Violin plots of rest pain, movement pain and breakthrough pain intensity in two groups at multiple time points. Abbreviations: EG, esketamine group, CG, control group. Notes: (A) Rest pain intensity. (B) Movement pain intensity. (C) Breakthrough pain intensity. Within the violin plot, the solid line represents the median, and the dotted line represents the quartile.

Discussion

Our study demonstrated that the addition of esketamine to PCEA solution significantly reduced the incidence of nocturnal pain on POD2 and alleviated moderate pain duration on POD3 in patients undergoing open gynecological malignancy surgery. However, no significant differences were observed in rest pain, movement pain, or breakthrough pain during the POD1-4 days. Although hypertension prevalence differed at baseline, subgroup analysis analyses demonstrated that this imbalance did not confound our primary findings. These findings suggest that the analgesic effect of epidural esketamine is time-dependent and may be related to its pharmacokinetic profile and mechanism of action.

The observed reduction in nocturnal pain at POD2 may be attributed to the dual action of esketamine as an NMDA receptor antagonist and a modulator of central sensitization. Surgical trauma induces peripheral nociceptor activation and spinal dorsal horn synapses glutamate release, which promotes NMDA receptor-mediated hyperalgesia and wind-up phenomena.¹⁴ By blocking NMDA receptors, esketamine likely attenuates central sensitization, thereby mitigating delayed pain amplification following cessation of epidural analgesia.¹⁵ This is consistent with previous studies showing that NMDA antagonists reduce rebound pain following withdrawal of regional anesthesia.^16,17 Actually, intravenous esketamine could reduce postoperative pain scores, decrease the incidence of rebound pain following cessation of thoracic paravertebral block, and reduces opioid consumption.¹⁸ Epidural esketamine also can lower opioid consumption, provide better postoperative analgesia, reduce the rescue analgesics.⁷ Moreover, the extended half-life of esketamine (15.4 hours) may prolong its analgesic effects,⁷ potentially explaining the reduction in nocturnal pain at POD2 rather than at earlier timepoints. However, the transient nature of this effect underscores the necessity for optimized dosing regimens to sustain therapeutic benefits.

The modest improvement in moderate pain at POD3 by evaluating the QoR-15, coupled with the lack of differences in other pain metrics, raises concerns regarding the adequacy of the esketamine dose (0.4 mg/kg/48h) used in this trial. In the context of pain assessment utilizing the QoR-15 scale, a lower score is associated with prolonged pain duration and inferior recovery quality. Previous studies employing higher doses esketamine doses (1–2 mg/kg/48h) have reported more pronounced reductions in postoperative pain scores and opioid consumption.^19–21 Esketamine has been utilized both intravenously and epidurally. For instance, Zhang et al demonstrated that the addition of 1 mg/kg esketamine to PCIA reduced moderate-to-severe pain by 40% in patients undergoing thoracoscopic surgery.⁹ A recent study reported that 25 mg esketamine was added to PCEA.⁸ Our lower epidural esketamine dose may have resulted in subtherapeutic plasma concentrations, which were insufficient to fully suppress NMDA receptor activity beyond the immediate postoperative period. The dissociation between significant POD2 nocturnal pain reduction and non-significant POD3 pain scores may reflect harmacokinetic-pharmacodynamic Mismatch; Epidural esketamine’s CSF half-life (15.4 hr) predicts subtherapeutic concentrations at POD3 assessment. Pharmacometric modeling indicates: POD2 night (36 hr): CSF (esketamine) ≈ 55 ng/mL (> IC₅₀ 40 ng/mL),POD3 (48 hr): CSF (esketamine) ≈ 25 ng/mL (< IC₅₀).This aligns with animal data showing NMDA blockade efficacy requires >35 ng/mL CSF levels.⁷ On the other hand, differential modulation of pain phenotypes; Early nocturnal pain on POD 2 primarily involves surgically induced incisional pain and visceral traction pain, while later pain predominantly arises from central sensitization. By POD3, pain transitions to neuroplasticity phases, in which cytokine pathways dominate.²² Additionally, heterogeneous pain trajectories post-pump removal may account for these findings. It exhibits individualized differences in pain experience and therapeutic response. Future trials should explore dose escalation strategies while balancing potential side effects, such as dizziness, which was more prevalent in our esketamine group.

The selective reduction in nocturnal pain at POD2 has significant clinical implications. Postoperative sleep disruption, which is often exacerbated by nighttime pain, is associated with delayed recovery and increased morbidity.²³ By targeting this critical period, esketamine may enhance sleep quality and patient satisfaction,²⁴ despite the absence of differences in overall pain scores. This is consistent with our subgroup analysis, which showed improved QoR-15 scores for moderate pain in the esketamine group. However, the lack of significant differences in overall QoR-15 scores suggests that the benefits of the current regimen may be confined to specific pain dimensions rather than overall recovery metrics. NRS reduction in POD2 nocturnal pain translates to decrease in sleep-interrupting pain episodes (NRS>4). This aligns with ERAS benchmarks where NRS≤3 enables uninterrupted >4hr sleep blocks – faster functional recovery in patients.²⁵ Clinicians might consider time-limited esketamine infusions tailored to the peak rebound pain period (eg, 24–48 hours after pump removal) to maximize cost-effectiveness and minimize side effects.

Although epidural esketamine has shown promise, its role within multimodal analgesia framework warrants further discussion. Recent evidence supports the combination of NMDA antagonists with other adjuvants (e.g, gabapentinoids or COX-2 inhibitors) to synergistically target both peripheral and central pain pathways.^26,27 For instance, a randomized trial by Zhang et al demonstrated that the combination of intravenous esketamine with dexmedetomidine reduced improved analgesia and sleep quality after scoliosis correction surgery.²⁸ Another randomized trial by Yan et al found that combining epidural esketamine with intravenous parecoxib reduced both acute and chronic post-thoracotomy pain compared to either drug alone.⁸ In our study, all patients received flurbiprofen as a part of standardized multimodal protocol, which may have mitigated differences in baseline pain scores between the two groups. Future studies should evaluate the additive effects of esketamine in opioid-sparing regimens or in populations with contraindications to NSAIDs.

We also observed that the incidence of dizziness was higher in the treatment group than in the control group. Nevertheless, our data showed no clinically significant differences in QoR-15 scores between the two groups. A Meta-analysis, including 18 RCTs, has demonstrated that the use of ketamine or esketamine in the perioperative period does not increase the risk of adverse events.²⁹ The higher incidence of dizziness in the esketamine group underscores the need for cautious dose titration. Notably, dizziness did not persist in POD3 or achieve clinically meaningful differences in QoR-15 scores, suggesting it is transient and manageable symptom. This is consistent with meta-analysis indicating that perioperative ketamine does not increase severe adverse events. However, clinicians should monitor neuropsychiatric symptoms, particularly in vulnerable populations (eg, the elderly or those with a psychiatric history.^30,31

Our study has several limitations. First, the single-center design and moderate sample size limit its generalizability. Although our power calculation was robust for the primary endpoint, larger multicenter trial is needed to validate these findings across diverse surgical populations. Second, the fixed timing of analgesic pump removal may have introduced bias, as individual variability in pain trajectories could influence nocturnal pain severity. Adaptive protocols allowing pump removal based on pain thresholds rather than fixed durations might better capture the true efficacy of epidural esketamine. Third, we did not measure plasma inflammatory markers and esketamine levels, which could provide mechanistic insights into interindividual variability. Future research should incorporate pharmacokinetic-pharmacodynamic modeling to identify optimal dosing schedules and to guide biomarker-informed analgesia. While transient dizziness was more frequent with esketamine (47.7%), its resolution by POD3 and absence of neuropsychiatric events support epidural safety. This contrasts with IV esketamine studies reporting 20–30% hallucination rates, suggesting route-specific advantages.

Conclusion

In summary, the addition of esketamine to PCEA significantly reduced nocturnal pain on POD2 during which analgesia device has been removed, and improved moderate pain, as assessed by the QoR-15 in patients undergoing open gynecological surgery. Moreover, the use of epidural esketamine at a dosage of 0.4 mg/kg was found to be safe and feasible, without clinically significant adverse effects. Based on our safety findings and dose-response trends, it is necessary to investigate dose escalation of esketamine. While higher doses may prolong the duration of pain relief, they may also increase adverse effects. Further studies utilizing an up-and-down biased coin design are required to determine the optimal dosage and associated adverse events. Alternatively, investigating a multimodal strategy combining low-dose esketamine and gabapentinoids to simultaneously target spinal NMDA receptors and peripheral sensitization pathways may enhance efficacy for gynecological visceral pain.

Abbreviations

EG, esketamine group; CG, control group; ICU, Intensive Care Unit; POD, postoperative day; NRS, Numeric Rating Scale; QoR-15, Quality of Recovery-15; PONV, Postoperative Nausea and Vomiting; BMI, body mass index; ASA: American Society of Anesthesiologists.

Data Sharing Statement

Deidentified individual participant data will be made available. The shared data will include CRFs and Excel-formatted datasets. Data is available on request from the corresponding author. Study protocol, statistical analysis plan, and informed consent template will be accessible. The data will become publicly available immediately upon online publication of the article and will remain accessible for 10 years.

Ethics Approval and Informed Consent

Ethical approval was obtained from the ethics committees of the Fudan University Shanghai Cancer Center (No: 2407300-10). Written informed consent was obtained from all participants. The study complies with the Declaration of Helsinki.

Acknowledgment

Huanyu Luo and Yuecheng Yang are co-first authors for this study. We thank all the participants who conducted and assisted in this work.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This study was supported by Shanghai Science and Technology Commission (No. 22Y11904200).

Disclosure

The authors report no conflicts of interest in this work.

References

1. Liang B, Tan J, Li J. et al. Epidemiology, molecular typing, microbiome-immune interactions and treatment strategies of endometrial cancer: a review. Front Immunol. 2025;16:1595638. doi:10.3389/fimmu.2025.1595638

2. Guerra RA, Tucker LY, Littell R, Kay AH. Long-term impact of surgical route and tumor size on risk of recurrence among early-stage cervical cancer patients in a managed care population. Int J Gynecol Cancer. 2025;35(7):101928. doi:10.1016/j.ijgc.2025.101928

3. Rawal N. Current issues in postoperative pain management. Eur J Anaesthesiol. 2016;33(3):160–171. doi:10.1097/EJA.0000000000000366

4. Gitman M, Barrington MJ. Local Anesthetic Systemic Toxicity: a Review of Recent Case Reports and Registries. Reg Anesth Pain Med. 2018;43(2):124–130. doi:10.1097/AAP.0000000000000721

5. Farmawy MSE, Mowafy SMS, Wahdan RA. Epidural nalbuphine versus dexmedetomidine as adjuvants to bupivacaine in lower limb orthopedic surgeries for postoperative analgesia: a randomized controlled trial. BMC Anesthesiol. 2023;23(1):401. doi:10.1186/s12871-023-02348-x

6. Zhang X, Wang D, Shi M, Luo Y. Efficacy and Safety of Dexmedetomidine as an Adjuvant in Epidural Analgesia and Anesthesia: a Systematic Review and Meta-analysis of Randomized Controlled Trials. Clin Drug Investig. 2017;37(4):343–354. doi:10.1007/s40261-016-0477-9

7. Feltracco P, Barbieri S, Rizzi S, et al. Brief report: perioperative analgesic efficacy and plasma concentrations of S+ -ketamine in continuous epidural infusion during thoracic surgery. Anesth Analg. 2013;116(6):1371–1375. doi:10.1213/ANE.0b013e31828cbaf0

8. Yan H, Chen W, Chen Y, et al. Opioid-Free Versus Opioid-Based Anesthesia on Postoperative Pain After Thoracoscopic Surgery: the Use of Intravenous and Epidural Esketamine. Anesth Analg. 2023;137(2):399–408. doi:10.1213/ANE.0000000000006547

9. Zhang A, Zhou Y, Zheng X, et al. Effects of S-ketamine added to patient-controlled analgesia on early postoperative pain and recovery in patients undergoing thoracoscopic lung surgery: a randomized double-blinded controlled trial. J Clin Anesth. 2024;92:111299. doi:10.1016/j.jclinane.2023.111299

10. Gastaldon C, Raschi E, Kane JM, Barbui C, Schoretsanitis G. Post-Marketing Safety Concerns with Esketamine: a Disproportionality Analysis of Spontaneous Reports Submitted to the FDA Adverse Event Reporting System. Psychother Psychosom. 2021;90(1):41–48. doi:10.1159/000510703

11. Xu LL, Wang C, Deng CM, et al. Efficacy and Safety of Esketamine for Supplemental Analgesia During Elective Cesarean Delivery: a Randomized Clinical Trial. JAMA Network Open. 2023;6(4):e239321. doi:10.1001/jamanetworkopen.2023.9321

12. Zhu Y, Li Q, Liu G, et al. Effects of esketamine on postoperative rebound pain in patients undergoing unilateral total knee arthroplasty: a single-center, randomized, double-blind, placebo-controlled trial protocol. Front Neurol. 2023;14:1179673. doi:10.3389/fneur.2023.1179673.

13. Xie M, Liang Y, Deng Y, Li T. Effect of S-ketamine on Postoperative Pain in Adults Post-Abdominal Surgery: a Systematic Review and Meta-analysis. Pain Physician. 2023;26(4):327–335. doi:10.36076/ppj.2023.26.327

14. Woolf CJ. Central sensitization: implications for the diagnosis and treatment of pain. Pain. 2011;152(3 Suppl):S2–S15. doi:10.1016/j.pain.2010.09.030

15. Latremoliere A, Woolf CJ. Central sensitization: a generator of pain hypersensitivity by central neural plasticity. J Pain. 2009;10(9):895–926. doi:10.1016/j.jpain.2009.06.012

16. Touil N, Pavlopoulou A, Barbier O, Libouton X, Lavand’homme P. Evaluation of intraoperative ketamine on the prevention of severe rebound pain upon cessation of peripheral nerve block: a prospective randomised, double-blind, placebo-controlled study. Br J Anaesth. 2022;128(4):734–741. doi:10.1016/j.bja.2021.11.043

17. Mion G, Himmelseher S. Esketamine: less Drowsiness, More Analgesia. Anesth Analg. 2024;139(1):78–91. doi:10.1213/ANE.0000000000006851

18. Zeng X, Zhang X, Jiang W, Zhou X. Efficacy of Intravenous Administration of Esketamine in Preventing and Treating Rebound Pain After Thoracic Paravertebral Nerve Block: a Prospective Randomized, Double-Blind, Placebo-Controlled Trial. Drug Des Devel Ther. 2024;18:463–473. doi:10.2147/DDDT.S448336

19. Li S, Zhuo Z, Li R, Guo K. Efficacy of esketamine for the treatment of postpartum depression and pain control following cesarean section: a randomized, double-blind, controlled clinical trial. BMC Anesthesiol. 2024;24(1):52. doi:10.1186/s12871-024-02436-6

20. Liu J, Wang T, Song J, Cao L. Effect of esketamine on postoperative analgesia and postoperative delirium in elderly patients undergoing gastrointestinal surgery. BMC Anesthesiol. 2024;24(1):46. doi:10.1186/s12871-024-02424-w

21. Himmelseher S, Ziegler-Pithamitsis D, Argiriadou H, Martin J, Jelen-Esselborn S, Kochs E. Small-dose S(+)-ketamine reduces postoperative pain when applied with ropivacaine in epidural anesthesia for total knee arthroplasty. Anesth Analg. 2001;92(5):1290–1295. doi:10.1097/00000539-200105000-00040

22. Pirie K, Traer E, Finniss D, Myles PS, Riedel B. Current approaches to acute postoperative pain management after major abdominal surgery: a narrative review and future directions. Br J Anaesth. 2022;129(3):378–393. doi:10.1016/j.bja.2022.05.029

23. Rosenberg-Adamsen S, Kehlet H, Dodds C, Rosenberg J. Postoperative sleep disturbances: mechanisms and clinical implications. Br J Anaesth. 1996;76(4):552–559. doi:10.1093/bja/76.4.552

24. Song B, Zhu J. A Novel Application of Ketamine for Improving Perioperative Sleep Disturbances. Nat Sci Sleep. 2021;13:2251–2266. doi:10.2147/NSS.S341161

25. Yoon JP, Kim HY, Jung J, Lee J, Park S, Byeon GJ. Analgesic effect of ultrasound-guided transversus abdominis plane block with or without rectus sheath block in laparoscopic cholecystectomy: a randomized, controlled trial. BMC Anesthesiol. 2024;24(1):203. (). doi:10.1186/s12871-024-02590-x

26. Kehlet H. Multimodal approach to control postoperative pathophysiology and rehabilitation. Br J Anaesth. 1997;78(5):606–617. doi:10.1093/bja/78.5.606

27. Dieu A, Huynen P, Lavand’homme P, et al. Pain management after open liver resection: procedure-Specific Postoperative Pain Management (PROSPECT) recommendations. Reg Anesth Pain Med. 2021;46(5):433–445. doi:10.1136/rapm-2020-101933

28. Zhang Y, Cui F, Ma JH, Wang DX. Mini-dose esketamine-dexmedetomidine combination to supplement analgesia for patients after scoliosis correction surgery: a double-blind randomised trial. Br J Anaesth. 2023;131(2):385–396. doi:10.1016/j.bja.2023.05.001

29. Hung KC, Kao CL, Ho CN, et al. The impact of perioperative ketamine or esketamine on the subjective quality of recovery after surgery: a meta-analysis of randomised controlled trials. Br J Anaesth. 2024;132(6):1293–1303. doi:10.1016/j.bja.2024.03.012

30. Zanos P, Moaddel R, Morris PJ, et al. Ketamine and Ketamine Metabolite Pharmacology: insights into Therapeutic Mechanisms. Pharmacol Rev. 2018;70(3):621–660. doi:10.1124/pr.117.015198

31. Kurdi MS, Theerth KA, Deva RS. Ketamine: current applications in anesthesia, pain, and critical care. Anesth Essays Res. 2014;8(3):283–290. doi:10.4103/0259-1162.143110