Role of N6-methyladenosine Regulated lncRNAs in Prognosis, Tumor Microenvironment and Immune Checkpoints of Pancreatic Ductal Adenocarcinoma


 BackgroundThe purpose of this study was to explore the correlation between N6-methyladenosine (m6A)-regulated lncRNAs and tumor prognosis, immune infiltration, immune checkpoints (ICPs) expression in pancreatic ductal adenocarcinoma (PDAC).MethodsWe downloaded the raw RNA-sequence data and clinical data of PDAC from https://xenabrowser.net/ (cohort: TCGA Pancreatic Cancer) and Genotype-Tissue Expression project (GTEx). The m6A-regulated lncRNA was obtained by co-expression analysis. After that, lncRNA profiles and PDAC survival information were merged, and m6A-regulated multi-lncRNA prognostic model was constructed through least absolute shrinkage and selection operator (LASSO) analysis. Through consensus clustering algorithm analysis, PDAC samples were divided into C1 and C2 groups. The downstream pathway signals of the two groups were constructed by Gene set enrichment analysis (GSEA) analysis. Finally, we detect the links between m6A regulated lncRNAs, immune infiltration, immune checkpoint gene expression.ResultsA total of 28 differential expressed m6A-regulated lncRNAs were identiﬁed, and based on this, a total of two subtypes of PDAC were obtained. A risk score nomogram consist of 11 m6A-regulated lncRNAs was constructed based on LASSO regression analysis. PDAC patients were divided into low-risk and high-risk groups based on risk scores. In addition to that, we identified IDO1 as a potential novel ICPs in PDAC.ConclusionThis study demonstrates an indispensable role for m6A-regulated lncRNAs in the tumor microenvironment and immune infiltration. We could screen patients suitable for immunotherapy. Long term survival of PDAC patients can be predicted by 11 m6A regulated lncRNAs superiorly. The immune infiltration and ICPs expression were further explored in both groups.


Introduction
Pancreatic ductal adenocarcinoma (PDAC) is a notorious tumor of extremely high malignancy, which is the seventh leading cause of cancer death by 2020, with the highest incidence in Europe, North America, and Australia / New Zealand [1] . It has been pointed out that PDAC is projected to surpass breast cancer as the third most tumor lethal disease by 2025 [2] . The only curable treatment of PDAC is surgical resection, but this results in a 5-year overall survival rate of only 9% due to the fact that most patients have lost surgical opportunity when diagnosed with PDAC [3] . The advancement of adjuvant chemotherapy has modestly improved the long-term survival of these patients. However, mutations like KRAS and unique heterogeneous tumor microenvironment (TME) characterized by a bro-in ammatory stroma that contributes to disease progression, a part of PDAC patients suffered from chemotherapy resistance [4,5] . In recent years, immunotherapy has received extensive attention and become a promising star in tumor therapy. Immune checkpoint blockade (ICBs) targeting ICPs molecules such as ipimab, pabolizumab and atzumab have shown encouraging effects in many solid tumors. However, the e cacy of PD-1/ PD-L1 based immunotherapy for PDAC still needs further exploration, and it is urgent to exploit new immune checkpoints for PDAC.
RNA modi cation confers upon tumor cells the ability to rapidly and inversely alter after transcription to allow them to adapt and survive in the rapidly changing tumor microenvironment (TME) [6] . Of more than 170 RNA modi cations, N6-methyladenosine (m6A) modi cations are the most common and most extensively studied post-transcriptional modi cations of RNA [6] . Long non coding RNAs (lncRNAs) are non-coding RNAs greater than 200 nucleotides in length and are vital for pathogenesis of cancers [7] .
LncRNAs play essential roles at the transcriptional / post transcriptional level and chromatin modi cation by regulating gene expression [8] .M6A can regulate lncRNA structure through binding sites for m6A readers, this may allow speci c RNA binding proteins to access m6A residues. Similarly, m6A can also affect the function of lncRNAs through the ceRNA network. For example, FAM225a, a lncRNA modi ed by m6A regulation, acts as a sponge for mir-590-3p and mir-1275 in nasopharyngeal carcinoma [9] . A study concluded that m6A reader IGF2BP2 can promote stemness-like properties of PDAC through regulates lncRNA DANCR [10] . However, the full function of m6A and its regulated lncRNAs in PDAC still needs further elaboration.
The PDAC tumor microenvironment contains a large proportion of dense desmoplasia and immunosuppressive cell populations, which limit cytotoxic T responses [11] . Cancer associated broblasts (CAFs) contribute to the broblast proliferative and generate immunosuppressive TME by producing extracellular matrix proteins and cytokines, as well as promote tumor cell growth [12] . PD-L1 on the surface of PDAC parenchymal tumor cells and PD-1 of T cells can be important mechanism for immunotherapy. However, ICBs targeting PD-L1 have limited effects in PDAC due to an attenuated tumor immunogenic and immunosuppressive TME [13] .
Yuan and Hu [14,15] previously constructed a prognostic signaling pathway based on m6A regulated lncRNAs in pancreatic ductal adenocarcinoma, but in Yuan's study, the shortcoming is lack of detailed exploration of immune microenvironment in PDAC and further exploration of immune in ltration. Hu's study started from m6A regulated lncRNAs shared by TCGA and ICGC, PDAC samples were not clustered and potential therapeutic targets were not further explored. In addition, some samples were additionally removed, which has the potential to increase bias. In our study, we promulgate the association of m6Aregulated lncRNAs with tumor prognosis, TME, immune in ltration, and expression of related ICPs genes in PDAC patients. we de ned a 11 m6A regulated lncRNAs that effectively predicted the prognosis of patients with PDAC.

Results
Expression of m6A RNA methylation regulators and m6A-regulated lncRNAs in PDAC The m6A RNA methylation regulators play an indispensable role in tumor occurrence, proliferation, invasion and metastasis. We probe and systematically download the 23 m6A regulator genes expression and lncRNA expression of 178 tumor samples and 4 adjacent paracancerous samples from TCGA database and 167 normal samples from GTEx database. All the expressions of 23 m6A RNA methylation regulator genes differed between tumor and normal samples (P<0.01).
Correlation between consensus clustering analysis of m6A-regulated lncRNAs and survival of PDAC patients Concordance clustering analysis was performed on PDAC samples, which were identi ed in the cumulative distribution function (CDF), k takes values ranging from 2-9 and k denotes cluster counting ( Figure S2). According to similarity of m6A-regulated lncRNAs expression and ambiguous proportion for clustering measurements, the optimal clustering parameter was k = 2 ( Figure 3A) We merged the survival information and target m6A-regulated lncRNAs expression data of PDAC patients. A total of 177 PDAC patients were grouped into Cluster 1 (n=87) and Cluster 2 (n=90) based on the m6A-regulated lncRNAs expression. Comparing the clinical variables and demographic variables between Cluster 1 and Cluster 2 according to 28 m6A-regulated lncRNAs, it was found that there was no difference, but there was a signi cant difference in the overall survival between the two groups. The overall survival of Cluster 1 was superior to that of Cluster 2 (P<0.001) ( Figure 3B).
Correlation between PD-L1 and other immune checkpoints gene expression and m6A-related lncRNAs in Cluster 1 and Cluster 2 To explore the relationship between m6A-regulated lncRNAs and immune checkpoints (ICPs) gene expression, we compared the distribution of ICPs gene between tumor group and paracancerous cohort and between Cluster 1 and Cluster 2( Figure 7A-H). There was no difference in PD-L1 gene expression between Cluster 1 and Cluster 2 as well as tumor tissues and normal tissues. (Figure 7A, B). However, the expression of the CD47 (P<0.05) ( Figure 7C), CD276 (P<0.001) ( Figure 7E), IDO-1 (P<0.001) ( Figure 7G) genes that potential be ICPs was statistically different between the tumor and normal groups, showing signi cantly upregulated expression in the tumor tissue compared to the normal tissue. Meanwhile, the CD276 (B7-H3) gene also showed expression discrepancy in Cluster 1 and Cluster 2 (P<0.001) ( Figure   7F). Co-expression analysis of the genes showed a positive correlation between PD-L1 expression and m6A-regulated lncRNA U62317.1 and LINC01705 ( Figure 7I).

Strati cation analysis of the m6A-related lncRNAs in immune in ltration and tumor microenvironment
The tumor microenvironment (TME) is complex and constantly evolving. In addition to stromal cells, broblasts, and endothelial cells, the TME contains immune cells [24] . Stromal cells in the TME may secrete growth factors, which may not only promote tumor cell growth, but may also chemoattract other cells to migrate into the TME. Immune cells like natural killer and T cells, which can eliminate tumor cells, portend a promising outcome. Researchers have endeavored to understand the characteristics of the

TME and the status of immune cells and their interactions with tumor cells. Certain antibody therapies
have demonstrated having the ability to activate a patient's own immune system to kill tumor cells [25] . We further evaluate the association between m6A-regulated lncRNAs and tumor immune in ltration.
Stromal and in ltration level of 22 immune cells was evaluated by calculating the stromal, immune and ESTIMATE score using CIBERSORT. In the subtyping of PDAC samples based on 28-m6A regulated lncRNAs, it can be seen that B cells naive, B cells memory, Plasma cells, CD8 T cells, T cells CD4 memory resting, Monocytes, Macrophages M0, and Dendritic cells activated in Cluster 1 and Cluster 2 have differential distributions ( Figure 8A). In Cluster 1, B cells naive, plasma cells, CD8 T cells and Monocytes were in ltrated at higher levels in the TME, whereas in Cluster 2, expression levels of B cells memory, T cells CD4 memory resting and Macrophages M0 were upregulated ( Figure 8A). The Cluster 1 had higher stromal score, immune score and estimate score, implicated lower tumor purity ( Figure 8B, C, D). All of these scores have statistically signi cant (P<0.001).

Gene Set Enrichment Analyses (GSEA)
The Gene Set Enrichment Analyses (GSEA) was performed to explore the signature of m6A regulated lncRNAs in Cluster 1 and Cluster 2. Filtering criteria were FDR < 0.05 to nd notable enriched pathways. As shown in Figure 9, the obvious enriched signature pathway in Cluster 2 was "cell cycle" in Gene Set Enrichment Anaylses (GSEA). Unfortunately, no signi cant signature enrichment was found in Cluster 1.

Construction of m6A-regulated lncRNAs prognostic nomogram
The 177 PDAC samples with survival time and status in TCGA were divided into a train set (N=108) and a test set (N=69) according to 6:4 proportion. Based on the 28 m6A-regulated lncRNAs as well as the LASSO regression analysis ( Figure S3), we nally obtained 11 m6A-regulatedlated lncRNAs correlated with prognosis in train set. The risk score for each patient in train and sets was calculated as follows: risk score = U62317.1 * 0.340993840 + LINC00941 * 0.120489300 + AL139287.1 * (-0.173219714) + AC012615.1 * (-0.075084474) + CASC8 * 0.229254017 + Z97832.2 * (-0.279993778) + AC099778.1 * (-0.037908532) + AC245041.2 * 0.058969635 + AC091057.1 * 0.580745126 + TSPOAP1-AS1 * (-1.587783887) + AC005332.3 * (-0.324442618). After that, the respective patients were divided into high and low risk groups based on the median risk score of the training set in both the training and testing sets. Multivariate Cox analysis in the train set showed that age, tumor grade, and the risk score were independent risk factors associated with prognosis ( Figure 4A). Multivariate Cox analysis in the test set also indicated that the risk score was an independent risk factor associated with prognosis ( Figure 4B).
The distributions of the risk score and survival status of training set and testing set are presented in Figure 5(training set: A, C, E and G; testing set: B, D, F and H). The Kaplan-Meier curve of OS and log-rank t test was obtained in training set and testing set, patient's survival probabilities in the two datasets were statistically differed (P<0.001 in training set and P=0.008 in testing set) ( Figure 5E, F). The ROC predicting 5-year survival probabilities in training set and testing set were presented in Figure 5G and Figure 5H. The AUC for 1-year, 3-year and 5-year OS in training set was 0.861, 0.829 and 0.885 and in testing set was 0.708, 0.785 and 0.818. Therefore, we constructed a prognostic risk nomogram based on the risk score ( Figure 6).

The clinicopathological factors, clusters and immune-scores associated with risk scores in PDAC
Based on the median score of 11 prognostic m6A-regulated lncRNAs, 177 PDAC patients were divided into low-risk and high-risk groups. We analyzed the clinicopathological variables, cluster analysis, gene mutations and immune scores in the low-risk and high-risk groups in both the training and testing sets. The visualized clinical relevant heatmap based on 11-m6A regulated lncRNAs between low-risk and highrisk group is shown in Figure 10. With a difference in distribution between the low -and high-risk groups were the immune scores (P<0.05) and cluster analysis (P<0.001)( Figure 10A), which were higher in the low-risk group, and cluster analysis, which showed a denser distribution of Cluster 2 in the high-risk group. Patients with a low immune score, have an elevated risk and it follows that the immune score is favorable for prognosis. In terms of the distribution of clinicopathological characteristics and risk scores, patients with T3-T4 stage had an elevated risk compared with T1-T2 (P=0.016) ( Figure 10B). PD-L1 expression was upregulated in the high-risk compared with the low-risk group (P=0.037).
We parallel the Kaplan-Meier curve of OS among different age, gender, clinicopathological grade, TNM stage. As shown in Figure11, the OS curves showed distribution discrepancy in terms of all clinical variables and demographic characteristics except M1 and stage III-IV variables. Patients in the low-risk group of these variables all had a superior overall survival than those in the high-risk group ( Figure 11A-F) (P<0.01).

Association of m6A-related lncRNAs with immunocytes
We performed correlation analysis between the risk score and the 22 immune cells in ltrated in the TME, deleting samples with P < 0.05 for immune cells in TME. A total of 8 immune cells associated with the risk score were obtained, which were B cells naive, Plasma cells, CD8 T cells, T cells CD4 memory resting, Monocytes, Macrophages M0, Macrophages M1, and Dendritic cells activated. Positively correlated with the risk score were Dendritic cells activated, Macrophages M0, Macrophages M1, CD4 T cells ( Figure 12A-D). However, B cells naive, plasma cells, CD8 T cells, Monocytes were negatively correlated with the risk score ( Figure 12E-H).

Discussion
Pancreatic ductal adenocarcinoma remains one of the most fatal malignancies, and the morbidity and mortality rates are increasing worldwide. Advances in surgical technology have provided opportunities for patients with localized resectable PDAC, while patients with metastatic urgently need more effective comprehensive therapy [26] . More distinctive, PDAC has a highly desmoplastic and immunosuppressive TME [27] . With the development of high-throughput sequencing, breakthroughs in the treatment of PDAC have been pursued at the molecular and cellular level, such as to explore novel molecular therapeutic targets [15] . Studies have found that the role of aberrant lncRNAs in the carcinogenicity of PDAC cannot be ignored [28,29] . M6A modi cation is the most common epigenetic methylation modi cation of lncRNAs [30] . Previous studies have shown that m6A methylation modi ed lncRNAs can signi cantly affect the function of targets through RNA-protein interactions in a variety of tumor genes [31,32] . However, the function and mechanism of m6A-regulated lncRNAs in PDAC remain to be explored.
In this study, by performing gene differentially expression analysis in PDAC tumor tissues and normal tissues and co-expression analysis between 23 m6A RNA methylation regulators and lncRNAs by combined TCGA and GTEx databases, we identi ed 28 m6Aregulated lncRNAs correlated with prognosis. Based on this 28 m6A regulated lncRNAs, PDAC samples were categorized into Cluster 1 and Cluster 2 by consensus clustering analysis. In Cluster 1, the KM curves showed superior overall survival than Cluster 2. However, in Cluster 2, CD276 expression was higher than in cluster 1. After nally screening 11 optimal prognosic correlated m6A regulated lncRNAs by LASSO regression analysis of 28 m6A-regulated lncRNAs, PDAC patients were divided into high-risk and low-risk groups according to the risk scores. Univariate and multivariate cox regression analysis of risk score, TNM stage, pathological grade, age, gender indicated that only risk scores were independent risk factors associated with prognosis. A prognostic risk nomogram was constructed based on the risk scores of 11 m6A-regulated lncRNAs, and the AUC curve indicated the preferable accuracy and usefulness. Yuan and Hu constructed m6A-modi ed lncRNAs prognostic signature associated with overall survival, respectively [14,15] . Yuan's study enrolled only ve m6A-modi ed lncRNAs, and Hu proposed a prognostic pathway based on four m6A-related lncRNAs. In our study, 11 m6A regulated lncRNAs that were the most relevant with prognosis were included according to LASSO regression analysis, based on which the prognostic model was constructed, which can reduce bias and exhibit the strongest correlation with prognosis.
Immunotherapy has triggered a shift in treatment patterns in many solid tumors [33] . Due to the low mutational burden of PDAC tumor genes as well as the dense, inaccessible TME characteristic for brotic, hypoxic, and immunosuppressive properties renders it notoriously resistant to immunotherapy [34] . Anti-PD-1 and anti-PD-L1 therapy has shown promising results in many solid tumors, but its e cacy in PDAC is limited [35,36] . A clinical trial study showing a 9.1 months median progression free survival (mPFS) and 15 months median OS in metastatic PDAC patients receiving concurrent anti-PD-1 therapy with gemcitabine / nab paclitaxel chemotherapy, which are the relatively desirable achieved by current anti-PD-1 therapies in PDAC [37] . In the present study, although patients in the high-risk group had a higher level of PD-L1 gene expression than those in the low-risk group (P=0.037), PD-L1 expression in PDAC tumor tissues and normal tissues did not differ. Therefore, it is urgent to exploit new immune checkpoints for PDAC to develop sensitive ICIs or combination chemotherapy to achieve satisfactory treatment e cacy. The outcome that CD47, CD276 and IDO1 genes were upregulated in PDAC tumor tissues compared with paracancerous tissues and had distribution differences (P<0.05). Akane [38] concluded that colorectal cancer (CRC) patients with upregulated CD47 expression have a poor prognosis and that antibodies against the CD47 immune checkpoint are effective for CRC. B7 homologous protein 3 (B7-H3, also known as CD276), a newly identi ed immunoregulatory protein member of the B7 family, which is also a novel target for tumor immunotherapy. CD276 is overexpressed in tumor tissues while limited expression is observed in normal tissues, and in TME plays an integral role [39] . Studies have con rmed CD276was upregulated in breast cancer, non-small cell lung and ovarian cancer tumor tissues [39] . Indoleamine 2,3 dioxygenase 1 (IDO-1), encoded by the IDO1 gene on human chromosome 8p12, is the rate limiting enzyme in the conversion of tryptophan to kynurenine. Its expression is low in normal tissues and upregulated in a variety of tumors [40] . Overexpressed IDO1 is under interferon regulation inhibits apoptosis and promotes dormancy of tumor repopulating cells (TRCs) [41][42][43] . In our study, CD47, CD276 and IDO-1 were upregulated in PDAC tumor tissues compared with normal pancreatic tissues. It is therefore reasonable to suggest that the CD47, CD276 and IDO1 genes serve as potential immune checkpoints in PDAC. Targeting them with immune checkpoint inhibitors or combined chemoradiotherapy might bring a new dawn to PDAC therapy.
PDAC samples were classi ed into cluster 1 and cluster 2 according to m6A regulated lncRNAs, for which GSEA analysis showed that the tumor-related signaling pathway was cell cycle in Cluster 2, which indicated a higher risk of tumor progression in Cluster 2. Meanwhile, Cluster 2 had a worse prognosis compared to Cluster 1. This suggests a potential correlation between m6A regulated lncRNAs and tumor progression.
Based on 11 prognosis correlated m6A regulated lncRNAs, PDAC patients were divided into high-risk and low-risk groups, and patients in the high-risk group had a signi cantly worse long-term overall survival than those in the low-risk group. PD-L1 was signi cantly upregulated in the high-risk group. However, it was not possible to distinguish PDAC tumor tissues from normal tissues on PD-L1 expression. Therefore, surgery, chemoradiotherapy combined with potential ICIs should be sought to synergistically against PDAC with a promising outcome to achieve.

Conclusion
We exhaustively evaluated the m6A-regulated lncRNAs in relation to the PDAC tumor microenvironment, speci c genetic mutations, differences in immune checkpoint expression, and prognosis. Therefore, highrisk PDAC patients can be distinguished, and furthermore, potential novel immune checkpoints may guide clinical drug development and immunotherapy.

Data source
We download the RNA transcription pro les and clinical information of PDAC patients from The Cancer Genome Atlas (TCGA) of UCSC Xena (https://xena.ucsc.edu/). Because the sample of paracancerous tissues in TCGA database was limited, the Genotype-Tissue Expression project (GTEx) database as a supplement to TCGA PDAC database. We downloaded the 167 gene pro les of PDAC normal samples from website https://xena.ucsc.edu/, whose gene expression was corrected and stored in FPKM form, which is consistent with the gene expression pro le form in TCGA and processed in expression form. R package "sva" 1 was used to remove batch effects between TCGA and GTEx. The "sva" 1 package supports the use of the "sva" function for the estimation of proxy variables and the "Combat" and "fsva" functions for the adjustment of batch and latent variables from known batch effects and prediction problems [16] . The RNA matrix le was divided into mRNA le and lncRNA data according to human genome annotation data. The correlation between 23 m6A regulator genes and lncRNAs was determined by the Pearson correlation coe cient. Absolute correlation coe cient > 0.4 and P value < 0.05 was considered as m6A related lncRNAs.
Construction of a prognostic risk model and evaluation the e cacy Least absolute shrinkage and selection operator (LASSO) regression analysis was selected to construct the prognosis risk signature of m6A regulated lncRNAs. Then we eliminate collinearity of the m6A regulated lncRNAs and avoid over-tting of the constructed model. Subsequently, multivariate Cox regression analysis was utilized to construct m6A regulated lncRNAs prognosis signature. The risk score was calculated using the following formula: Riskscore= Where cod is the coe cient and Xi is the transformed relative expression value of each selected lncRNA.Using this formula, each patient's risk score can be calculated, then the PDAC cohort was categorized into low-risk and high-risk groups based on median risk scores of training set. Through univariate and multivariate Cox regression analysis in training set, the clinical variables and m6A regulated lncRNAs associated with prognosis screened to construct the nomogram. The prognostic effect was assessed by area under the time-dependent ROC (t-ROC) analysis in the training, and testing set. The clinicopathological features were compared between low-risk and high-risk group and visualized in heatmap.

Gene expression difference analysis and co-expression of PD-L1 and other ICs
We select PD-L1 and several other potential immune checkpoint genes such as CD47, CD276 and IDO1 [19][20][21] to visualize the distribution between normal and tumor, Cluster 1 and Cluster 2, low-risk and highrisk groups. The association between PD-L1 expression and m6A regulated lncRNAs was analyzed and visualized with R package "corrplot" 5 .

Comparation of immune in ltration and Gene Set Enrichment Analysis (GSEA)
Stromal, immune and estimate scores were calculated by using the ESTIMATE algorithm and the "estimate" 6 R package [22] . The differences between Cluster 1 and Cluster 2 were analysed. CIBERSORT deconvolution algorithm was used to obtain the abundance of tumor in ltrating 22 immune cells [23] . We use 1,000 permutations algorithm to calculate the nal scores, the distribution of the immune in ltration of 22 immune cell subtypes between Cluster 1 and Cluster 2 groups was compared. The Gene Set Enrichment Analyses (GSEA) 4.1.0 was utilized to elaborate downstream Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in Cluster 1 and Cluster 2 group based on m6A-regulated lncRNAs subtype data. The c2.cp.kegg.v7.4.symbols.gmt was selected as the reference gene set database. The criteria for meeting the conditions are as follows: P <0.05 and FDR<0.1. 1 https://bioconductor.org/packages/release/bioc/html/sva.html 2 https://cran.rproject.org/package=igraphdata Statistical analysis All statistical analysis is performed on R (version 4.0.4) software (https://www.r-project.org/) and multivariate Cox regression analysis was performed to determine the m6A-related lncRNAs associated with survival. The Kaplan-Meier curves were implemented to compare the different OS between the highrisk and low-risk groups according to distinct clinicopathological features. The Chi-square test or Fisher's exact test was used to analyse differences clinical information. The m6A lncRNA subtypes, clinicopathological features, risk scores, immune checkpoint expression, and immune in ltration were analyzed by a Pearson correlation test. A P value <0.05 was considered statistically signi cant. The ow chart can be seen in Figure 1. The ow chart of this study. Atlas (TCGA) and GTEx (Genotype-Tissue Expression) datasets. P < 0.05 ("*"), P < 0.01 ("**"), and P < 0.001 ("***").    The risk prognosis nomogram constructed by the risk score to predict the 1 -, 3 -, and 5-year overall survival (OS) probabilities of PDAC patients  Gene set enrichment analysis (GSEA) and genes mutation analysis were performed to predict the potential functions and pathways between the two clusters.