Transcriptomic insights into adenoid cystic carcinoma via RNA sequencing

Background: The aim of this study was to investigate the underlying mechanisms of adenoid cystic carcinoma (ACC) at the transcriptome level. Materials and methods: We obtained paired tumor and normal salivary gland tissues from 15 ACC patients, which were prepared for RNA sequencing. Results: Gene enrichment analysis revealed that the upregulated pathways were mainly involved in axonogenesis, and the downregulated pathways were mainly related to leukocyte migration, the adaptive immune response, lymphocyte-mediated immunity, and the humoral immune response. T-cells, B-cells and NK cells showed low infiltration in ACC tissues. In addition to the gene fusions MYB-NFIB and MYBL1-NFIB, a new gene fusion, TVP23C-CDRT4, was also detected in 3 ACC tissues. PRAME was significantly upregulated in ACC tissues, while antigen-presenting human leukocyte antigen (HLA) genes were downregulated. Conclusion: We found a new gene fusion, TVP23C-CDRT4, that was highly expressed in ACC. PRAME may be an attractive target for ACC immunotherapy.


Introduction
Adenoid cystic carcinoma (ACC) is a rare and aggressive malignancy, representing 1% of head and neck cancers (Sung et al., 2021). ACC is characterized by slow local growth, a high incidence of perineural invasion (PNI), infrequent regional metastases, frequent development of local recurrences and mostly slowly progressive and relatively indolent distant metastases (Szanto et al., 1984;Tang et al., 2010). One study suggested that the overall 5-year survival rate was 68%-90%, while the 10-year and 15-year survival rates were 52% and 28%, respectively, mainly due to PNI, local recurrence and distant metastasis (DM) (Atallah et al., 2020). DM is a common malignant feature of ACC, with an incidence rate of 16.1%-72.7%. The most common site is the lung. The overall 5-year, 10-year and 20-year survival rates of patients without DM were 85.6%, 67.4%, and 50.4%, and those of patients with DM were 69.1%, 45.7%, and 14.3%, respectively (Gao et al., 2013;Jang et al., 2017;Fang et al., 2022). The incidence of PNI in ACC was 48%-82%, and PNI was an independent prognostic factor negatively correlated with prognosis. The 5-year, 10-year and 15-year disease-specific survival (DSS) rates of patients with ACC of the head and neck and with PNI were 62%, 53%, and 27%, respectively, while the DSS rates of patients without PNI were as high as 90% (Amit et al., 2015;de Morais et al., 2021). The mainstay treatment for ACC is surgical resection in combination with adjuvant radiotherapy (Sood et al., 2016). For patients with recurrent, unresectable, or metastatic ACC, chemotherapy has limited efficacy (Papaspyrou et al., 2011).
Since its first appearance in (Margulies et al., 2005, highthroughput sequencing has become a useful approach for understanding biological processes at the molecular level and conducting detailed research to elucidate the genome and transcriptome. As an essential part of high-throughput sequencing, RNA sequencing (RNA-seq) has become a powerful tool for comprehensive characterization of the whole transcriptome at the gene and exon levels, with a unique ability to identify gene alterations, gene fusions, novel splicing variants, and transcripts with high resolution and efficiency (Morozova and Marra, 2008;Wang et al., 2009;Meyerson et al., 2010;Ren et al., 2012). RNA-seq studies have made great contributions in various fields, especially in the field of cancer research, including studies on differential gene expression analysis, cancer biomarkers, cancer heterogeneity and evolution, cancer drug resistance, the cancer microenvironment and immunotherapy, and neoantigens (Hong et al., 2020). Ferraroto and colleagues revealed two molecular subtypes: ACC-I (37%) and ACC-II (63%). MYC signaling drives ACC-I via NOTCH mutations or direct amplification, and the upregulation of TP63 and receptor tyrosine kinase represents a less aggressive clinical course of ACC-II (Ferrarotto et al., 2021). Through transcriptome sequencing, Ross and co-workers found that compared with the more aggressive tumor type salivary adenocarcinoma, ACC had a significantly lower TP53 mutation rate and tumor mutation burden (Ross et al., 2017).
However, because of the low incidence rate of ACC, only a limited number of scholars have studied ACC based on omics data. In this study, we explored 15 ACC tissues and their matched normal salivary gland (NSG) tissues through RNA-seq, attempting to discover novel markers in ACC.

Materials and methods
We declare that all methods were carried out in accordance with the World Medical Association Declaration of Helsinki. The study was approved by the Ethics Committee and was conducted under the guidance of international ethical standards (IRB number: PKUSSIRB-202165081). Written informed consent documents were obtained from all of the patients in our study.

Sample collection
A total of 15 pairs of ACC tumor and normal salivary gland (NSG) tissues were collected from January 2020 to March 2021. The inclusion criteria were as follows: 1. The patients with primary tumors had not received preoperative chemotherapy, radiotherapy, or any other anticancer treatment prior to surgery. 2. The patients received surgical treatment for primary tumors in our hospital. 3. The tumors were confirmed as ACC by two pathologists from the Department of Diagnostic Pathology in our hospital. The exclusion criteria were as follows: 1. Patients who declined to participate in the study; 2. Recurrent or metastatic ACC samples. NSG tissues were collected from the submandibular glands of patients with primary ACC undergoing neck dissection. Tissue specimens were freshly frozen immediately after surgery and stored at −80°C.

RNA extraction and qualification
RNA was extracted from tissues using TRIzol (Invitrogen, California, United States) according to the manufacturer's instructions. Fresh samples were ground to fine powder in liquid nitrogen, transferred to a tube, mixed with 1 mL of TRIzol, and then kept at room temperature for 10 min. Then, 200 μL of chloroform was added, and the sample was centrifuged (12,000 rpm, 10 min) at 4°C. An equal volume of phenol: chloroform (v:v 25:24) was added to the supernatant and centrifuged for 10 min at 12,000 rpm and 4°C. The supernatant was added to an equal volume of chloroform and centrifuged for 10 min at 12,000 rpm and 4°C. An equal volume of isopropyl alcohol was added to the supernatant; the sample was kept at −20°C for 1 h, and then centrifuged for 10 min at 12,000 rpm at 4°C; the supernatant was discarded. The pellet was washed with 1 mL of 75% ethanol and centrifuged (X 8000 rpm, 5 min) at 4°C, and the supernatant was discarded. The residual ethanol was removed by pipetting and dried under a vacuum. Finally, 50 μl of RNase-free water was added to dissolve the RNA (Wei et al., 2022). RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using a NanoPhotometer ® spectrophotometer (IMPLEN, CA, United States). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States).

Library preparation for transcriptome sequencing
A total of 1 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using the NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, United States) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized using random hexamer primers and Frontiers in Genetics frontiersin.org 02 M-MuLV Reverse Transcriptase (RNase H). Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of the 3′ ends of DNA fragments, NEBNext adaptors with hairpin loop structures were ligated in preparation for hybridization (Zhu et al., 2018). To preferentially select cDNA fragments 250-300 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, United States). Then, 3 µL USER Enzyme (NEB, United States) was incubated with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. Then, PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and Index (X) Primer (Zhu et al., 2018). Finally, PCR products were purified (AMPure XP system), and library quality was assessed on the Agilent Bioanalyzer 2100 system (Wei et al., 2022).

Clustering and sequencing
Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina NovaSeq platform, and 150 bp paired-end reads were generated (Wei et al., 2022).

Analysis of RNA-seq data
Quality control analyses of raw reads were first performed by FastQC and MultiQC, followed by read trimming by FastP. Filtered data were then mapped to GENCODE v37. The reads of the transcripts were further aligned and quantified against GENCODE human transcriptome v37 using Salmon (v1.4.0), imported and summarized in R by tximeta. The median number of reads per sample was 20.4 million. DESeq2 Love et al., 2014) was then used to perform exploratory analysis and differential expression analysis. PCA suggested that two normal tissue samples were grouped with tumor tissues, suggesting that they were outliers of the normal tissues. These two normal tissue samples were excluded from the following analysis. p values calculated by Wald statistics were corrected by Benjamini-Hochberg methods. Adjusted p values lower than 0.05 were considered significant. Functional enrichment of differentially expressed genes (DEGs) using gene sets from Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) and MSigDB hallmark gene sets was performed with clusterProfiler. For gene fusion analysis, prefiltered reads were first aligned against the genome reference from CTAT_GENOME_LIB (GRCh38_ gencode_v37) by STAR. Candidate fusion transcripts were identified by STAR-Fusion and further visualized by visualization tools provided by Arriba (Haas et al., 2019;Uhrig et al., 2021). Identification and analysis of isoform switches were performed by IsoformSwitchAnalyzeR using default settings (Vitting-Seerup and Sandelin, 2019).

Immunohistochemistry
The tissue samples for immunohistochemistry were collected from the tumors and the submandibular glands of the patients with primary ACC of head and neck undergoing neck dissection. Twelve pairs of paraffin-embedded ACC and NSG tissues were sectioned coronally to generate 3 μm slices. The sections were incubated with the following primary antibodies as indicated: CD3 (Gene Tech, GA045207) at a 1:20 dilution ratio; CD4 (Gene Tech, GT219107) at a 1:20 dilution ratio; CD8 (Gene Tech, GT211207) at a 1:20 dilution ratio and CD19 (Gene Tech, GT212802) at a 1:50 dilution ratio. The samples were stained with horseradish peroxidase-labeled secondary antibody (Bond Polymer Refine Detection, Leica; RS9800) and hematoxylin. The slides were photographed using a Leica microscope. IHC images were analyzed by 2 pathologists. They randomly selected five images of the peritumoral and intratumoral regions of ACC and NSG under ×400 magnification. Then, Image-Pro Plus 6.0 (Media Cybernetics, Rockville, MD, United States) was used to automatically count the number of positive cells and manually remove false-positive cells. The cell density = the total number of positive cells in 5 pictures/0.125 mm 2 .

Clinicopathological characteristics
To gain insight into the molecular pathogenesis of ACC, we collected 15 ACC tissues and 15 paired NSG tissues. The clinicopathological information of patients with ACC is shown in Table 1, including age, sex, tumor location, pathological type, PNI status, and TNM classification. Figure 1A illustrates the histopathological features of ACC. The tumor tissue presents cribriform or tubular changes.

Exploratory analysis and DEG analysis of ACC samples
After RNA-seq and preliminary filtration, a total of 15 ACC and 13 NSG samples were subjected to downstream analysis. NSG tissues from two cases failed the quality inspection. DESeq2 Love et al., 2014) was then used to perform exploratory analysis and differential expression analysis. PCA and sample distances were used to assess the differences between ACC tissues and NSG tissues. PCA plots showed that the first principle component (PC1) accounted for 60% of the overall variance of the data, and the second principle component (PC2) accounted for 12% ( Figure 1B), suggesting an obvious difference between the two groups. The heatmap of intersample correlation is shown in Figure 1C. The results showed that tumor and normal tissues were clustered separately, indicating distinct changes in the transcriptome in tumor samples.

Functional annotation and pathway enrichment analysis
GO and KEGG pathway enrichment analyses were performed to further understand the biological significance of the DEGs. The significantly enriched GO biological process (BP) terms of the upregulated pathways in the tumor samples were primarily related to extracellular structure organization, skeletal system development, nuclear division, and axonogenesis ( Figure 2A), indicative of changes in cell proliferation and the tumor microenvironment. The significantly enriched GO BP terms of the downregulated pathways were mainly involved in leukocyte migration, the adaptive immune response, lymphocyte-mediated immunity, and the humoral immune response, suggesting a lack of immune response in local ACC tissues ( Figure 2B).
Gene set enrichment analysis (GSEA) of KEGG pathways, including both activated and suppressed pathways, was also performed ( Figure 2C). In our analysis, activated pathways in the ACC groups were mainly related to the terms EMC-receptor interaction, microRNAs in cancer, focal adhesion, and neuroactive ligand-receptor interaction, consistent with the GO BP analysis. The suppressed pathways were mainly related to oxidative phosphorylation, salivary secretion, pancreatic secretion, and several metabolic pathways, suggesting loss of normal salivary Further signaling pathway analysis by ACC subtype revealed that some classical pathways, the VEGF, MAPK, PI3K, TGFb, WNT, estrogen, and p53 pathways, were activated ( Figure 2H).
Furthermore, as seen from the signaling pathway results, most of these pathways were inactivated in solid-type ACC, suggesting that therapy in this subtype may be different from that in other pathological types (cribriform-tubular or cribriform).

Screening of hub genes suggests dynamic changes in the immune microenvironment
According to the GO BP enrichment analysis, genes involved in axonogenesis, synapse organization, extracellular structure organization, and mesenchyme development were enriched, as shown in Figures  Some investigators have demonstrated that ACC exhibits low infiltration of antitumor immune cells (being a cold tumor). In our study, deconvolution of the cellular composition data was also performed. The R package MCP-counter (18) was applied to the normalized log2-transformed FPKM expression matrix to produce the absolute abundance scores for eight major immune cell types (CD3 + T-cells, CD8 + T-cells, cytotoxic lymphocytes, natural killer (NK) cells, B lymphocytes, monocytic lineage cells, myeloid dendritic cells and neutrophils), endothelial cells and fibroblasts. The deconvolution profiles were then hierarchically clustered and compared across the control and ACC groups ( Figure 4A).
The downregulated genes involved in leukocyte migration and the adaptive immune response are shown in Figures  Regarding M1/M2 macrophages, most M1 macrophage-related genes were downregulated, whereas 4 M2 macrophage-related genes were upregulated and 1 was downregulated, suggesting that overall inflammation was downregulated ( Figures 5A-K).  To confirm the changes in immune infiltration in ACC tissues, we further performed immunohistochemistry staining for CD3, CD4, CD8, and CD19 in the ACC tissues. We used the staining intensity of cells to assess the immune infiltration level of these markers and analyzed their differences by one-way ANOVA-Tukey's multiple comparisons test. The Supplementary Materials show these statistical data. We first compared the expression differences of these markers in the normal and peritumoral and intratumoral regions ( Figure 6A). Compared with that in normal tissues, the expression of CD3, CD4, CD8, and CD19 in ACC tissues decreased, and the expression in the intratumoral region decreased more significantly than that in the peritumoral region. In addition, the expression levels of these markers were obviously different among the three pathological types of ACC ( Figure 6B). Lower expression levels of CD3, CD4, and CD8 were observed in cribriform ACCs, while the expression level of CD19 was the lowest in solid ACCs. Finally, there was a significant difference in the expression level of these markers in different tumor invasion types ( Figure 6C). The lowest expression level of CD3 was observed in the bone-invasive type, and CD4 and CD8 were observed in the non-invasive type. The expression level of CD19 was generally low in all types. Consistent with the above RNA-seq results, the immunohistochemistry results showed decreased infiltration of T-cells and B-cells in the tumor microenvironment compared to the normal tissues ( Figure 6D).   Frontiers in Genetics frontiersin.org 07

Detection of gene fusions
In addition to DEG analysis, we also performed analysis of gene fusions. In our 15 ACC samples, a total of 29 possible gene fusions were detected ( Figure 7A). The breakpoints of different types of gene fusions are shown in Table 3. Among them, 11 NFIB-related gene fusions were detected in 15 samples ( Figure 7B), and there were different subtypes of fusions (Table 4). The gene fusions involving the transcription factors MYB and NFIB consisted of different regions in chromosomes 6 and 9, as shown in Figures 7C-E, while parts of chromosomes 8 and 9 fused to form the fusion gene of the transcription factors MYBL1 and NFIB ( Figures 7F-H). The fusion of genes may change gene functions via different regulatory mechanisms. The gene expression of NFIB and MYB was significantly upregulated in most samples, while that of MYBL1 was also upregulated in 5 samples ( Figures 7I-K). Moreover, the gene fusion TVP23C-CDRT4 was also detected in three ACC tissues, and different types of fusions were detected, as shown in Table 4. TVP23C-CDRT4 gene fusions include different regions of chromosome 17, and the fusion structures are shown in Figures 7L, M. In addition, the expression of TVP23C and CDRT4 was significantly upregulated in 2 ACC samples, and the pathological subtype of these 2 samples was solid type (Figures 7N, O).

Expression of cancer germline antigens
Cancer germline antigens play an important role in the occurrence and progression of tumors. To find novel markers, we further screened cancer germline antigens. The top 35 antigens expressed differently in the ACC tissues vs. NSG tissues were detected. Among them, PRAME was the most significantly upregulated in all ACC samples. The results of the analysis of cancer germline antigen expression in different patients are shown in Figure 8A.

Different isoform analysis of mRNAs
In addition to differential expression of genes, mRNA isoform switching may also lead to functional changes in genes, even without significant changes in total gene expression. Therefore, we also performed isoform switching analysis in our samples. We identified 1,690 isoformEs for all significant genes: 1,318 pairs of isoforms had significantly changed usage in opposite directions, and 1,162 genes had significantly differential transcript usage. Among them, 1,336 isoforms were identified with potential functional consequences: 1,063 pairs of isoforms had significantly changed usage in opposite directions, and 922 genes had significantly differential transcript usage ( Figures 9A, B).
Genome-wide analysis of isoform switches revealed that, compared to normal tissues, ACCs had a significantly higher fraction of switches resulting in complete open reading frame (ORF) loss, domain loss, intrinsically disordered region (IDR) loss, and shortened ORFs (Figures 9C-E). In the genome-wide analysis of alternative splicing and changes, ACCs had a significantly higher fraction of switches, resulting in loss of alternative 3′ acceptor sites (A3s), alternative transcription termination sites (ATTSs), and alternative 5ʹ acceptor sites (A5s) and gain of multiple exon skipping (MES) and alternative transcription start sites (ATSSs) (Figures 9F-H). The difference between the isoforms involved in an isoform switch may arise from changes in ATSSs, alternative splicing, and ATTSs. The relative importance of alternative splicing, ATSSs, and ATTSs is shown in Figure 9I. The top genes with isoform switching, including LMNB2, RHBDL1, TRAF3IP3, VGLL4, and NFIB, are shown in Figures 9J-N.

Discussion
ACC has been traditionally subdivided into 3 histological groups (cribriform, tubular, and solid) based on the cellular components of the tumor (Perez et al., 2006). Because of its rarity and distinctive clinical features, further exploration of the underlying mechanisms in ACC is imperative. With the rapid development of nextgeneration high-throughput sequencing technology (Shendure, 2008), RNA-seq has become a new and important means of gene expression and transcriptome analysis, and it can quickly and comprehensively reveal the sequence information of almost all RNA transcripts in tissues. In our study, we assessed multiple transcriptomic factors of ACC, including DEGs, gene fusions, cancer germline antigens, and different isoforms of mRNA, by analyzing 15 ACC tissues and 13 paired normal tissues.
In our study, the top 10 upregulated DEGs were COL27A1, APBA2, TUBA1A, TBX1, FNDC1, PLSCR3, FABP7, TP53, ABCC1, and EN1, while the top 10 downregulated DEGs were CST2, SMR3B, KLK1, CTBS, DHRS2, STATH, SH3BGRL2, HTN1, SLC13A5, and SLC9A1. Research has shown that EN1 and FNDC1 are significantly  highly methylated genes in ACC (Bell et al., 2011). In addition, EN1 has been confirmed to be involved in cell-substrate adhesion and neural guidance and is a potential biomarker for a poor prognosis in ACC (Bell et al., 2012;Lin et al., 2022). Supplementary Table S1 shows the relevant functions of these top 20 upregulated and downregulated DE mRNAs in detail. Studies have confirmed that VEGF could promote the progression of tumors by inducing the formation of tumor vessels, and the high expression of VEGF was related to a poor prognosis in ACC patients . The MAPK/ERK and Wnt/β-catenin signaling pathways have been proven to be related to the progression and metastasis of ACC, and blocking these pathways regulates the migration and invasion of ACC cells. MicroRNA-181a suppresses salivary adenoid cystic carcinoma metastasis by targeting the MAPK-Snai2 pathway (Ji et al., 2020). In ACC, the mutation rate of PIK3CA is 5%-22%, and this mutation or activation of the PI3K/AKT pathway by upstream signals (EGFR, HIF-1 α, etc.) could promote the invasion, progression and metastasis of ACC cells (Sun et al., 2021;Han et al., 2022). The mutation rate of p53 in ACC was 8%, and its expression was strongly correlated with solid ACC subtypes. The mutant p53 protein not only lost its original function of inhibiting cancer but also induced excessive cell proliferation, leading to cellular malignant transformation and increased tumor invasion. Compared to ACC patients without p53 mutations, those with p53 mutations had a higher probability of metastasis, local recurrence and nerve infiltration and a shorter median survival period (Adderley et al., 2021). Functional annotation analysis showed that the upregulated DEGs were significantly enriched in the GO BP term axonogenesis. The genes involved in axonogenesis are mainly those encoding neurotrophins and their receptors, including NGF, BDNF, NTF4, and NTRK3. In previous studies, NGF,    Frontiers in Genetics frontiersin.org BDNF, and NTRK3 were found to be overexpressed in ACC tissues and were significantly correlated with invasion, metastasis and a poor prognosis. NGF is a member of the NGF-beta family and encodes a secreted protein that stimulates nerve growth activity and is involved in the regulation of growth and the differentiation of sympathetic and certain sensory neurons. NGF was significantly overexpressed in ACC tissues and was more strongly stained in tumors with PNI than in those without PNI. In addition, high NGF expression leads to a worse prognosis (Hao et al., 2010). BDNF encodes a protein belonging to the nerve growth factor family. BDNF is highly expressed in ACC, and the BDNF/TkB axis promotes PNI and a poor prognosis by inducing EMT progression (Jia et al., 2015). NTRK3 encodes TrkC, a member of the neurotrophic tyrosine receptor kinase family, and blocking the autocoline NT-3/TrkC signaling pathway can suppress ACC nerve invasion and growth (Ivanov et al., 2013). Neurotrophin-4 (NTF4) belongs to the family of neurotrophic factors (NTFs), which are most commonly known for their roles in the nervous system (Shen et al., 2010). BDNF and NT-4/5 can stimulate resistance to apoptosis in breast cancer cells through p75NTR and TrkB-T1 and then promote tumor progression (Vanhecke et al., 2011). NTF4 promotes cell invasion and a poor prognosis in colorectal cancer through autophagy regulation (Yang et al., 2020). However, the function of NTF4 in ACC is unclear. The importance of NTF4 in ACC needs to be further explored. Furthermore, leukocyte migration has been found to be downregulated in ACCs, and the downregulation of leukocyte migration has been found to be significantly correlated with the downregulation of chemokines, including CCL5, CCL19, CCL20, CCL21, CXCL6, and CXCL12. It has been reported that the CCL5/CCR5 axis could promote the migration, invasion and PNI of ACC . CXCL12/ CXCR4 participate in the course of EMT and Schwann-like differentiation in PNI and can promote ACC cell invasion and PNI through the Twist/S100A4 axis . In addition, CCL20 enhances the anticancer response of the immune system by recruiting dendritic cells (Bonnotte et al., 2004). The increased expression of inflammatory CCL19 and CCL21 leads to the infiltration of tumor infiltrating lymphocytes (TILs) and improves the anticancer effect and prognosis of multiple types of cancer, including liver cancer, breast cancer, and colorectal cancer (Sharma et al., 2000;Liang et al., 2007;Wu et al., 2011). In summary, these downregulated chemokines are related to PNI, tumor progression and a type IV (immune tolerance) immune microenvironment in ACC. In addition, the downregulated DEGs in ACCs were mainly enriched in BP terms related to the adaptive immune response, lymphocyte-mediated immunity, and the humoral immune response. The downregulated DEGs associated with the humoral immune response mainly included B-cell-related and immunoglobulin genes. The downregulated DEGs associated with the adaptive immune response mainly included T-cell-related genes, especially CD8 + T-cell-related genes. These results suggest that lymphocytes play a crucial role in ACC. In addition, in our analysis, eight major immune cell types (CD3 + T-cells, CD8 + T-cells, cytotoxic lymphocytes, NK cells, B lymphocytes, monocytic lineage cells, myeloid dendritic cells and neutrophils) mostly showed low infiltration, and M1/ M2 macrophage analysis suggested downregulation of overall inflammation in ACC tissue. A microenvironment with low CD8 + tumor-infiltrating lymphocyte (TIL) levels was considered a type IV (immune tolerant) microenvironment, which was consistent with a previous report (Mosconi et al., 2019). The immune environment in ACC tissue is an important factor determining cell proliferation and the efficacy of immunotherapy, chemotherapy and radiotherapy (Demaria et al., 2015;Galluzzi et al., 2015;Weichselbaum et al., 2017). ACC recurrence was found to be associated with low CD8 + TIL counts (Mosconi et al., 2019). In our immunohistochemistry analysis, samples from patients with ACC were analyzed for the expression of CD3, CD4, CD8, and CD19, and the results suggested that the ACC microenvironment exhibits low immunogenicity, as represented by low T and B-cells. The type IV microenvironment may help the tumor escape from the immune system and partly explain the poor prognosis of ACC patients.
In ACC, the recurrent t (6; 9) (q22-23; p23-24) (Nordkvist et al., 1994) translocation in adenoid cystic carcinoma results in a novel fusion of the MYB proto-oncogene with the transcription factor gene NFIB (Persson et al., 2009). In contrast, the MYB-NFIB fusion is not expressed in non-adenoid cystic carcinoma neoplasms of the head and neck, confirming the high specificity of the MYB-NFIB fusion, which has been used as a diagnostic marker for ACC (Toper and Sarioglu, 2021). Moreover, in 2016, a translocation involving MYBL1 and NFIB on chromosomes 8 and 9 [t (8; 9)] was found in most ACC patients without MYB-NFIB fusion (Mitani et al., 2016). The activation of MYB is considered a potential therapeutic target for ACC. Tumor progression can be prevented by silencing MYB, targeting the downstream effector of this gene, and blocking the protein-protein interaction in the transcriptional complex (Liu et al., 2018). In our study, the MYB-NFIB fusion gene appeared in 8 of 15 (53.3%) ACC tissues, and the MYBL1-NFIB fusion gene appeared in 3 of 15 (20%) ACC tissues; the rates were slightly higher than those in the literature [45% (9/ 20) and 15% (3/20), respectively, as determined by RNA-seq using formalin-fixed, paraffin-embedded (FFPE) tissues] (Shibata et al., 2021). It is worth mentioning that a new gene fusion, the TVP23C-CDRT4 fusion, was found in our transcriptome sequencing study. TVP23C-CDRT4 gene fusions include different anatomical regions of chromosome 17 and were found in 3 of 15 (20%) ACC tissues. TVP23C is involved in protein secretion and vesicle-mediated transport. For patients with colorectal cancer, a high level of TVP23C in plasma indicates a longer survival period (Holm et al., 2019). CDRT4 is upregulated in breast cancer cells (Yellapu et al., 2022). However, research on TVP23C-CDRT4 gene fusion is very rare. This study is the first report of the fusion gene in ACCs. In addition, the expression of TVP23C and CDRT4 is upregulated mainly in solid-type ACCs (2/3), which indicates that TVP23C and CDRT4 may be useful for distinguishing solid-type ACC from the cribriform and tubular types of ACC. The TVP23C-CDRT4 fusion gene may promote the development of malignancy, and our results provide new targets for personalized treatment, although they need to be confirmed by further research.
Cancer cells can evade immune surveillance, which allows them to grow and metastasize. One mechanism of immune escape during cancer development is the downregulation of HLA class I molecule expression (Marincola et al., 2000;Seliger et al., 2002;Raffaghello et al., 2005;Lopez-Albaitero et al., 2006). HLA class I molecules play a central role in cell-mediated immunity, especially as antigen-presenting molecules for cytotoxic T lymphocytes (CTLs), which recognize tumor antigen-bound peptides presented on the cell surface through HLA class I molecules and kill the target cancer cell (Khong and Restifo, Frontiers in Genetics frontiersin.org 2002; Bubenik, 2004). HLA class I expression appears to be lost or downregulated on the tumor cell surface, which could represent a mechanism by which neoplastic cells escape killing by CTLs, allowing tumor dissemination and metastasis (Pistillo et al., 2000). HLA class I antigen downregulation has been used as a prognostic biomarker in esophageal cancer (Hosch et al., 1997), melanoma (Kageshita et al., 1999), lung cancer (Kikuchi et al., 2007), ovarian cancer (Vitale et al., 2005), renal cell carcinoma (Kitamura et al., 2007), bladder cancer (Homma et al., 2009), and oral squamous cell carcinoma (Koike et al., 2020) and has improved the survival of patients. In our study, the expression levels of class I HLA molecules, which are responsible for antigen presentation to tumor-killing T-cells, and the expression levels of class II HLA molecules, which are responsible for antigen presentation to helper T-cells, were decreased in some patients. HLA class I molecules play a central role in anticancer immunity, but their value in ACC remains unclear. Our data suggest that tumor cells may escape tumor immunity by downregulating HLA genes. \PRAME, first described as a melanoma antigen (Ikeda et al., 1997), is known to repress retinoic acid receptor signaling and to repress the transcription of genes involved in growth arrest, differentiation, and apoptosis in melanoma models (Epping et al., 2005). Subsequently, numerous studies reported that overexpression of PRAME was significantly correlated with clinicopathological features in malignant cancers, including medulloblastoma (Orlando et al., 2018), lung adenocarcinoma (Pan et al., 2017), uveal melanoma (Gezgin et al., 2017), high-grade serous cancer (Zhang et al., 2016), myxoid liposarcoma (Iura et al., 2015), osteosarcoma (Tan et al., 2012), bladder cancer (Dyrskjot et al., 2012), breast carcinoma (Epping et al., 2008), and solitary fibrous tumors . In our research, PRAME was obviously upregulated in ACC tissues. However, the role of PRAME in the occurrence and development of ACC remains unknown. This is the first time that PRAME has been found to be upregulated in ACC. It is common knowledge that the efficacy of chemotherapy and targeted drugs is poor in recurrent, unresectable, or metastatic ACC (Papaspyrou et al., 2011;Hanna et al., 2020). In addition, ACC tends to exhibit low levels of TILs (Mosconi et al., 2019). These features suggest that ACC may not be an ideal candidate for a lymphocyte immune checkpoint inhibitor-based approach to immunotherapy. Therefore, PRAME may be an attractive target for immunotherapy against ACC.
Alternative splicing is a ubiquitous regulatory mechanism of gene expression that allows the generation of more than one unique mRNA species from a single gene (Baralle and Giudice, 2017). Alternative splicing provides transcriptional plasticity by controlling which RNA isoforms are expressed at a given time in a given cell type. Cancer cells subvert this process to produce isoforms that benefit cell proliferation or migration or enable escape from cell death (Urbanski et al., 2018). Compared to normal tissues, ACCs have a significantly higher fraction of switches resulting in loss of complete ORFs, domains, and IDRs; intron retention; and shortened ORFs and have a significantly higher fraction of switches resulting in loss of A3s, ATTSs, and A5s and gain of MES and ATSSs. The differences between the isoforms involved in an isoform switch can arise via changes in three distinct biological factors: ATSSs, alternative splicing, and ATTSs. The top genes with isoform switching were LMNB2, RHBDL1, TRAF3IP3, VGLL4, and NFIB. Among them, LMNB2 is significantly upregulated in ACCs, and it is probably upregulated by mRNA isoforms with coding functions. VGLL4 is a transcription factor and is upregulated in ACCs. It has three isoforms with coding functions and a noncoding isoform. One of the isoforms with coding function is specifically upregulated in ACCs, and one domain of this isoform is different from that in other subtypes, which may have an impact on its transcriptional regulation function.
In our study, neurotrophins and their receptors were found to be closely related to the characteristics of PNI. The downregulated genes were mainly related to the antitumor immune response in the tumor microenvironment. Eight major immune cell types, CD3 + T-cells, CD8 + T-cells, cytotoxic lymphocytes, NK cells, B lymphocytes, monocytic lineage cells, myeloid dendritic cells and neutrophils, mostly showed low infiltration, and M1/M2 macrophage analysis suggested downregulation of overall inflammation in ACC tissue. Furthermore, in our preliminary study, T-cells and B-cells also showed low infiltration in ACC tissues by immunohistochemistry. Moreover, we found a new fusion gene, TVP23C-CDRT4, that may be a specific marker for solidtype ACC. We also assessed the expression of cancer germline antigens. PRAME was significantly upregulated in ACC tissues and may be an attractive target for immunotherapy against ACC. We also performed isoform and DEG analyses.

Conclusion
Through RNA-seq of 15 ACC tissues and 13 paired normal tissues, we found significantly upregulated and downregulated genes and pathways in ACC, and T-cells, B-cells and NK cells showed low infiltration in ACC tissues. In addition to the classic gene fusions MYB-NFIB and MYBL1-NFIB, a new gene fusion, TVP23C-CDRT4, was also detected in 3 ACC tissues, and the expression of TVP23C and CDRT4 was significantly upregulated in solid-type ACC. Screening of cancer germline antigens showed that PRAME was significantly upregulated, while antigen-presenting human leukocyte antigen (HLA) genes were downregulated.

Future perspectives
We found a new fusion gene, TVP23C-CDRT4, that may be a specific marker for solid-type ACC. PRAME is significantly upregulated in ACC tissues and may be an attractive target for immunotherapy against ACC. This study provides an accurate transcriptome analysis of ACC and a basis for the study of the molecular mechanism of ACC.

Data availability statement
Readers can obtain the data through NCBI database (BioProject accession number PRJNA900701).

Ethics statement
The studies involving human participants were reviewed and approved by Medical Ethics Committee of Peking University School Frontiers in Genetics frontiersin.org and Hospital of Stomatology. The patients/participants provided their written informed consent to participate in this study.

Author contributions
W-JW and JZ contributed to the concept and designed the research. Quality control of data and algorithms was performed by W-JW and P-GA. Data analysis and statistical analysis were performed by Y-FT, P-GA, B-XG, SY, and XH. P-GA and Y-FT commented on the manuscript preparation and editing. W-JW and JZ performed the manuscript review critically and approved the final version.

Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's note All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.