Comprehensive analysis of tRNA-derived fragment expression in endometriosis using PANDORA-seq technology

Endometriosis (EMs) is a chronic systemic disease characterized by the growth of endometrioid tissue outside the uterus with tumour-like characteristics, affecting 5%-10% of women of reproductive age worldwide [ 1 ]. The gold standard for the diagnosis of endometriosis is still laparoscopic observation of the lesion and histopathological confirmation [ 2 ]. Usually, the time between the onset of symptoms and the definitive diagnosis can reach 7–9 years [ 3 ]. This diagnostic delay is related to the invasiveness and complexity of the laparoscopic surgery required for diagnosis, as well as to the variability of disease symptoms and the overlap of symptoms with other diseases and the lack of specific markers, among other factors [ 4 ]. Therefore, it is urgent to explore specific and sensitive non-invasive biomarkers that can predict or diagnose EMs early. This biomarker research has been identified as a research focus by the World Endometriosis Society (WES) and the World Endometriosis Research Foundation (WERF) [ 5 ]. Previous studies have shown that noncoding RNAs (ncRNAs) play a key role in the occurrence and development of endometriosis [ 6 , 7 ]. With the development of next-generation sequencing technology, tRNA-derived tsRNAs have attracted increasing attention. Recent studies have shown that tsRNAs are a class of small noncoding RNAs with posttranscriptional gene regulation functions, mainly including tRNA-derived fragments (tRFs) and tRNA-derived stress-induced RNAs (tiRNAs) [ 8 , 9 ]. Many studies have shown that tsRNAs are multifunctional regulators related to various physiological and pathological states and can even serve as the aetiology of diseases such as cancer and inherited metabolic disorders [ 9 , 10 ]. Numerous studies have demonstrated that the dysregulation of tsRNAs is involved in various pathological processes, including cancer, metabolic disorders, and neurological diseases, positioning tsRNAs as ideal biomarkers for the diagnosis and prognosis of these conditions [ 11 , 12 ]. Specifically, the use of tRFs as biomarkers in cancer was first reported in 2015, with further developments highlighted in 2019 [ 13 , 14 ]. These contributions have significantly advanced the field. For example, tRNA-derived fragment dysregulation has been shown to be related to the pathogenesis and progression of ovarian endometriosis (EMs) [ 15 ], and tRFs/tiRNAs have proven to be potential biomarkers for the diagnosis and treatment of EMs in the future [ 16 ]. However, despite these advances, the function of tRFs in endometriosis remains largely unexplored. Compared to miRNAs, tsRNAs and rsRNAs are rich in RNA modifications due to the various RNA modifications of their precursor RNAs, such as internal modifications (e.g. m1G, m1A, m3C, N2-dimethylguanosine, and m22G) [ 17–21 ] and terminal modifications (e.g. 3′-phosphate and 2′,3′-cyclic phosphate) [ 22 , 23 ]. When standard RNA-seq analysis for highly modified sncRNAs, some RNA modifications interfere with library construction processes such as adaptor ligation and reverse transcription, which prevents discovery of highly modified tsRNAs that may be responsible for endometriosis [ 24 ]. The removal of internal RNA modifications prior to sequencing is crucial for efficient RNA analysis [ 22 , 25 , 26 ]. PANDORA-seq addresses the sequencing challenges arising from RNA modifications by employing T4PNK and AlkB enzymatic treatment on small RNAs (15-50nt). This process rectifies ligation issues due to modifications, converting 3’phosphates or 2’,3’−cyclic phosphates to 3“-OH and 5”−OH to 5’ phosphates, and removes specific methylation modifications (m1A, m1G, m3C, m22G) that impede reverse transcriptase activity [ 20 ]. In this study, PANDORA-Seq technology was used to reveal the tRF expression profiles of the ectopic endometrium, eutopic endometrium and normal eutopic endometrium in the cyst wall of patients with ovarian endometriosis ( Figure 1 ). Through screening and verification, it was found that differential tRF is expected to become a new biomarker for endometriosis.These findings may provide new insights into the molecular mechanism of EMs and future therapeutic approaches. Figure 1. Technology roadmap. EN: eutopic endometrial tissue of the normal uterus; EC: EMs ectopic endometrial tissue; EU: EMs eutopic endometrial tissue; tRF: tRNA-derived fragment; FC: Fold change; qRT-PCR: quantitative real-time polymerase chain reaction. Technology roadmap. EN: eutopic endometrial tissue of the normal uterus; EC: EMs ectopic endometrial tissue; EU: EMs eutopic endometrial tissue; tRF: tRNA-derived fragment; FC: Fold change; qRT-PCR: quantitative real-time polymerase chain reaction.

This study was approved by the Ethics Committee of the People’s Hospital of Guangxi Zhuang Autonomous Region, and all participants provided informed consent. In the experimental group, patients who were diagnosed with ovarian endometriosis by biopsy after laparoscopic surgery at the People’s Hospital of Guangxi Zhuang Autonomous Region between November 2021 and December 2022 were included. The control group included patients who underwent hysteroscopy and laparoscopic surgery due to fallopian tube obstruction at the People’s Hospital of Guangxi Zhuang Autonomous Region during the same period. Eutopic endometrial tissues (EUs) from patients with ovarian endometriosis (25–35 years old), matched ectopic endometrial tissues (ECs) from the same patient and eutopic endometrial tissues (ENs) from control patients (25–35 years old) were collected. According to the American Society for Reproductive Medicine (ASRM) classification of endometriosis, patients were classified as moderate to severe (III-IV) stage. All patients in the EMs group were diagnosed with grade III-IV endometriosis, and all were in the endometrial proliferative phase when the tissues were collected. All patients had a normal menstrual cycle (21–35 days) and had not received gonadotropin-releasing hormone analogs or other hormonal medications for at least 6 months before surgery. All tissue samples were removed from the body, placed in RNAlater solution within 20 minutes ex vivo, stored at 4°C overnight, and then frozen at −80°C until use. Clinical and demographic data of EU, EC, and EN donors are provided in Supplementary Table S1. For RNA isolation, 1 mL of TRIzol reagent (Invitrogen; catalogue number 15,596,018) was added to a microtube containing 50–100 mg of finely pulverized tissue. The mixture was then allowed to solubilize for 5 minutes at ambient temperature. Subsequently, 200 μL of chloroform was added per mL of sample, followed by vigorous vortexing for 15 seconds. The samples were incubated for 2–3 minutes at room temperature before being centrifuged at 10,000 g for 10 minutes at 4°C. The supernatant was carefully transferred to a new tube and mixed gently with an equal volume of isopropanol. After a 10-minute incubation, the mixture was centrifuged at 10,000 g for 10 minutes at 4°C. The supernatant was aspirated, and the RNA pellet was washed with 1 mL of 75% ethanol, followed by centrifugation at 7,500 g for 5 minutes at 4°C. The ethanol was decanted, and the pellet was air-dried for 5 minutes before being resuspended in nuclease-free water. For the enzymatic treatment, 200 ng of RNA was incubated in a 50 μL reaction mixture composed of 50 mM HEPES (pH 8.0), 75 μM ferrous ammonium sulphate (pH 5.0), 2 mM sodium ascorbate, 1 mM α-ketoglutaric acid, 50 mg/L bovine serum albumin, 4 μg/mL AlkB (Epibiotek; catalogue number R1822), and 2,000 U/mL RNase inhibitor at 37°C for 30 minutes. Following this incubation, the RNA isolation was performed using 500 μL of TRIzol reagent. The RNA was then incubated in a fresh 50 μL reaction mixture containing 5 μL of 10× PNK buffer (New England Biolabs; catalogue number B0201S), 10 mM ATP (New England Biolabs; catalogue number P0756S), and 10 U of T4PNK (New England Biolabs; catalogue number M0201L) at 37°C for 20 minutes. The RNA was subsequently isolated once more using 500 μL of TRIzol reagent. Small RNA library construction was initiated by resolving RNA segments through polyacrylamide gel electrophoresis (PAGE), followed by the precise isolation of a 15- to 45-nucleotide band. Adapters were procured from the QIAseq® miRNA Library Kit (QIAGEN; catalogue number 331,505) and sequentially ligated to the purified RNA fragments. Quality control (QC) of the library was conducted using the Bioptic Qsep100 Analyzer to ensure the integrity and purity of the library. Subsequently, the amplified flow cell was subjected to single-end 75-cycle (SE75) sequencing on an Illumina platform by Epibiotek (Guangzhou Epibiotek Co., Ltd., Guangzhou, China) to achieve deep sequencing coverage. The raw sequencing reads underwent rigorous quality control processing to meet the following criteria: (1) exclusion of reads with ambiguous nucleotides (N); (2) removal of reads with more than four consecutive bases having a Phred quality score < 10; (3) elimination of reads with more than six consecutive bases having a Phred quality score < 13; (4) filtering out reads containing 5’primer contaminants or lacking a 3’ primer; (5) discarding reads that do not possess the insert tag; (6) removal of reads with poly-A sequences; and (7) exclusion of reads shorter than 15 nucleotides or longer than 44 nucleotides. The subsequent sequencing data analysis was conducted on the high-quality, filtered reads to ensure the accuracy and reliability of downstream bioinformatics analyses. Trimmomatic, as described by Bolger et al. (2014), was employed to excise adapter sequences and to discard reads of substandard quality. The resultant high-fidelity reads were then mapped against the reference sequences in miRBase 22 (Kozomara et al., 2019). Unmapped reads within the size range of 24–33 bp were subsequently aligned to the sequences in NCBI and the GtRNA database (Chan & Lowe, 2016), which is dedicated to tRNAs derived from piRNAs. After the removal of rRNA, the remaining unmapped reads were aligned to the sequences present in tsRBase (Zuo et al., 2021). For the alignment and annotation of tsRNAs, SPORTS 1.1 (Shi et al., 2018) was utilized, leveraging its precompiled annotation databases, with a one mismatch tolerance (SPORTS1.1 parameter setting: -M 1). The R package DESeq2 was used for gene differential expression analysis, and tRFs with a |Fold Change| ≥2 and p  < 0.05 were defined as significantly differentially expressed tRFs. R software was used to perform hierarchical clustering and correlation analysis of the differentially expressed tRFs, and hierarchical cluster plots and volcano plots were drawn to facilitate visual identification. To validate the results of PANDORA-seq, we selected an additional 15 tissue samples from the EN, EU, and EC groups, with 5 independent samples per group, to perform qRT-PCR experiments on 13 differentially expressed tRFs. The inclusion and exclusion criteria for sample selection, as well as the methods of sample collection and preservation, are consistent with those used for the PANDORA-seq samples. U6 was used as the reference gene, and Primer5.0 software was used to design primers ( Table 1 ). All primers were synthesized by Sangon Biotech (Shanghai, China). For the reverse transcription of RNA in preparation for qPCR, we used the QIAseq® First-Strand cDNA Synthesis Kit from Arraystar. Starting with 5 micrograms of total RNA, we performed a series of pretreatment steps including 3’ end deacetylation, removal of 3’-cP, addition of 5’-P, demethylation treatment, and ligation of 3’ and 5’ adapters. The reverse transcription reaction was conducted at 50°C for 1 hour to yield cDNA. The 10 μL PCR system consisted of t2 μL cDNA, 5 μL 2x Master premix, 1 μL qPCR primer mixture, and 2 μL water without ribonuclease. For qPCR, we used 2 microlitres of cDNA as the template, which corresponds to a small fraction of the initial RNA quantity to ensure the efficiency and sensitivity of the PCR reaction.The qRT-PCR conditions were as follows: predenaturation at 95 °C for 10 min, followed by 40 cycles of amplification at 95 °C for 30 secures, followed by amplification at 95 °C for 10 sec and 60 °C for 60 sec. Table 1. Primer sequences. Gene name Primer sequence Product length (bp) U6 F:5’GCTTCGGCAGCACATATACTAAAAT3’ 89   R:5’CGCTTCACGAATTTGCGTGTCAT3’   tDR-42:56-Val-AAC-1-M9 F:5’ CTACAGTCCGACGATCAAAGGTC 3’ 45   R:5’ TGCTCTTCCGATCTGAACCG3’   tRF-16-F0S7RND F:5’ GTCCGACGATCAGAAGGTCCT3’ 41   R:5’ TGCTCTTCCGATCTGAACCC 3’   tRF-32-I2VHJYOZ4E8B4 F:5’ TCTACAGTCCGACGATCATGG3’ 55   R:5’ CGATCTCCTTTGAATTAGAAGTCC3’   tRF-19-O6M83OJX F:5’ ACAGTCCGACGATCGACCAC3’ 45   R:5’ CTCTTCCGATCTTCCATTAGGC3’   tRF-19-OR183OJX F:5’ ACAGTCCGACGATCGACCCA3’ 45   R:5’ CTCTTCCGATCTTCCATTAGGC 3’   tRF-17-QB1MK8Q F:5’ CTACAGTCCGACGATCGCG3’ 41   R:5’ TCCGATCTAACCGGGGACC3’   tRF-19-R118LOJX F:5’ CTACAGTCCGACGATCGGG3’ 39   R:5’ ATCTTCCATTGCGCCACTG3’   tRF-20-R2IP4OI3 F:5’ TCCGACGATCGGGGAATTA3’ 44   R:5’ TGCTCTTCCGATCTACCATTTGA3’   tDR-42:56-Arg-CCG-2 F:5’ TCCGACGATCGGGGATTGT3’ 45   R:5’ GACGTGTGCTCTTCCGATCTGA 3’   tRF-25-R8YP9LON4V F:5’ ACGATCGGTGGTTCAGTGGT 3’ 44   R:5’ GCTCTTCCGATCTGCGAGAA3’   tRF-20-UDJKSV5N F:5’ CCGACGATCTAAACCAGGGG3’ 47   R:5’ CGTGTGCTCTTCCGATCTAACTC3’   tRF-19-86J8WP1Z F:5’ AGTCCGACGATCTCCCACAT3’ 44   R:5’ GCTCTTCCGATCTACCGCTA3’   tRF-21-W3FJHPEZE F:5’ CGACGATCTCGGTAGAGCATC 3’ 46   R:5’ GTGTGCTCTTCCGATCTAAAAGTC3’ Primer sequences. The TargetScan and miRanda algorithms were used to predict the target genes of the tRFs, and the common genes were used for subsequent analysis [ 27 ]. We employed miRanda V3.3 and TargetScan 6.0 for the identification of potential target genes. The specific parameters for the human species used in these software versions are as follows: For miRanda, we set the parameters to miranda_score >140 and miranda_energy < −10. For TargetScan 6.0, we used the parameter context_plus_score < −0.4. The Visualization and Integrated Discovery database (Metascpe) was used for annotated gene ontology (GO) functional analysis. For the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, we utilized Kobas software and set the P -value cut-off at 0.05 to identify significantly enriched pathways. The functions of the target genes were inferred based on these enriched pathways. Using Fisher’s exact test to evaluate the results, p  < 0.05 considered as significant. The 2^(-ΔΔct) method was used to calculate the gene expression levels in the endometrium of each group. All experimental data were analysed using Prism 9.0 statistical analysis software (GraphPad, USA) and are presented as the mean ± standard error (SEM) of at least three independent experiments. Unpaired t tests were used to compare means from 2 independent groups, and one-way analysis of variance (ANOVA) was used to compare means from 2 or more samples. A p value less than 0.05 was considered to indicate statistical significance.

To clarify the expression profiles of tRFs in the eutopic and ectopic endometria of patients with EMs and to explore whether the differential expression levels of tRFs are related to the occurrence and progression of EMs. Using PANDORA-seq, we sequenced eutopic and ectopic endometrial tissues from patients with endometriosis, as well as normal tissues from controls. The illustration of tsRNA length distribution and sequencing read count across the samples is presented in the subsequent figure ( Figure 2a ). The base length of the tRFs fragment was enriched in the ranges of 15–21 nt and 32–34 nt, and the normal tissues and the eutopic and ectopic endometrial tissues of EMs patients showed two peaks, among which the EC group had the largest read count. The variance analysis showed significant differences between EC group and EU group at 16 nt ( Figure 2b ).The distribution of tRFs and tiRNAs (mature tRNA and precursor tRNA) in each group of samples ( Figure 2c ). The sequencing results of the three groups of samples showed that the expression of tRNA-Gly-GCC was the highest in the normal group (EN), the disease group (EU), and the EC group, and the read value was the highest in the ectopic endometrial tissue EC group ( Figure 2d ). In addition, sequencing revealed that the types of tRNA in the EN, EU, and EC groups were similar, and the expression level of tRNA-Gly was the highest among the three groups ( Figure 2e-g ). Figure 2. Length and species distribution characteristics of the tRFs. (a) Length distribution of the identified tRFs. X axis shows the length of the detected tRFs and Y axis shows the abundance of tRFs classified by different lengths. (b) The length of tRFs in endometrial tissues(EU 'EC) and normal tissues (EN) mainly ranged from 15 to 45 nt, ** p < 0.01. (c) Reads the source number of tRFs and tiRNAs (precursor and mature tRNA). (d) species distribution characteristics among the three groups. (e-g) radar chart of pairwise species comparisons among the three groups. EN group : normal eutopic endometrial tissue; EU group: eutopic endometrial tissue with endometriosis; EC group: endometriosis ectopic endometrial tissue. Length and species distribution characteristics of the tRFs. (a) Length distribution of the identified tRFs. X axis shows the length of the detected tRFs and Y axis shows the abundance of tRFs classified by different lengths. (b) The length of tRFs in endometrial tissues(EU 'EC) and normal tissues (EN) mainly ranged from 15 to 45 nt, ** p < 0.01. (c) Reads the source number of tRFs and tiRNAs (precursor and mature tRNA). (d) species distribution characteristics among the three groups. (e-g) radar chart of pairwise species comparisons among the three groups. Differential expression analysis by DEseq2, |Fold Change| ≥2 and p  < 0.05 were selected as the gene screening conditions. A total of 1045 tRFs were differentially expressed in EN and EU groups, of which 606 tRFs were upregulated and 439 were downregulated ( Figure 3a,d ); 1929 tRFs were differentially expressed in EN and EC groups, of which 1374 tRFs were upregulated and 555 were downregulated ( Figure 3b,e ). In the EU and EC groups, 1841 tRFs were differentially expressed, of which 987 tRFs were upregulated and 854 tRFs were downregulated ( Figure 3c,e ). Figure 3. Pairwise comparison map of differentially expressed tRFs among the EN, EU and EC groups. Hierarchical clustering analysis of differentially expressed tRFs between EN and EU (a), EN and EC (b), and EU and EC (c). Volcano plots of the differentially expressed tRFs in the EN vs EU (d), EN vs EC (e), and EU vs EC (f) comparisons. Red/blue dots indicate differentially expressed tRFs with |fold change| ≥ 2 and p < 0.05, and gray dots indicate tRFs that are not differentially expressed. EN group : normal eutopic endometrial tissue; EU group: eutopic endometrial tissue with endometriosis; EC group: endometriosis ectopic endometrial tissue; tRF: tRNA-derived fragment. Pairwise comparison map of differentially expressed tRFs among the EN, EU and EC groups. Hierarchical clustering analysis of differentially expressed tRFs between EN and EU (a), EN and EC (b), and EU and EC (c). Volcano plots of the differentially expressed tRFs in the EN vs EU (d), EN vs EC (e), and EU vs EC (f) comparisons. Red/blue dots indicate differentially expressed tRFs with |fold change| ≥ 2 and p < 0.05, and gray dots indicate tRFs that are not differentially expressed. Under the conditions of a |Fold Change| ≥2 and padj  < 0.05, a total of 13 differentially expressed tRFs were screened compared with those in the normal EN group, of which 1 tRFs was downregulated and 12 tRFs were upregulated, as visualized by heatmap using the tDRnamer or ‘license plate’ nomenclature [ 28 , 29 ] ( Figure 4 ). 13 differentially expressed tRFs sequences, length, fold change ( Table 2 ). Secondary structure prediction of 13 differentially expressed tRFs ( Figure 5 ). Figure 4. Cluster plots and microarray data of 13 differentially expressed tRFs in the EU and EC groups of patients with EMs compared with the control group. The heatmap on the left shows the 13 differentially expressed tRFs clusters and their names, the colours in the heatmap indicate the relative expression levels, and the colour bars at the bottom show the different sample groups. The colour bars in the upper right corner indicate the different groups in the microarray data. Compared with that in the normal EN group, the expression of tRF-32-I2VHJYOZ4E8B4 was downregulated in the EU and EC groups. CPM: counts per million; tRF: tRNA-derived fragment. Figure 5. Secondary structure prediction of 13 differentially expressed tRFs. Predicting the secondary structure of tRNA fragment by RNAfold web server. Table 2. 13 differentially expressed tRFs sequences, length, Fold change table. tRF ID Fragment sequence,5–3 nt EN VS EU log2Fold Change EN VS EC log2Fold Change tDR-42:56-Val-AAC-1-M9 AAAGGTCCCCGGTTC 15 4.780447832 6.431769589 tRF-16-F0S7RND AGAAGGTCCTGGGTTC 16 3.348778916 3.450797131 tRF-19-O6M83OJX GACCACGTGGCCTAATGGA 19 6.082307411 4.190300591 tRF-19-OR183OJX GACCCAGTGGCCTAATGGA 19 2.895833683 3.078061994 tRF-17-QB1MK8Q GCGAAAGGTCCCCGGTT 17 3.14329591 3.092659508 tRF-19-R118LOJX GGGCCAGTGGCGCAATGGA 19 6.760140778 2.998554924 tRF-20-R2IP4OI3 GGGGAATTAGCTCAAATGGT 20 4.182494372 3.216985697 tDR-42:56-Arg-CCG-2 GGGGATTGTGGGTTC 15 5.926475143 3.406450013 tRF-25-R8YP9LON4V GGTGGTTCAGTGGTAGAATTCTCGC 25 4.21348835 2.974723376 tRF-20-UDJKSV5N TAAACCAGGGGTCGCGAGTT 20 4.062751167 3.489898714 tRF-19-86J8WP1Z TCCCACATGGTCTAGCGGT 19 2.923786557 3.219418297 tRF-21-W3FJHPEZE TCGGTAGAGCATCAGACTTTT 21 4.635682851 4.076338412 tRF-32-I2VHJYOZ4E8B4 ATGGATAGCGCATTGGACTTCTAATTCAAAGG 32 −3.59053481 −2.561654689 Cluster plots and microarray data of 13 differentially expressed tRFs in the EU and EC groups of patients with EMs compared with the control group. The heatmap on the left shows the 13 differentially expressed tRFs clusters and their names, the colours in the heatmap indicate the relative expression levels, and the colour bars at the bottom show the different sample groups. The colour bars in the upper right corner indicate the different groups in the microarray data. Compared with that in the normal EN group, the expression of tRF-32-I2VHJYOZ4E8B4 was downregulated in the EU and EC groups. Secondary structure prediction of 13 differentially expressed tRFs. Predicting the secondary structure of tRNA fragment by RNAfold web server. 13 differentially expressed tRFs sequences, length, Fold change table. To verify the sequencing results, qRT-PCR was performed on the 13 selected DE tRFs. The qRT-PCR results showed that 11 tRFs were upregulated and 1 tRF was downregulated in the EU group compared with the EN group, which was consistent with the sequencing results. Compared with those in the EN group, 9 tRFs were upregulated and 1 tRF was downregulated in the EC group, which was consistent with the sequencing results ( Figure 6 ). Figure 6. qRT-PCR verification of the abnormal expression of 13 differentially expressed tRFs in the EU and EC groups of EMs patients and the control group of EN patients. * indicates a significant difference between the two groups (* p < 0.05, ** p < 0.01, *** p < 0.001); ns indicates no significant difference between the two groups. EN group : normal eutopic endometrial tissue; EU group: eutopic endometrial tissue with endometriosis; EC group: endometriosis ectopic endometrial tissue; tRF: tRNA-derived fragment. qRT-PCR verification of the abnormal expression of 13 differentially expressed tRFs in the EU and EC groups of EMs patients and the control group of EN patients. * indicates a significant difference between the two groups (* p < 0.05, ** p < 0.01, *** p < 0.001); ns indicates no significant difference between the two groups. The tsRNAs and miRNAs have similar functions. TsRNAs regulate the stability of mRNA by binding to it, thereby inhibiting translation and regulating multiple mechanisms, such as gene expression, the cell cycle, chromatin and epigenetic modifications [ 20 ]. The aim of this study was to further explore the role of DE tRFs in endometriosis.TargetScan and miRanda were used to predict the target genes of the 11 differentially expressed tRFs with lengths ≥ 16 nucleotides (nt), and the common target genes predicted by the two methods were used for subsequent functional analysis. A total of 1036 target genes were predicted by 11DE tRFs ( Figure 7 ). Figure 7. Venn diagram showing the predicted target genes of the differentially expressed tRFs. The two prediction software programs simultaneously predicted 1036 mRNA targets of 11 differentially expressed tRFs. Venn diagram showing the predicted target genes of the differentially expressed tRFs. The two prediction software programs simultaneously predicted 1036 mRNA targets of 11 differentially expressed tRFs. The GO enrichment results showed that the molecular functions of the 11 DE tRFs target genes were mainly related to positive regulation of apoptotic signalling pathway, positive regulation of apoptotic signalling pathway, Wnt signalling pathway, ubiquitin protein ligase binding, acyltransferase activity, transcription coregulator activity, etc ( Figure 8 ). KEGG analysis revealed that the spliceosome, mitophagy, MAPK, Ras and p53 signalling pathways were significantly enriched. These pathways play important roles in the occurrence and development of EMs ( Figure 9 ) . Figure 8. GO analysis of 11 differentially expressed DE tRFs target genes; bubble plots are aggregated by molecular function, cellular component, and biological process. The vertical axis shows the annotated functions of the target genes. The horizontal axis shows the proportion of enriched genes. The right sbar SampleGroup shows that BP, CC, and MF represent molecular function, cellular component, and biological process, respectively. Count shows the number of genes enriched in each cluster, and use color swatches to show enrichment scores (-log10 transformed q-value). Figure 9. Bubble plot formed by KEGG analysis aggregation of target genes of 11 differentially expressed tRFs. GO analysis of 11 differentially expressed DE tRFs target genes; bubble plots are aggregated by molecular function, cellular component, and biological process. The vertical axis shows the annotated functions of the target genes. The horizontal axis shows the proportion of enriched genes. The right sbar SampleGroup shows that BP, CC, and MF represent molecular function, cellular component, and biological process, respectively. Count shows the number of genes enriched in each cluster, and use color swatches to show enrichment scores (-log10 transformed q-value). Bubble plot formed by KEGG analysis aggregation of target genes of 11 differentially expressed tRFs.

In summary, this study used PANDORA-seq technology to analyse the tRF expression profiles of EU, EC and EN tissues, which is more in line with the different stages of the disease and restores the process of endometrial tissue metastasis and invasion to the greatest extent. In the future, the biological function of differential tRF will be further verified by functional experiments, and the pathogenesis of endometriosis will be further explored. It provides a new strategy for the study of non-invasive markers and new therapeutic targets for the diagnosis and prognosis evaluation of endometriosis.

Endometriosis is a systemic heterogeneous disease. It is urgent to explore specific and sensitive non-invasive biomarkers for the early prediction or diagnosis of EMs. TRFs play important roles in a variety of biological regulatory pathways and provide new ideas for the diagnosis and treatment of diseases. In this study, PANDORA sequencing technology was used to analyse the tRF expression profiles of EU, EC and EN tissues, which are more consistent with different stages of the disease, to identify new biomarkers or therapeutic targets for early clinical prediction and improve the specificity of clinical detection. A number of studies have shown that dysregulated tRFs are involved in the occurrence and development of cancers (breast cancer, ovarian cancer, lung cancer, prostate cancer, colorectal cancer, etc.), nervous system diseases, viral infectious diseases, metabolic diseases, and angiogenesis-related diseases [ 30 , 31 ]. In the reproductive system, dysregulated tsRNAs can affect gametocyte maturation, zygote activation, and early embryonic development, transmit epigenetic information, and alter the phenotype of offspring [ 32 ]. These studies reveal the promise of tsRNAs as disease diagnostic and prognostic biomarkers [ 33 ]. The tsRFun database evaluated the expression patterns and prognostic value of tRF in 32 human cancers, and demonstrated their biological functions using high-throughput CLASH/CLEAR or CLIP sequencing data [ 34 ]. By comparing the 13 differentially expressed tRFs sequences screened in this study with the tsRFun database, it was found that 4 differentially expressed tRFs (tRF-20-R2IP4OI3, tRF-19-O6M83OJX, tRF-19-OR183OJX, tRF-19-R118LOJX) were validated in various human tumour diseases ( Figures 10a-d ). In particular, tRF-17-QB1MK8Q was verified to be differentially expressed in human endometrial cancer. Figure 10. Box expression of four tRFs (tRF-20-R2IP4OI3, tRF-19-O6M83OJX, tRF-19-OR183OJX, tRF-19-R118LOJX) among the 13 differentially expressed tRFs in various human tumours ( Figure 10a-d ). Figure 10. (Continued). Box expression of four tRFs (tRF-20-R2IP4OI3, tRF-19-O6M83OJX, tRF-19-OR183OJX, tRF-19-R118LOJX) among the 13 differentially expressed tRFs in various human tumours ( Figure 10a-d ). (Continued). Through the analysis of the biological functions and pathways enriched by the target genes of the DE tRFs, it was found that the target genes of the DE tRFs were mainly enriched in positive regulation of apoptotic signalling pathway, positive regulation of apoptotic signalling pathway, Wnt signalling pathway. KEGG pathway analysis revealed that DE tRF was enriched mainly in the MAPK signalling pathway. Several studies have shown that the MAPK signalling pathway plays a key role in the pathophysiology of endometriosis [ 35 ]. In particular, lncRNA BANCR inhibitors may inhibit the development of ectopic endometrial tissue by inhibiting the ERK/MAPK signalling pathway and subsequently inhibiting the production of angiogenic factors in EMs lesions [ 36 ]. The overexpression of TGF-β also enhances the migration and invasion of ectopic ESCs through the ERK/MAPK signalling pathway [ 37 ]. The DE tRFs identified in this study may also participate in the occurrence and development of EMs diseases by participating in the upstream MAPK signalling pathway, which provides the latest research direction for the diagnosis and treatment of EMs. There are some limitations in this study. The tissue samples of the experimental group only included endometriosis ASRM III/IV cases. It is necessary to study the differential expression profiles of tRF in specimens before and after different stages of disease, surgery or drug inhibition treatment. In addition, due to the differences in the naming of tRF sequencing data in different studies, RNA sequence alignment has not been fully carried out, and the advantages of PANDORA-seq technology compared with traditional methods and newly discovered tRFs have not been shown. We identified two 15-nucleotide (nt) tsRNAs. Although 79% of tRNA-derived 15-mers map to non-canonical genomic loci [ 38 ], we retained these sequences in differential expression analysis given PANDORA-seq’s capability to resolve modified sequences. However, to mitigate analytical bias, they were excluded from target gene prediction and GO/KEGG enrichment analyses. Future studies should elucidate their biogenesis and functional roles. In conclusion, in this study, PANDORA-seq was used to detect the tRF expression profiles of eutopic and ectopic endometria from patients with endometriosis and of eutopic endometria from healthy controls, which best restored the process of endometrial metastasis and invasion, and four potential tRFs that could be used as biomarkers were screened. This study provides a new research strategy for the diagnosis, treatment and prognosis of endometriosis. In the future, the biological function of DE tRFs will be further verified by functional experiments, and the pathogenesis of endometriosis will be further explored, which will provide a new perspective for the development of non-invasive diagnostic markers and the elucidation of disease mechanisms of endometriosis.

Ems Endometriosis qRT-PCR Quantitative Real-time Polymerase Chain Reaction GO Gene Ontology KEGG Kycrto Encyclopedia of Genes and Genomes ncRNA Non-Coding RNA miRNA microRNA lncRNA Long noncoding RNA mRNA Messenger RNA sncRNA small ncRNA nt Nucleotide ASRM American Society for Reproductive Medicine Endometriosis Quantitative Real-time Polymerase Chain Reaction Gene Ontology Kycrto Encyclopedia of Genes and Genomes Non-Coding RNA microRNA Long noncoding RNA Messenger RNA small ncRNA Nucleotide American Society for Reproductive Medicine