Veterinary and Comparative Oncology
Transcriptome sequencing reveals two subtypes of cortisol-secreting adrenocortical tumours in dogs and identifies CYP26B1 as a potential new therapeutic target
Karin Sanders 1,2
İD
| Hans S. Kooistra 1
İD Marieke van den Heuvel1 | |
Michal Mokry 3,4
Guy C. M. Grinwis 5
İD
| Noortje A. M. van den Dungen 4,6
Frank G. van Steenbeek 1,6
İD
|
Sara Galac 1
İD
1Department of Clinical Sciences, Faculty of Veterinary Medicine, Utrecht University, Utrecht, The Netherlands
2Princess Máxima Center for Pediatric Oncology, Utrecht, The Netherlands
3Central Diagnostics Laboratory, University Medical Center Utrecht, Utrecht, The Netherlands
4Laboratory of Experimental Cardiology, Department of Cardiology, University Medical Center Utrecht, Utrecht, The Netherlands
5Department of Biomolecular Health Sciences, Faculty of Veterinary Medicine, Utrecht University, Utrecht, The Netherlands
‘Department of Cardiology, Division of Heart and Lungs, University Medical Center Utrecht, Utrecht, Netherlands
Correspondence
Sara Galac, Department of Clinical Sciences, Faculty of Veterinary Medicine, Utrecht University, Utrecht, The Netherlands. Email: s.galac@uu.nl
Funding information
Nederlands Kankerfonds voor Dieren
Abstract
Cushing’s syndrome (CS) is a serious endocrine disorder that is relatively common in dogs, but rare in humans. In ~15%-20% of cases, CS is caused by a cortisol-secreting adrenocortical tumour (csACT). To identify differentially expressed genes that can improve prognostic predictions after surgery and represent novel treatment targets, we performed RNA sequencing on csACTs (n = 48) and normal adrenal cortices (NACs; n = 10) of dogs. A gene was declared differentially expressed when the adjusted p-value was <. 05 and the log2 fold change was >2 or < - 2. Between NACs and csACTs, 98 genes were differentially expressed. Based on the principal compo- nent analysis (PCA) the csACTs were separated in two groups, of which Group 1 had significantly better survival after adrenalectomy (p = . 002) than Group 2. Between csACT Group G1 and Group 2, 77 genes were differentially expressed. One of these, cytochrome P450 26B1 (CYP26B1), was significantly associated with survival in both our canine csACTs and in a publicly available data set of 33 human cortisol-secreting adrenocortical carcinomas. In the validation cohort, CYP26B1 was also expressed sig- nificantly higher (p = . 012) in canine csACTs compared with NACs. In future studies it would be interesting to determine whether CYP26B1 inhibitors could inhibit csACT growth in both dogs and humans.
KEYWORDS
adrenocortical carcinoma, canine diseases, Cushing’s syndrome, RNA-seq
|
1 INTRODUCTION
Hypercortisolism, also known as Cushing’s syndrome (CS), is a serious endocrine disorder that results in significant morbidity and mortality if left untreated.1 In dogs, CS is a relatively common disorder with a
Frank G. van Steenbeek and Sara Galac should be considered joint senior author.
@ 2022 The Authors. Veterinary and Comparative Oncology published by John Wiley & Sons Ltd.
prevalence of 1 per 400.2 In humans, CS is a rare disorder with a prev- alence of 39-79 per million.3 In both humans and dogs, ~15%-20% of CS cases are caused by a cortisol-secreting adrenocortical tumour (csACT).4,5 If metastases are not detected, the treatment of choice for a csACT is surgical removal (adrenalectomy).6
A csACT can be classified as benign (adrenocortical adenoma; ACA) or malignant (adrenocortical carcinoma; ACC).7 Although the histopathological classification of ACA or ACC has been well- described in humans,7-9 this is less the case in dogs. For both humans and dogs, the histopathological parameters that are commonly used to make this classification are known to have high interobserver vari- ability.8,10,11 Moreover, using this traditional histopathological classifi- cation of ACA or ACC, several studies observed no significant differences in survival times of dogs after adrenalectomy.12,13 We therefore recently developed a new histopathological scoring system, the Utrecht score, to predict the prognosis of dogs with a csACT after adrenalectomy.14 The Utrecht score is based on histopathological parameters, including the Ki67 proliferation index, that had low intra- and interobserver variability and were associated with survival times after adrenalectomy,14 and is comparable to the Helsinki score that is used in human ACTs.8 In addition to histopathology, gene expression profiling has been shown to be related to patient outcome in many solid tumour types, 15 including ACTs in humans,16 and could therefore further improve prognostic predictions for dogs with csACTs.
In addition to improving prognostic predictions, gene expression profiling could also provide more insight into potential therapeutic tar- gets. We previously identified three genes of which the expression was significantly associated with survival after adrenalectomy in dogs with csACTs: Steroidogenic factor-1, topoisomerase II alpha and pitui- tary tumour-transforming gene-1.17 However, in that study, we used a candidate gene approach with only 14 genes. We hypothesized that many more genes are associated with canine csACTs and survival after adrenalectomy. These genes or their products could represent potential novel therapeutic targets in dogs. Because dogs with sponta- neously occurring csACTs represent a potential animal model for humans, the identified genes would also be interesting to study in human csACTs. In this study, we performed RNA sequencing (RNA- seq) on 48 csACTs and 10 normal adrenal cortices (NACs) of dogs.
2 METHODS
|
|
2.1 Patients
For the discovery cohort (RNA-seq), 49 canine csACTs were collected between 2002 and 2017. For the validation cohort (RT-qPCR), 25 canine csACTs were collected between 2017 and 2020. Permis- sion to use the csACT tissues for research purposes was obtained from all dog owners. The suspicion of CS was based on clinical signs and routine laboratory findings that were consistent with hypercorti- solism. The presence of non-suppressible hypercortisolism was detected with the low-dose dexamethasone suppression test, or with urinary corticoid: creatinine ratios (UCCRs) in combination with a
high-dose dexamethasone suppression test.5 The presence of an adre- nal mass was visualized using abdominal ultrasonography, computed tomography, or both. All dogs underwent adrenalectomy, after which the tissues were formalin-fixed and paraffin-embedded for histopa- thology, and either snap-frozen or first stored in RNAlater™ Stabiliza- tion Solution (Invitrogen™, ThermoFisher Scientific). The tissues were stored at -70℃ until RNA isolation. Histopathology confirmed the adrenocortical origin of the adrenal masses in all cases. The Utrecht score was assessed as previously described (Ki67 proliferation index +4 if ≥33% of cells have clear/vacuolated cytoplasm +3 if necrosis is present).14 Of the included csACT samples, 36 (of which one was excluded from further analyses) were also included in our previously published studies on prognostic markers in canine csACTs.14,17 Adre- nal cortices of healthy dogs were used as reference material, these dogs were euthanized for reasons unrelated to this study, approved by the Ethical Committee of Utrecht University conform Dutch legis- lation. All methods were performed in accordance with the relevant guidelines and regulations.
2.2 Extraction, library preparation and RNA sequencing |
RNA was isolated from the csACTs or from the adrenal cortex of nor- mal adrenal glands using the RNeasy Mini Kit (Qiagen) according to the manufacturer’s instructions. RNA concentrations were measured with the Qubit™ RNA high sensitivity Assay Kit (ThermoFisher Scien- tific). The quality of RNA samples was assessed using the 2100 Bioa- nalyzer Instrument (Agilent Technologies) using RNA Nano Chips (Agilent Technologies).
CEL-Seq218 was used to generate sequencing libraries. The meth- odology captures 3’-end of polyadenylated RNA species and includes unique molecular identifiers (UMIs), which allow direct counting of unique RNA molecules in each sample; 20 ng of total RNA was precip- itated using isopropanol and washed with 75% ethanol. After remov- ing ethanol and air-drying the pellet, primer mix containing 5 ng primer per reaction was added, initiating primer annealing at 65℃ for 5 min. Subsequent RT reaction was performed; first strand reaction for 1 h at 42℃, heat-inactivated for 10 min at 70℃, second strand reaction for 2 h at 16℃, and then put on ice until proceeding to sam- ple pooling. The primer used for this initial reverse-transcription (RT) reaction was designed as follows: an anchored polyT, a unique 6 bp barcode, a unique molecular identifier (UMI) of 6 bp, the 5’ Illu- mina adapter and a T7 promoter. Each sample now contained its own unique barcode due to the primer used in the RNA amplification, mak- ing it possible to pool together up to 10 cDNA samples per pool. The cDNA was cleaned using AMPure XP beads (Beckman Coulter), washed with 80% ethanol, and resuspended in water before proceed- ing to the in vitro transcription (IVT) reaction (AM1334; Thermo- Fisher) incubated at 37℃ for 13 h. Next, primers were removed by treating with Exo-SAP (Affymetrix, Thermo-Fisher) and amplified RNA (aRNA) was fragmented and then cleaned with RNAClean XP (Beckman-Coulter), washed with 70% ethanol, air-dried, and
resuspended in water. After removing the beads using a magnetic stand, RNA yield and quality in the suspension were checked by Bioa- nalyzer (Agilent). The cDNA library construction was then initiated by performing an RT reaction using SuperScript II reverse transcriptase (Invitrogen/Thermo-Fisher) according to the manufacturer’s protocol, adding randomhexRT primer as a random primer. Next, PCR amplifica- tion was done with Phusion High-Fidelity PCR Master Mix with HF buffer (NEB) and a unique indexed RNA PCR primer (Illumina) per reaction, for a total of 11-15 cycles, depending on aRNA concentra- tion, with 30 s elongation time. PCR products were cleaned twice with AMPure XP beads (Beckman Coulter). Library cDNA yield and quality were checked by Qubit fluorometric quantification (Thermo- Fisher) and Bioanalyzer (Agilent), respectively.
2.3 Sequencing read mapping and quality filtering |
Libraries were sequenced on the Illumina Nextseq500 platform; a high output paired-end run of 2 x 75 bp was performed (Utrecht Sequenc- ing Facility). The reads were demultiplexed and aligned to canine cDNA reference build CanFam3.1 using the BWA (0.7.13) by calling ‘bwa aln’ with settings -B 6 -q 0 -n 0.00 -k 2 -1 200 -t 6 for R1 and -B 0 -q 0 -n 0.04 -k 2 -| 200 -t 6 for R2, ‘bwa sampe’ with settings -n 100 -N 100. Multiple reads mapping to the same gene with the same unique molecular identifier (UMI, 6 bp long) were counted as a single read. The raw read counts were corrected for UMI sampling (corrected_count = - 4096*(In(1- (raw_count/4096)))).
2.4 RNA-seq data analysis |
The raw data files are available in the Gene Expression Omnibus (GEO) repository, accession number GSE196108. TMM normalized counts and fold changes were calculated in edgeR. P-values were calculated with the Mann-Whitney U test, which were corrected for multiple comparisons using the FDR Benjamini-Hochberg method.
ToppFun was used for functional enrichment analysis based on functional annotations and protein interactions networks,19 including genes with a log2 fold change of >1.5 or < - 1.5 and FDR-corrected p- values of <. 05. Vulcano-plots were generated using EnhancedVulcano (R package version 1.12.0).20
2.5 RT-qPCR validation cohort |
After isolation of RNA, cDNA was synthesized with the iScript cDNA Synthesis Kit (Bio-Rad, Veenendaal, The Netherlands) according to the manufacturer’s instructions, and subsequently diluted to 1 ng/uL. RT-qPCR reactions were performed using SYBR-green Supermix (Bio- Rad) in a CFX384 Touch Real-Time PCR Detection System (Bio-Rad). Primers for CYP26B1 (forward: TCA CCT CTT CGA AGT CTA CC; reverse: AGT AGT CCT TGC CCT GG) were designed with PerlPrimer
software,21 checked for secondary structures with the Mfold web server,22 and ordered from Eurogentec (Maastricht, The Netherlands). The mRNA expression levels of four reference genes were analysed for data normalization: signal recognition particle receptor, succinate dehydrogenase complex subunit A, ribosomal protein S19, and tyro- sine 3-monooxygenase/tryptophan 5-monooxygenase activation pro- tein zeta.23,24 All reactions were performed in duplicate. The geNorm25 method was used to analyse the pairwise variance and sta- bility of reference gene expression, which justified their use. The rela- tive expression of CYP26B1 was calculated using the 2-44CT method.26
2.6 Statistical analyses |
In all dogs, assessment of UCCRs was routinely performed 2 months after surgery to assess remission. In the follow-up, endocrine testing (either UCCRs or low-dose dexamethasone suppression test) was performed when hypercortisolism-related clinical signs were present. Recurrence was confirmed when endocrine testing confirmed the clinical suspicion. In some dogs with recurrence, metastases or regrowth of the csACT were confirmed with diagnostic imaging and/or fine needle aspiration biopsy. The owners of all dogs that were still alive at the end of the study were asked to submit a urine sample of their dog for UCCR assessment. For the survival analyses, dogs were considered to have had an event when they were eutha- nized because of recurrence of hypercortisolism, resulting from either metastases or regrowth of the csACT. Dogs that died from an unrelated cause, were still alive at the end of the study, or were lost to follow-up, were censored at the timepoint they were last known to be alive. Survival times were recorded as the time between adre- nalectomy and censoring or death due to recurrence. Continuous variables were analysed using the Cox proportional hazards model, while bivariate variables were analysed using the Kaplan-Meier product-limit method with the log-rank test. To determine optimal cut-off values for group classification, receiver operating characteris- tic curves were used to find the value with highest Youden’s index (sensitivity + specificity -1).
For comparisons between groups, the Mann-Whitney U test was used. Survival analyses and Mann-Whitney U tests were performed using IBM SPSS Statistics (Version 25, IBM Corp.). p-values of <. 05 were considered significant.
3 RESULTS
|
|
3.1 ACTs versus NACs
For the initial analyses, we compared gene expression profiles between NACs (n = 10) and csACTs (n = 49). The clinical data of all included dogs can be found in Data S1. The RNA isolated from all samples was analysed with a Bioanalyzer, which gave a mean RNA integrity number (RIN) value of 8.7 (SD±0.8). The number of
(A) 30
(C)
. Normal adrenal cortex
· Adrenocortical tumor
20
PC2: 15% variance
10
0
-10
-20
-30
-20
-10
0
10
20
30
PC1: 20% variance
(B)
NS
. Log_ FC >2 or ← 2
Adj_P < 0.05
. Adj_P < 0.05 and Log2 FC >2 or ← 2
6
VAT1L
EFHD1
CNKSR2
AQP6
CHAC1
TDRD9
SQLE
S100A12
CHGB BCHE
CCDC33
ILSS
-Log P
ALDH1A1
BMP1
LON2
IRGS2
FAM20A
CALYO i
GPR22
BHLHE40
CYP26B1
PPBP
BMP7
MARCO
SERPINA5
HBM
3
GAL
C29H80rf34
ATP10B CD5L
HK2
INPP5J
IGRIK2 ECRG4
AHSP
KLK2
0
-5
0
5
Log2 fold change
| GO:Biological | Process Name | pValue | FDR B&H |
|---|---|---|---|
| GO:1901617 | organic hydroxy compound biosynthetic process | 1,56E-10 | 5,82E-07 |
| GO:0046165 | alcohol biosynthetic process | 1,46E-09 | 2,73E-06 |
| GO:1901615 | organic hydroxy compound metabolic process | 9,35E-09 | 1,06E-05 |
| GO:0006695 | cholesterol biosynthetic process | 1,42E-08 | 1,06E-05 |
| GO:1902653 | secondary alcohol biosynthetic process | 1,42E-08 | 1,06E-05 |
| GO:0016126 | sterol biosynthetic process | 2,80E-08 | 1,75E-05 |
| GO:1902930 | regulation of alcohol biosynthetic process | 3,47E-08 | 1,86E-05 |
| GO:0031667 | response to nutrient levels | 2,32E-07 | 1,08E-04 |
| GO:0007267 | cell-cell signaling | 2,59E-07 | 1,08E-04 |
| GO:0006066 | alcohol metabolic process | 3,45E-07 | 1,29E-04 |
| GO:0016125 | sterol metabolic process | 4,37E-07 | 1,49E-04 |
| GO:0019218 | regulation of steroid metabolic process | 5,14E-07 | 1,60E-04 |
| GO:0009991 | response to extracellular stimulus | 5,87E-07 | 1,69E-04 |
| GO:0033273 | response to vitamin | 1,33E-06 | 3,52E-04 |
| GO:0008203 | cholesterol metabolic process | 1,54E-06 | 3,52E-04 |
| GO:0090181 | regulation of cholesterol metabolic process | 1,59E-06 | 3,52E-04 |
| GO:0015837 | amine transport | 1,70E-06 | 3,52E-04 |
| GO:0009713 | catechol-containing compound biosynthetic process | 1,79E-06 | 3,52E-04 |
| GO:0042423 | catecholamine biosynthetic process | 1,79E-06 | 3,52E-04 |
| GO:0015842 | aminergic neurotransmitter loading into synaptic vesicle | 2,28E-06 | 4,27E-04 |
| GO:0015850 | organic hydroxy compound transport | 2,43E-06 | 4,33E-04 |
| GO:0008202 | steroid metabolic process | 3,02E-06 | 4,92E-04 |
| GO:1902652 | secondary alcohol metabolic process | 3,02E-06 | 4,92E-04 |
| GO:0022008 | neurogenesis | 3,74E-06 | 5,83E-04 |
| GO:0009719 | response to endogenous stimulus | 3,96E-06 | 5,93E-04 |
FIGURE 1 Differences in transcriptomic profiles of canine cortisol-secreting adrenocortical tumours (csACTs; n = 48) compared with normal adrenal cortices (NACs; n = 10). Principal component analysis (A) shows that the majority of csACTs cluster apart from the NACs. The volcano plot (B) shows the differentially expressed genes with FDR-corrected p-values <. 05 and log2 fold changes of >2 (n = 50 genes) or < - 2 (n = 48 genes) in red. The enrichment analysis (C) based on Gene Ontology: Biological Processes shows the top 25 of most affected pathways
reads per sample after quality control ranged from 993 628 to 5 385 270 (mean = 2 078 649), on which a correction for library size was performed. The Trimmed Mean of M-values (TMM) nor- malized counts can be found in Data S2. One csACT was deter- mined an outlier based on the principal component analysis (PCA) plot and was therefore excluded from further analyses, leaving a total of 48 csACTs. In the PCA plot, the csACT cohort largely clus- tered apart from the NAC cohort (Figure 1A). After false discovery rate (FDR) correction of the obtained P-values, 1452 genes were differentially expressed between csACTs and NACs (Data S3). Considering not only p-values but also (stringent) effect size mea- sures, 98 of these differentially expressed genes had a log2 fold change of >2 or < - 2, of which 50 were upregulated in ACTs com- pared with NACs, and 48 were downregulated (Figure 1B;
Data S4). The top 25 pathways that were affected in csACTs com- pared with NACs are shown in Figure 1C.
3.2 csACT classification |
Of the 48 included csACTs, follow-up information was available for 35 dogs (Data S1), which were also included in our previous stud- ies.14,17 We divided these 35 csACTs in two groups: those with a his- topathological Utrecht score of <6 (n = 13; low risk of recurrence tumours) and those with a score of ≥6 (n = 22; moderate-to-high risk of recurrence tumours).17 As shown in the PCA plot in Figure 2A, the samples with Utrecht scores <6 showed a potential tendency to form a different cluster from those with scores of ≥6, but the groups
WILEY
Veterinary and Comparative Oncology
(A) 20
(B) 1.0
10
0.8
Utrecht score <6
PC2: 14% variance
Cumulative survival
P = 0.014
0
0.6
-10
0.4
-20
. Normal adrenal cortex
0.2
· Utrecht score <6
Utrecht score ≥6
. Utrecht score ≥6
-20
-10
0
10
20
0
20
40
60
(C)
PC1: 22% variance
(D)
Survival (months)
20
1.0
10
0.8
PC2: 14% variance
Cumulative survival
P = 0.002
0
0.6
-10
0.4
csACT group 1
-20
. Normal adrenal cortex
0.2
· csACT group 1
- csACT group 2
. csACT group 2
-20
-10
0
10
20
0
20
40
60
PC1: 22% variance
Survival (months)
(E) 30-
. Normal adrenal cortex
· csACT group 1
· csACT group 2
20
· csACT unknown
PC2: 15% variance
10
0
-10
-20
-30
-20
-10
0
10
20
30
PC1: 20% variance
(A) 10 -
· NS
. Log, FC >2 or ← 2
Adj_P < 0.05
AGRN
· Adj_P < 0.05 and Log2 FC >2 or ← 2
PLAAT1
C4BPA
INHA
IGNA14
PLA2G5
NRCAM
TNFAIP6
ATP13A5
SH3RF3 PIP5K1B
-Log P
LSAMP
GG
SPOCK1
CNTFR
DUSP26
ECRG4
5
GREB1
CCDC33
CYP26B1
GAD1 IFITM10
TSPAN5
FFAR4
RDH8
F6F1
PIWIL24CYB5A
RRDM16
SLO1A5
RTN1
ADRA2A
LRRC4C
ATP 10B
PGPR85
SLC11A14CDH7
PDZK1IP1
FAM189A2
SLC35F3 NPP5J
SERPINF1
LMO4
SHH
PERP
HBM
ADGRB1
CHI3L1
ATP2C2
C8A
EGFR
TNMD>SLC2A3 AHSPCOL20A1
PLAUR OS100A12
ITGAZ_SLC35F1
LTF TSPAN33ªAJAP1
0
-4
0
4
Log2 fold change
(B)
(C)
Upregulated
Tumor vs normal
Group 2 vs group 1
42
8
38
Downregulated
Tumor vs normal
Group 2 vs group 1
48
0
31
| GO:Biological | Process | Name | pValue | FDR B&H |
|---|---|---|---|---|
| GO:0007267 | cell-cell signaling | 3.57E-12 | 1.41E-08 | |
| GO:0046903 | secretion | 5.12E-10 | 1.01E-06 | |
| GO:0140352 | export from cell | 3.75E-09 | 4.95E-06 | |
| GO:0032940 | secretion by cell | 5.06E-09 | 5.01E-06 | |
| GO:0050804 | modulation of chemical synaptic transmission | 3.77E-08 | 2.56E-05 | |
| GO:0099177 | regulation of trans-synaptic signaling | 3.88E-08 | 2.56E-05 | |
| GO:0098916 | anterograde trans-synaptic signaling | 5.17E-08 | 2.56E-05 | |
| GO:0007268 | chemical synaptic transmission | 5.17E-08 | 2.56E-05 | |
| GO:0099537 | trans-synaptic signaling | 6.29E-08 | 2.77E-05 | |
| GO:0099536 | synaptic signaling | 9.67E-08 | 3.83E-05 | |
| GO:0022610 | biological adhesion | 1.50E-07 | 5.42E-05 | |
| GO:0007155 | cell adhesion | 4.16E-07 | 1.37E-04 | |
| GO:0010038 | response to metal ion | 4.91E-07 | 1.50E-04 | |
| GO:0045055 | regulated exocytosis | 9.44E-07 | 2.67E-04 | |
| GO:0023061 | signal release | 2.01E-06 | 5.30E-04 | |
| GO:0000902 | cell morphogenesis | 3.65E-06 | 9.03E-04 | |
| GO:0016050 | vesicle organization | 4.40E-06 | 1.02E-03 | |
| GO:0051240 | positive regulation of multicellular organismal process | 6.14E-06 | 1.35E-03 | |
| GO:0006906 | vesicle fusion | 7.25E-06 | 1.51E-03 | |
| GO:0090174 | organelle membrane fusion | 8.00E-06 | 1.54E-03 | |
| GO:0006869 | lipid transport | 8.85E-06 | 1.54E-03 | |
| GO:0006887 | exocytosis | 9.36E-06 | 1.54E-03 | |
| GO:0140029 | exocytic process | 9.36E-06 | 1.54E-03 | |
| GO:0099500 | vesicle fusion to plasma membrane | 9.36E-06 | 1.54E-03 | |
| GO:0022600 | digestive system process | 9.69E-06 | 1.54E-03 |
FIGURE 3 Differences in transcriptomic profile between csACT Group 1 (n = 22) and csACT Group 2 (n = 26). The volcano plot (A) shows the differentially expressed genes with FDR-corrected p-values <. 05 and log2 fold changes of >2 (n = 46 genes) or < - 2 (n = 31 genes) in red. The enrichment analysis (B) based on Gene Ontology: Biological Processes shows the top 25 of most affected pathways in csACT Group 1 compared with csACT Group 2. The pie-charts in (C) show the overlap in genes that were either up- or downregulated both in csACTs compared with NACs (left) as well as in csACT Group 2 compared with csACT Group 1 (right). csACT: cortisol-secreting adrenocortical tumour.
overlapped to a large extent. As expected, based on their inclusion in our previous study concerning the Utrecht score, these samples had significantly different (p = . 014) survival times after adrenalectomy when classified as Utrecht score of <6 or ≥6 (Figure 2B).
Because the classification based on Utrecht scores <6 or ≥6 did not result in clearly separated clusters, we divided these 35 samples into two groups based on their PCA profile: csACT group 1 (n = 17) and Group 2 (n = 18; Figure 2C). In survival analysis, these groups had
significantly different (p = . 002) survival times (Figure 2D), with a p- value that was lower than when separating the samples based on hav- ing an Utrecht score of <6 or ≥6.
We subsequently added the 13 additional samples of which we had little to no follow-up information on survival in the PCA. Based on their positions in the PCA plot, we allocated these samples to either csACT group 1 (n = 5, group total n = 22) or csACT Group 2 (n = 8, group total n = 26; Figure 2E).
1.0
0.8
CYP26B1 low
Cumulative survival
0.6
P = 0.013
0.4
0.2
- CYP26B1 high
0
24
48
72
96
120
144
Survival (months)
3.3 csACT Group 2 versus Group 1 |
When comparing csACT Group 2 with csACT Group 1, a total of 2277 genes were differentially expressed after FDR correction (Data S3). Of these, 77 had a log2 fold change of >2 (n = 46) or < - 2 (n = 31) in csACT Group 2 compared with csACT Group 1 (Figure 3A; Data S4). The top 25 pathways that were affected in csACT Group 2 compared with csACT Group 1 are shown in Figure 3B.
As a therapeutic target, a gene would be most interesting when it is up- or downregulated in csACTs compared with NACs, as well as in csACT Group 2 (poor prognosis) compared with csACT Group 1 (favourable prognosis). We found eight genes that were upregulated with log2 fold change >2 in csACTs compared with NACs, as well as in csACT Group 2 compared with csACT Group 1, but no genes that were downregulated with log2 fold change ← 2 in both comparisons (Figure 3C). Of these eight upregulated genes, three were significantly associated with survival after surgery in our cohort of 35 patients with follow-up data: CD5 antigen-like (CD5L; p = . 026, hazard ratio [HR] = 1.030), Cytochrome P450 26B1 (CYP26B1; p =. 031, HR = 1.010), and Oesophageal Cancer-Related Gene 4 (ECRG4; p = . 013, HR = 1.046).
3.4 Comparison with human ACCs |
To determine whether CD5L, CYP26B1 and ECRG4 could also be rele- vant in human ACTs, we consulted a RNA-seq data set of The Cancer Genome Atlas (TCGA) including 79 human ACCs, as described by
Zheng et al. (2016)27 and uploaded in the web-based genomic analysis and visualization platform R2 (http://r2platform.com).28 To make the human data set more comparable to our canine data set, we selected only those ACCs with known excessive cortisol secretion (either alone or in combination with other hormones, n = 33). When dividing these patients in R2 based on high or low expression of the respective genes in the ACCs, calculated using their optimal cut-off values, there were no significant differences in overall survival for CD5L (p = . 174) or ECRG4 (p = . 091; Figure S1), but there was a significant difference (p = . 013) in survival of patients with high (n = 14) or low (n = 19) CYP26B1 expression (Figure 4).
Because of its potential as a therapeutic target based on its func- tion29 and expression profile in both human and canine csACTs, and the availability of CYP26B1 inhibitors,30 we decided to more closely study CYP26B1 expression.
3.5 Canine csACT validation cohort CYP26B1 |
When comparing samples with high CYP26B1 expression (n = 13) to those with low CYP26B1 expression (n = 23) in our cohort of 35 dogs with follow-up data, we found that those with high CYP26B1 expres- sion had a significantly worse prognosis than those with low expres- sion (p = . 004; Figure 5A). To validate our RNA-seq findings in the discovery cohort (10 NACs and 48 csACTs), we performed RT-qPCR for CYP26B1 mRNA expression on an independent patient cohort (8 NACs and 25 csACTs). As in the discovery cohort, CYP26B1 mRNA was expressed significantly higher in csACTs than in NACs in the vali- dation cohort (Figure 5B).
|
4 DISCUSSION
This is the first study in which RNA-seq was performed on canine csACTs to identify new prognostic markers and potential therapeutic targets for adrenal-dependent CS. We have found many genes and biological processes that were dysregulated in csACTs compared with NACs. Additionally, we identified two transcriptionally distinct groups of csACTs that had significantly different survival times after adrenal- ectomy and identified CYP26B1 as a potential new therapeutic target for both dogs and humans with csACTs.
Based on our previously developed histopathological Utrecht score,14 we initially classified the tumours based on having an Utrecht score of <6 or ≥6. Although this classification showed a potential ten- dency for separation of samples in the PCA plot, there was substantial overlap between the groups. We therefore classified the samples based on their position in the PCA plot as either csACT Group 1 or csACT Group 2 and found that this classification improved the signifi- cance of the survival analysis. This suggests that there is still room for improvement in the histopathological analysis of csACTs. In future studies it would be interesting to determine whether assessment of additional markers that are upregulated in csACT Group 2 compared with csACT Group 1, such as CYP26B1, CD5L and ECRG4, could
(A)
1.0
0.8
Cumulative survival
P = 0.004
0.6
CYP26B1 low
0.4
0.2
- CYP26B1 high
0
20
40
60
Survival (months)
(B)
Discovery cohort
Validation cohort
1000
Normal adrenal cortex
100
Normal adrenal cortex
Adrenocortical tumor
Adrenocortical tumor
P = 0.0003
P = 0.012
CYP26B1 expression RNA-seq Trimmed mean of M-values
CYP26B1 expression RT-qPCR Relative expression (fold change)
100
10
10
1
1
0
0
improve the prognostic prediction of dogs with csACTs after adrenal- ectomy. The survival analyses, however, do have some limitations. The follow-up monitoring after surgery was not standardized for all dogs, because some were monitored in our University clinic, while some were monitored in external specialist clinics. Although assess- ment of UCCRs was routinely performed 2 months after surgery for all dogs, the follow-up period after this was not standardized. There- fore, it is possible that cases with recurrence but with limited clinical signs of hypercortisolism could have been missed. In addition, for the
dogs with recurrence, it was not always determined whether this was caused by metastases or by regrowth of the tumour.
While csACT Group 1 has a more favourable prognosis compared with csACT Group 2, it does not cluster more closely with the NACs in the PCA plots. In human ACTs, transcriptomic analyses showed that ACAs cluster closely together with NACs, while ACCs cluster apart from the NAC/ACA groups.16 CsACT Group 1 and csACT Group 2 might therefore not represent the classical ACA versus ACC classifi- cation. The underlying cause for the different gene expression
patterns in these groups is currently still unknown. It would be inter- esting to perform, for example, whole genome sequencing on these samples to determine whether the different gene expression patterns and survival times could be explained by different tumour-driving mutations.
In the normal adrenal cortex, progenitor cells are located in the sub- capsular region within the zona glomerulosa.31 After differentiating into functional hormone-producing zona glomerulosa cells, these cells migrate inward over time, transdifferentiate into cells of the zona fasci- culata and subsequently of the zona reticularis, and will eventually undergo apoptosis at the cortical-medullary boundary.31 Of note is that several zona glomerulosa markers, such as SHH31 and DAB,32 are expressed significantly higher in csACT Group 2 than in csACT Group 1, while the opposite is true for several zona reticularis markers, such as CYB5A33 and FGG34 (Figure 2, Data S3 and S4). These different expres- sion profiles could suggest that the group classification is related to the location in the adrenal cortex from which the csACTs arise. For exam- ple, csACT Group 2 tumours could arise from progenitor cells located in the zona glomerulosa, while csACT Group 1 tumours could arise from further differentiated cells near or in the zona reticularis. Alternatively, the csACTs could either de-differentiate towards a more progenitor-like phenotype (csACT Group 2) or further differentiate towards their end- station (csACT Group 1). In both hypotheses, the progenitor-like pheno- type of csACT Group 2 compared with the further differentiated phe- notype of csACT Group 1 could explain the observed differences in survival times after adrenalectomy.
As a therapeutic target, a gene would be most interesting when it is expressed higher (or lower) in csACTs compared with NACs, as well as in csACTs with poor prognosis compared with those with favour- able prognosis. Theoretically, pharmaceutically targeting a product of such a gene could then work in (almost) all csACTs, but even better in those with poor prognosis (which is the most important group in terms of requiring additional therapeutics). We therefore looked at which genes were differentially expressed in both comparisons. We found eight genes that fulfilled this criterium, of which three were sig- nificantly associated with survival times after surgery. One of these genes, CYP26B1, was not only prognostically relevant in our data set of canine csACTs, but also in a publicly available data set of 33 human cortisol-secreting ACCs. Additionally, CYP26B1 expression was also increased in an independent validation cohort of 25 csACTs compared with 8 NACs of dogs, as assessed with RT-qPCR.
In addition to its prognostic value, the function of the CYP26B1 protein also makes it an interesting treatment target. CYP26B1 is an enzyme that metabolizes all-trans-retinoic acid (ATRA), the most biolog- ically active metabolite of vitamin A (retinol).29 ATRA is an important regulator of cell differentiation, proliferation and apoptosis35 and has been reported to inhibit cancer development by inducing differentiation and/or apoptosis.36 CYP26B1 can metabolize ATRA into less biologi- cally active intermediates,29 and therefore reduces ATRA-induced dif- ferentiation or apoptosis.35 Previous studies have shown that increased expression of CYP26A1, an enzyme that is closely related to and has similar functions as CYP26B1, increases the resistance of cells to apop- togenic factors.35 Increasing CYP26 expression could therefore be a protection mechanism of tumour cells against the influence of ATRA.
Indeed, studies in other tumour types such as neuroblastoma have shown that in response to ATRA, cells can increase their CYP26A1 expression.37 Subsequently, inhibition of CYP26B1 activity could be a valuable treatment strategy. CYP26B1 has also been identified as a prognostic marker in other cancer types, such as colorectal cancer.29
In human ACTs, a previous meta-analysis study identified dis- turbed retinoic acid (RA) signalling as a major pathogenetic path- way.38 To our knowledge, this has not previously been linked to CYP26B1. In a subsequent study, these researchers tried to exploit this finding by treating xenograft mice inoculated with an ACC cell line with 9-cis retinoic acid (9-cisRA), an isomer of ATRA. They found that 9-cisRA significantly reduced the tumours’ Ki67 prolifera- tion index.39 Although this was promising, the usefulness of ATRA therapy is limited because resistance can rapidly emerge, partly due to accelerated RA metabolism.40 A strategy to overcome this is to inhibit the enzymes that metabolize RA, also referred to as retinoic acid metabolism blocking agents (RAMBAs).40 Especially in tumour types where CYP26B1 seems to be involved in tumorigenesis, such as in canine and human csACTs, using RAMBAs specifically aimed at CYP26B1 seem particularly useful to inhibit tumour growth. CYP26B1 inhibitors are available and are under investigation for several diseases. 30,41,42
5 CONCLUSION |
With this study, we have gained important insights into the molecular pathways that underly tumorigenesis in canine csACTs. We have identified two transcriptionally distinct groups of csACTs. These groups do not seem to represent the classical ACA versus ACC classi- fication, but rather a more progenitor-like phenotype versus a further differentiated phenotype. It would be interesting to determine whether the classification of these groups is related to their muta- tional background. We have identified CYP26B1 as a novel prognostic marker and as a promising potential therapeutic target in both canine and human csACTs. This represents an important proof-of-concept of research in spontaneously-occurring ACTs of dogs that could lead to novel discoveries in ACTs of humans, highlighting the importance of this comparative model. In future in vitro and in vivo studies, it would be interesting to determine whether CYP26B1 inhibitors can inhibit tumour growth in both dogs and humans with ACTs.
ACKNOWLEDGEMENTS
We would like to thank the Dutch Cancer Fund for Animals (Nederlands Kankerfonds voor Dieren; NKFD) for their financial con- tribution. We would like to thank Adri Slob for technical support; and Drs. Jenny Buijtels, Sylvie Daminet, Federico Fracassi, Rob Gerritsen, Dieneke Jongepier, Bart Sjollema, Marjanne Zaal and Giora van Stra- ten for contributing samples to our csACT data set.
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.
DATA AVAILABILITY STATEMENT
The raw data files are available in the Gene Expression Omnibus (GEO) repository, accession number GSE196108. The analysed data can be found in the Data S1.
ORCID
Karin Sanders DD https://orcid.org/0000-0001-9634-9853
Hans S. Kooistra ID https://orcid.org/0000-0002-8058-0492
Guy C. M. Grinwis DD https://orcid.org/0000-0002-1583-9648
Frank G. van Steenbeek (D https://orcid.org/0000-0001-5460-6540 Sara Galac DD https://orcid.org/0000-0002-4831-4995
REFERENCES
1. Creemers SG, Hofland LJ, Lamberts SW, Feelders RA. Cushing’s syn- drome: an update on current pharmacotherapy and future directions. Expert Opin Pharmacother. 2015;16(12):1829-1844.
2. O’Neill DG, Scudder C, Faire JM, et al. Epidemiology of hyperadreno- corticism among 210,824 dogs attending primary-care veterinary practices in the UK from 2009 to 2014. J Small Anim Pract. 2016; 57(7):365-373.
3. Lacroix A, Feelders RA, Stratakis CA, Nieman LK. Cushing’s syndrome. Lancet. 2015;386(9996):913-927.
4. Sharma ST, Nieman LK, Feelders RA. Cushing’s syndrome: epidemiol- ogy and developments in disease management. Clin Epidemiol. 2015; 7:281-293.
5. Galac S, Reusch CE, Kooistra HS, Rijnberk A. Adrenals. In: Rijnberk A, Kooistra HS, eds. Clinical Endocrinology of Dogs and Cats. 2nd ed. Schlütersche; 2010:93-154.
6. Sanders K, Kooistra HS, Galac S. Treating canine Cushing’s syndrome: current options and future prospects. Vet J. 2018;241:42-51.
7. Papotti M, Libè R, Duregon E, Volante M, Bertherat J, Tissier F. The Weiss score and beyond-histopathology for adrenocortical carcinoma. Horm Cancer. 2011;2(6):333-340.
8. Pennanen M, Heiskanen I, Sane T, et al. Helsinki score - a novel model for prediction of metastases in adrenocortical carcinomas. Hum Pathol. 2015;46(3):404-410.
9. Renaudin K, Smati S, Wargny M, et al. Clinicopathological description of 43 oncocytic adrenocortical tumors: importance of Ki-67 in histoprognostic evaluation. Mod Pathol. 2018;31:1708- 1716.
10. Aubert S, Wacrenier A, Leroy X, et al. Weiss system revisited: a clini- copathologic and immunohistochemical study of 49 adrenocortical tumors. Am J Surg Pathol. 2002;26(12):1612-1619.
11. Labelle P, Kyles AE, Farver TB, de Cock HEV. Indicators of malignancy of canine adrenocortical tumors: histopathology and proliferation index. Vet Pathol. 2004;41(5):490-497.
12. Anderson CR, Birchard SJ, Powers BE, Belandria G, Kuntz CA, Withrow SJ. Surgical treatment of adrenocortical tumors: 21 cases (1990-1996). J Am Anim Hosp Assoc. 2001;37:93-97.
13. Schwartz P, Kovak JR, Koprowski A, Ludwig LL, Monette S, Bergman PJ. Evaluation of prognostic factors in the surgical treatment of adrenal gland tumors in dogs: 41 cases (1999-2005). J Am Vet Med Assoc. 2008;232(1):77-84.
14. Sanders K, Cirkel K, Grinwis GCM, et al. The Utrecht score: a novel histopathological scoring system to assess the prognosis of dogs with cortisol-secreting adrenocortical tumours. Vet Comp Oncol. 2019; 17(3):329-337.
15. Bucca G, Carruba G, Saetta A, Muti P, Castagnetta L, Smith CP. Gene expression profiling of human cancers. Ann N Y Acad Sci. 2004;1028:28-37.
16. Giordano TJ, Kuick R, Else T, et al. Molecular classification and prog- nostication of adrenocortical tumors by transcriptome profiling. Clin Cancer Res. 2009;15(2):668-676.
17. Sanders K, van Staalduinen GJ, Uijens MCM, et al. Molecular markers of prognosis in canine cortisol-secreting adrenocortical tumours. Vet Comp Oncol. 2019;17(4):545-552.
18. Hashimshony T, Senderovich N, Avital G, et al. CEL-Seq2: sensi- tive highly-multiplexed single-cell RNA-seq. Genome Biol. 2016; 17(1):77.
19. Chen J, Bardes EE, Aronow BJ, Jegga AG. ToppGene suite for gene list enrichment analysis and candidate gene prioritization. Nucleic Acids Res. 2009;37(SUPPL. 2):W305-W311.
20. Blighe K, Rana S, Lewis M. EnhancedVolcano: publication-ready vol- cano plots with enhanced colouring and labeling. 2021.
21. Marshall OJ. PerlPrimer: cross-platform, graphical primer design for standard, bisulphite and real-time PCR. Bioinformatics. 2004;20(15): 2471-2472.
22. Zuker M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 2003;31(13):3406-3415.
23. Schlotter YM, Veenhof EZ, Brinkhof B, et al. A GeNorm algorithm- based selection of reference genes for quantitative real-time PCR in skin biopsies of healthy dogs and dogs with atopic dermatitis. Vet Immunol Immunopathol. 2009;129(1-2):115-118.
24. Stassen QEM, Riemers FM, Reijmerink H, Leegwater PAJ, Penning LC. Reference genes for reverse transcription quantitative PCR in canine brain tissue. BMC Res Notes. 2015;8(1):761.
25. Vandesompele J, De Preter K, Pattyn I, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of mul- tiple internal control genes. Genome Biol. 2002;3(711):34-31.
26. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-4ACT method. Methods. 2001;25(4):402-408.
27. Zheng S, Cherniack AD, Dewal N, et al. Comprehensive pan-genomic characterization of adrenocortical carcinoma. Cancer Cell. 2016;29(5): 723-736.
28. . R2: Genomics Analysis and Visualization Platform. http:// r2platform.com
29. Brown GT, Cash BG, Blihoghe D, Johansson P, Alnabulsi A, Murray GI. The expression and prognostic significance of retinoic acid metabolising enzymes in colorectal cancer. PLoS One. 2014;9(3):90776.
30. Veit JGS, De Glas V, Balau B, et al. Characterization of CYP26B1-selective inhibitor, DX314, as a potential therapeutic for keratinization disorders. J Invest Dermatol. 2021;141(1):72- 83.e6.
31. Hammer GD, Basham KJ. Stem cell function and plasticity in the nor- mal physiology of the adrenal cortex. Mol Cell Endocrinol. 2021;519: 111043.
32. Sanders K, Mol JA, Kooistra HS, Slob A, Galac S. New insights in the functional zonation of the canine adrenal cortex. J Vet Intern Med. 2016;30(3):741-750.
33. Rainey WE, Nakamura Y. Regulation of the adrenal androgen biosyn- thesis. J Steroid Biochem Mol Biol. 2008;108(3-5):281-286.
34. Rege J, Nakamura Y, Wang T, Merchen TD, Sasano H, Rainey WE. Transcriptome profiling reveals differentially expressed transcripts between the human adrenal zona fasciculata and zona reticularis. J Clin Endocrinol Metab. 2014;99(3):E518-E527.
35. Osanai M, Petkovich M. Expression of the retinoic acid-metabolizing enzyme CYP26A1 limits programmed cell death. Mol Pharmacol. 2005;67(5):1808-1817.
36. Das BC, Thapa P, Karki R, et al. Retinoic acid signaling pathways in development and diseases. Bioorg Med Chem. 2014;22(2): 673-683.
37. Armstrong JL, Ruiz M, Boddy AV, CPF R, ADJ P, Veal GJ. Increasing the intracellular availability of all-trans retinoic acid in neuroblastoma cells. Br J Cancer. 2005;92(4):696-704.
38. Szabó PM, Tamási V, Molnár V, et al. Meta-analysis of adrenocortical tumour genomics data: novel pathogenic pathways revealed. Onco- gene. 2010;29(21):3163-3172.
39. Nagy Z, Baghy K, Hunyadi-Gulyás É, et al. Evaluation of 9-cis retinoic acid and mitotane as antitumoral agents in an adrenocortical xeno- graft model. Am J Cancer Res. 2015;5(12):3645-3658.
40. Njar V. Cytochrome P450 retinoic acid 4-hydroxylase inhibitors: potential agents for cancer therapy. Mini Rev Med Chem. 2005;2(3): 261-269.
41. Ocaya PA, Elmabsout AA, Olofsson PS, Törmä H, Gidlöf AC, Sirsjö A. CYP26B1 plays a major role in the regulation of all-trans-retinoic acid metabolism and signaling in human aortic smooth muscle cells. J Vasc Res. 2010;48(1):23-30.
42. Njar VCO, Gediya L, Purushottamachar P, et al. Retinoic acid metabo- lism blocking agents (RAMBAs) for treatment of cancer and dermato- logical diseases. Bioorg Med Chem. 2006;14(13):4323-4340.
SUPPORTING INFORMATION
Additional supporting information can be found online in the Support- ing Information section at the end of this article.
How to cite this article: Sanders K, Kooistra HS, van den Heuvel M, et al. Transcriptome sequencing reveals two subtypes of cortisol-secreting adrenocortical tumours in dogs and identifies CYP26B1 as a potential new therapeutic target. Vet Comp Oncol. 2023;21(1):100-110. doi:10.1111/vco.12871