Transcriptomic analysis by RNA sequencing characterises malignant progression of canine insulinoma from normal tissue to metastatic disease

Insulinomas (INS) are the most common human and canine functioning pancreatic neuroendocrine tumours. The long-term prognosis for malignant INS is poor, because micrometastases are frequently missed during surgery. As human and canine malignant INS share clinical and histopathological features, dogs have been proposed as models for INS research. Using RNA-sequencing, we conducted a pilot study to better understand the underlying molecular mechanisms of canine INS. Normal canine pancreas and lymph node control tissues were compared with primary INS and INS-metastatic lymph nodes, revealing more than 3,000 genes differentially expressed in normal pancreas compared to primary INS. Only 164 genes were differentially expressed between primary INS and INS-metastatic lymph nodes. Hierarchical clustering analysis demonstrated similar genetic profiles in normal pancreas and early clinical stage primary INS, whereas late clinical stage primary INS resembled the genetic profile of INS-metastatic lymph nodes. These findings suggest that markers of malignant behaviour could be identified at the primary site of the disease. Finally, using the REACTOME pathways database, we revealed that an active collagen metabolism, extracellular matrix remodelling, beta-cell differentiation and non-beta-cell trans-differentiation might cause disease progression and hyperinsulinism in INS, identifying major pathways worthy of future research in this currently poorly controlled disease.


Tissue homogenisation and RNA extraction
Pancreas tissue was perfused with RNAlater TM (Qiagen, UK) using a 5 mL syringe with a 26gauge needle in order to obtain a homogeneous distribution of the RNA later throughout the tissue, which was then cut in small pieces and snap-frozen with liquid nitrogen. Mesenteric lymph nodes were collected and immersed in RNA-later and subsequently sectioned and snapfrozen. Portions of INS tissue (5-15 mm diameter) were removed from the central portion of the identifiable mass, placed in cryovials and snap-frozen in liquid nitrogen. All tissues were stored at − 80 °C before RNA extraction. Tissue was homogenised using rotary homogeniser (#SHM2, Stuart, UK) with an additional first step of beta-mercaptoethanol addition to decrease RNases and optimise RNA quality. Total RNA was extracted using the RNeasy Mini Kit TM (Qiagen) To prevent contamination of the samples with genomic DNA, an on-column DNase treatment was performed (Qiagen). RNA concentrations were quantified by spectrophotometry (NanoDrop ND-1000 TM , Isogen Life Sciences, The Netherlands). RNA integrity was evaluated using an Agilent 2100 Bioanalyzer TM (Agilent Technologies, Santa Clara, CA) according to their RNA integrity number (RIN) values. RIN values above 6.5 were considered acceptable.

cDNA and primers preparation
RNA was transcribed into cDNA using the Omniscript reverse transcription kit (Qiagen) according to the manufacturer's instructions. Where available, previously published primers were checked using NCBI Primerblast (https://www.ncbi.nlm.nih.gov/tools/primer-blast). For primer design, the NCBI nucleotide database was used. For novel primer design, primer melting temperatures (57-63 degrees) were chosen to be similar (forward and reverse), while avoiding primer hairpins, self-dimers and heterodimerisation with the 3'-end of the sequence.
The primer sequences proposed were analysed using the IDT OligoAnalyzer online programme (https://www.idtdna.com/calc/analyzer). Primers were produced by Eurofins Genomics (UK). Primers are listed in Supplementary Table 3.

Statistical analysis
P-values were adjusted for multiple hypotheses corrections using the Benjamini-Hochberg approach, which controls the FDR. Probesets with FDR <0.05 and log2 fold change >2 were considered differentially expressed. All analyses were carried out using R programming language (www.R-project.org).

Supplementary data
Supplementary Figure 1 Differential gene expression (DEG) analysis between normal mediastinal lymphatic tissue and metastatic lymph nodes. Smear plots and vulcano plots of normal pancreatic tissues versus primary insulinomas (A,B) and primary insulinomas versus metastatic lymph nodes (C,D). Smear plots show gene differential expression comparing log ratio (log2FC) versus abundance (logCPM). Vulcano plots show the distribution of differentially expressed genes with log2FC >2 and <-2 based on their False discovery rate (FDR). FC= fold change; CPM= count per million

Supplementary Figure 3 Differential gene expression (DEG) analysis between normal mediastinal lymphatic tissue and metastatic lymph nodes.
Smear plots (A) using log ratio versus abundance and violin chart plot (B) based on log2FC display the DEG highlighting the set of altered genes. The heatmap analysis of the 6000DEG using heatmap2 (C) and unsupervised matrix clustering (D) shows that the samples separate in two superclusters: one made of normal lymph nodes (Ln) and the second cluster of metastatic lymph nodes (ML).

Supplementary Figure 5 Differentially expressed genes overlap between different pairwise comparisons.
Venn diagram showing number of overlapping differentially expressed genes between the three pairwise comparisons analysed in this study. NP= Normal pancreas; PI= Primary insulinoma; NL= Normal lymph nodes; ML=Metastatic lymph nodes.

Supplementary Table 1. Clinical features and RNA-quality of control tissues and insulinoma samples used in this study.
RIN, RNA integrity numbers; DOD, death of disease (0=alive; 1=dead, LTF=lost to follow-up).

Breed
Age ( Tumour limited to the pancreas and size<2cm T2 Tumour limited to the pancreas and size 2-4 cm N: regional lymph nodes NX Regional lymph nodes cannot be assessed N0 No regional lymph node metastasis N1 Regional