Genome wide exploration of the methylome in aggressive B-cell lymphoma in Golden Retrievers reveals a conserved hypermethylome

ABSTRACT

Few recurrent DNA mutations are seen in aggressive canine B cell lymphomas (cBCL), suggesting other frequent drivers. The methylated island recovery assay (MIRA-seq) or methylated CpG-binding domain sequencing (MBD-seq) was used to define the genome-wide methylation profiles in aggressive cBCL in Golden Retrievers to determine if cBCL can be better defined by epigenetic changes than by DNA mutations. DNA hypermethylation patterns were relatively homogenous within cBCL samples in Golden Retrievers, in different breeds and in geographical regions. Aberrant hypermethylation is thus suspected to be a central and early event in cBCL lymphomagenesis. Distinct subgroups within cBCL in Golden Retrievers were not identified with DNA methylation profiles. In comparison, the methylome profile of human DLBCL (hDLBCL) is relatively heterogeneous. Only moderate similarity between hDLBCL and cBCL was seen and cBCL likely cannot be accurately classified into the subtypes seen in hDLBCL. Genes with hypermethylated regions in the promoter-TSS-first exon of cBCL compared to normal B cells often also had additional hyper- and hypomethylated regions distributed throughout the gene suggesting non-randomized repeat targeting of key genes by epigenetic mechanisms. The prevalence of hypermethylation in transcription factor families in aggressive cBCL may represent a fundamental step in lymphomagenesis.

KEYWORDS:

• BCL
• DLBCL
• lymphoma
• methylome
• MIRA
• sequencing
• Golden Retriever
• canine
• cancer

Introduction

Diffuse Large B Cell Lymphoma (DLBCL) is the most common subtype of lymphoma in humans and in dogs, accounting for 21% and 52% of all lymphomas, respectively [1]. DLBCL is the most common subtype of aggressive B-cell lymphoma in dogs. Canine BCL (cBCL) shows remarkable similarities to human DLBCL (hDLBCL); both are treated with and respond to CHOP-like chemotherapy protocols, share a common cytologic and histologic appearance, and share common dysregulated pathways including NF-kB [2–5]. The importance of the canine model for hDLBCL has been proven with the completion of clinical trials in companion dogs that have provided the necessary data to bring drugs into human clinical trials, including a topoisomerase 1 inhibitor (LMP744) [6], an autologous vaccine with heat shock protein peptide complexes with hydroxylapatite [7], and a PI3Kδ inhibitor (RV1001) [8].

Genomic sequencing can define the underlying nucleic acid changes that lead to cellular transformation and the malignant biology of tumours. Understanding these changes can be used to create non-invasive diagnostics and assist in tumour classification for prognostication. Genetic studies have helped define the normal development of B cells and initiation and maintenance of cancer cells in hDLBCL [9]. More recent sequencing studies in hDLBCL and cDLBCL/cBCL have evaluated the methylome. The canine and human methylomes (autosomes) are composed of a similar number of CpGs, ~24 million and ~25 million CpGs [10], respectively. The canine methylome has a greater number of CGIs, ~45 thousand versus ~30 thousand CGIs (UCSC genome browser), but the CGIs are smaller on average than the human methylome (p < 2.2e-16, Wilcoxon rank sum). Epigenetic programming is precisely controlled in normal B-cell development, but aberrant DNA methylation can lead to cancers across different cell types and species [11]. Focal hypermethylation has been utilized for subtyping hDLBCL. Supervised hierarchical clustering of hDLBCL with the HELP (HpaII tiny fragment enrichment by ligation-mediated PCR) promoter microarray showed a 91% predictive accuracy in classification of patients into the activated B cell (ABC) and germinal centre B cell (GCB) subgroups [12]. A similar later study divided hDLBCL into six clusters based on degree and direction of methylation differences to normal germinal centre B-cells [13]. ABC and GCB subtypes were distributed throughout the clusters in this study but showed skewing in their distributions. Human DLBCL subgroups with different prognoses are more clearly separated on the transcriptome and genome but the methylome can also be used to highlight the differences between the subgroups.

Genome-wide DNA methylation in cBCL that have included cDLBCL have been evaluated in only four studies to date [14–17] and only one study, using microarrays, evaluated methylation from lymph node tissue in dogs with cBCL [14]; none of these previous studies used purified B cells as a control. Similarities between cBCL and hDLBCL were seen at the pathway level. Pathways that were commonly hypermethylated in cBCL and hDLBCL included the homeobox and developmental pathways [13,18]. In general, although similarities in dysregulated genes at the genetic and epigenetic level can be shared, the genetic landscape is marked with differences [2,3,19,20]. For example, gene expression profiling (GEP) can separate hDLBCL into three categories, ABC, GCB, and unclassified. This hDLBCL cell of origin (COO) classification scheme has not resulted in consistent separation in cBCL [19,21]. The CD10, BCL6, and MUM1/IRF4 IHC markers that are used in hDLBCL as surrogates for COO GEP also fail to separate cBCL into COO categories [19]. Unlike in hDLBCL, BCL6 does not appear to be commonly expressed [19], and BCL2 and MYC expressions appear to be uniformly high and fail to be prognostic in cBCL [22]. Recently, four prominent genetic subtypes, accounting for ~50% of hDLBCLs, termed MCD (based on the co-occurrence of MYD88L265P and CD79B mutations), BN2 (based on BCL6 fusions and NOTCH2 mutations), N1 (based on NOTCH1 mutations), and EZB (based on EZH2 mutations and BCL2 translocations) were described [23]. Somatic mutations in these genes were not reported in studies of WGS in 6 dogs with cBCL [24], RNA-seq in 62 dogs with cBCL, or WES in 64 Cocker Spaniels and Golden Retrievers with cBCL [2]. A CD79B SNV was found in a targeted sequencing and WES study in 63 dogs with cBCL. Despite the genetic differences, the clinical similarities, such as response to CHOP-based chemotherapy protocols, have shown the value of cBCL as a model for hDLBCL.

Methylated-CpG island recovery assay and next generation sequencing (MIRA-seq) or methylated CpG binding domain sequencing (MBD-seq), is an affinity enrichment-based method to identify highly methylated areas using a whole genome-based approach. MIRA-seq was used in the current study for the advantages that the number of evaluable loci is not limited by space restrictions on a chip and because the procedure exploits the stark differences in areas of high-CpG-density such as CGIs [11] in cancer and normal tissue methylomes and is thus well suited for identifying areas of interest that can be subsequently interrogated with methods that allow base pair resolution of methylation status [25]. MIRA-seq was selected over more commonly used methods in human oncology due to the lower cost and greater CpG coverage of MIRA-seq compared to other comparable methods such as reduced representation bisulphite sequencing (RRBS), familiarity with analysing MIRA-seq data at our institution [26,27], and concordant performance of MIRA-seq compared to other methods [28,29]. The objective of this study was to characterize the genome-wide methylation patterns in aggressive cBCL in Golden Retrievers. The American Golden Retriever breed was selected due to the relatively low genetic variation and the high risk of developing lymphoma, estimated to be an 11.5% lifetime risk [30] and an overall relative risk of 3.12 when compared to mixed breed dogs [31]. The secondary objectives of this study were to determine if the methylome in cBCL is consistent within Golden Retrievers and to identify key epigenetically dysregulated genes by analysis of genic and intergenic regions.

Results

Samples and quality control

Lymph node aspirates from 23 Golden Retrievers with BCLs were included in the study, Table S1. The control group consisted of immunophenotypically sorted B cells from the blood and lymph nodes from seven dogs with palpably normal lymph nodes and blood smears. Normal B cells were defined as CD21+ cells with low forward and side scatter. MIRA-seq quality control parameters of libraries made from these samples were met and are included in the Supplementary Information and in Table S2. Three hypermethylated differentially methylated regions (DMRs) in cBCL, HOXA13, PITX2, and WT1, were evaluated at a higher resolution, in previously untested cBCL and normal cohorts, with the EpiTYPER, Figures 1 and S1 and Table S3. Methylation fraction was higher at each CpG in the cBCL compared to the normal cohorts, thus validating MIRA-seq results. Methylation profiles between normal T and B cells were compared as well for further validation of the MIRA-seq methods used in this study. Normal T and B cells had distinct methylation signatures that centred around T-cell receptor (TCR) and B-cell receptor (BCR) signalling (see detailed discussion in the Supplementary Information, Tables S4-6 and Figure S2).

Figure 1. Validation of the hypermethylated regions identified with MIRA-seq in cBCL was done with the EpiTYPER. Methylation fraction as determined by the EpiTYPER in the evaluable CpGs in the differentially methylated regions are shown for HOXA13. The black dots and triangles represent EpiTYPER data for cBCL and normal lymph node aspirate samples respectively. The bar plots represent the interquartile range and the horizontal line in the bar plots represent the median. The normal population included 10 dogs that did not develop lymphoma within a 1-year period.

More hypermethylated than hypomethylated DMRs were identified by MIRA-seq in cBCL

Normal B cells from peripheral blood and lymph node (n = 7) were compared to cBCL (n = 23) for differential methylation, Table S1. 13,982 hypermethylated and 4467 hypomethylated DMRs were identified in cBCL compared to normal B cells, Table S7 and S8. Considering only those in genes with an Ensembl gene ID, this corresponds to 8468 hypermethylated and 2127 hypomethylated genic DMRs. These DMRs represent 4271 hypermethylated and 1218 hypomethylated unique protein-coding genes (Ensembl gene symbols), 601 and 230 lincRNAs, 22 and 1 miRNAs, 15 and 7 pseudogenes, 1 and 0 snRNA, 2 and 0 miscRNA, 14 cfa-miR, and 1 (miRNA), respectively. In a previous canine microarray study in cDLBCL, 1011 hypermethylated and 183 hypomethylated DMRs were found [14].

Differentially methylated regions (DMRs) in cBCL were preferentially located in CGIs and the first exon

The mean and median size of the DMRs were 424 bp and 298 bp (range 98–5398 bp) for those hypermethylated and 275 bp and 198 bp (range 98–2898 bp) for those hypomethylated, Figure S3. The hypermethylated and hypomethylated DMRs were enriched in the 5ʹUTR, promoter-TSS, exon, 3ʹUTR and transcription termination site (TTS) and were not enriched in the intron and intergenic space, Tables S9a and S9b. The hypermethylated cBCL DMRs were distributed most frequently around the TSS, Figure S4, more specifically 217 bp (median) downstream of the TSS. The median, mean and range for the length of the first exon (all canine genes from Ensembl) in Canis familiaris is 174 bp, 311 bp, 0–8007 bp. The hypermethylated DMRs were thus commonly located in the first exon and intron. This pattern follows the distribution of CGIs in the canine genome, Figure S4. Of the hypermethylated cBCL DMRs, 48.9% of the DMRs were in a CGI, 20.7% were in the CGI shores (within 2000 bp of the CGI) and 30.3% were in the CGI shelves and the open sea, Figure S5. In comparison, CGIs and shores represent only 1.69% and 8.27%, respectively, of the genomic space. DNA methylation in CGI shores can be strongly associated with gene expression, and can be highly conserved across species and tissue types, likely due to regulatory elements that are located in these regions [32]. Distribution of the hypomethylated DMRs is described in the Supplementary Information and Figures S4 and S5.

cBCL has a distinct methylation signature from normal B cells

All of the normal B-cell samples separated from the cBCL samples with hierarchical clustering of the hypermethylated and hypomethylated DMRs, Figure 2, and with a random forest model with all of the hypermethylated and hypomethylated DMRs, Figures S6-8. The cBCL samples separated into two branches in the hierarchical clustering of the hypermethylated and hypomethylated cBCL DMRs, Figure 2. The mean Euclidean distance and Spearman correlation between the samples are included in Table S10a and b. The frequency for each branch was 26% and 74%. The two branches are likely due to a batch effect, since all of the libraries in the first branch were produced and sequenced at a later date and with a different DNA polymerase (Group B) compared to most of those in the second branch. Samples MIRA_BCL9, MIRA_BCL16, MIRA_BCL18, and MIRA_BCL19 were sequenced with the samples from the first branch, Group B, but clustered in the second branch, with all of the libraries that were produced and sequenced at an earlier date (Group A); see Supplementary Information for a detailed discussion and Figure 2. Comparison of the MEDIPS QC metrics between Groups A and B, yielded significant differences in non-specific read pull-down p = 0.007 (% of reads that do not contain CG) and depth of CG coverage (% CG’s with >5x coverage) p = 0.018, Table S11. The libraries in Group B had more contaminating non-CG reads, a lower depth of CG coverage and the DMRs were less methylated than Group A. This could be due to technical variability in producing or sequencing the libraries or to real differences in the methylomes of the Group A and B samples. Batch correction of the cBCL data with Combat-seq was unable to completely correct the batch effect, Figures 2 and S9 [33]. This is likely because Combat-seq was optimized for RNA-seq data, which typically has higher magnitudes of variations in read counts than MBD-seq. The batch corrected data were interpreted with caution since batch correction has been shown to introduce error in methylation datasets (e.g., microarray), and Combat-seq has not previously been applied to MBD-seq data [34]. The difference in methylation between cBCL subgroups in Golden Retrievers is likely small and the batch effect seen in this study was larger than the difference in methylation between cBCL subgroups (if subgroups are present). All of the normal T and B cells were sequenced with the polymerase used in the group B samples and thus were not affected by the batch effect.

Figure 2. Batch corrected log2 transformed hypermethylated and hypomethylated rpkm data. The cBCL samples separated into 2 branches. These branches likely reflect a technical batch effect and not a biological difference. Samples MIRA_BCL9, MIRA_BCL16, MIRA_BCL18 and MIRA_BCL19 were sequenced with the first branch but clustered in the second branch and with all of the libraries that were produced and sequenced at an earlier date (Group A). DMRs were identified from the MIRA-seq of 23 cBCL and 7 normal B cell samples. The heatmap shows unsupervised hierarchical clustering of the 13,982 hypermethylated and 4467 hypomethylated DMRs in cBCL with a fold change value equal to or greater than 2. DMRs were identified with a q < 0.001 (Benjamin-Hochberg correction). Log2 of 1+ RPKM transformation of the values of each DMR was used for visualization purposes.

Pathways that were overrepresented in both the hypermethylome and hypomethylome included neurologic and PI3K-AKT pathways

KEGG pathways that were hypermethylated in cBCL compared to normal B-cells included neurologic and cardiac pathways and the PI3K-Akt and Wnt/Hippo signalling pathways, Table 12a. GO terms that were hypermethylated included neurologic, cardiac, and Wnt pathways, Table S12b. Hypermethylation occurred in genes with homeobox, tyrosine kinase (e.g., EGF), Wnt, G-protein coupled receptor, immunoglobulin, voltage-gated calcium channels, T-box transcription factor, and neurotransmitter-gated ion channel protein domains, Table 12c. The KEGG pathway that was also hypomethylated was the PI3K-Akt signalling pathway, Table S13a. GO terms that were also hypomethylated included neurologic pathways, more specifically pathways involved with microtubules (e.g., semaphorin, plexin, neuropilin, kinesin genes) and neurotransmitter activity, Table S13b. The specific neurologic GO terms that were enriched in the hypermethylated genes were different than for the hypomethylated genes. Hypomethylated protein domains in common with the hypermethylated gene list included tyrosine kinase (e.g., EGF), G-protein coupled receptor, immunoglobulin, plexin, semaphorin but also uniquely include tyrosine phosphatase, ephrin, FERM (4.1 protein, ezrin, radixin, and moesin), and cadherin domains, Table S13c. Many of the overrepresented pathways and protein domains were common to both hypermethylated and hypomethylated data sets and were also found when both hypermethylated and hypomethylated data sets were combined and analysed, Tables S12-14.

The HOX, TBX, NKX, FOX, and ONECUT transcription factor families and the Wnt pathway were overrepresented in the cBCL hypermethylome

The homeobox transcription factor families, HOX, NKX, and ONECUT, and TBX transcription factor families, and the homeobox domains were notably absent in the hypomethylated genes in cBCL compared to normal B-cells while 76.5% (26/34), 66.7% (10/15), 92.3% (12/13), and 100% (3/3) of the known Ensembl HOX, TBX, NKX, and ONECUT genes, respectively, in the dog were hypermethylated, Figures 3 and S10 and Table S15. Of the 32 canine FOX genes, 22 were hypermethylated and only one was hypomethylated. Seventy-eight per cent (15/19) of the known canine WNT genes had hypermethylated DMRs and only 15.8% (3/19) had a hypomethylated DMR. The negative regulators of the canonical Wnt pathway, the secreted frizzled-related proteins (SFRPs), dickkopf proteins (DKKs) and WIFI were also hypermethylated and none were hypomethylated. Likely multiple aspects of the Wnt pathway are dysregulated via differential methylation in cBCL. Promoter DNA hypermethylation in these transcription factor families, as also seen in human cancers, support a stem cell origin of cancer [35,36].

Figure 3. Focal areas of hypermethylation were distributed throughout the HOX clusters. The HOXA cluster on chromosome 14 is shown. The top track has the methylation profiles of 23 cBCL (red) and 7 normal B cell (blue) samples overlain. There were clusters of hypermethylated differentially methylated regions (DMRs) in the HOXA cluster. There were no hypomethylated DMRs in this region. A hypermethylated region in HOXA13 was selected for validation. CpG islands are shown in the bottom track.

A methylation signature of genes hypermethylated in cBCL is stable and conserved between dog breeds and geography

Upon intersection (with Bedtools) of the DMR locations in this study compared to the locations of the DMRs identified in the cDLBCL microarray study performed in an Italian population of dogs of various breeds diagnosed with cDLBCL, 94.3% (953/1011) of the hypermethylated, and 23.5% (43/183) of the hypomethylated DMRs were common [14]. The most recent study from Hsu et al. could not be compared since a complete list of hypermethylated genes was not provided and the data was not found in the provided GEO accession.

A 10 DMR signature was able to distinguish normal B cells from cBCL confirming the reliability of DNA methylation to differentiate normal B cells from cBCL

The top 10 DMRs (10 hypermethylated DMRs) identified by a random forest model were used for a MDS plot and AUC curve, Figures S6-7 and Table S16. These 10 DMRs were located in the intergenic (5), intronic (2), 3ʹUTR (1), TTS (1) and exonic (1) regions. The mean decrease in Gini coefficient and mean decrease accuracy values were similar, and the variable importance curve was steep, Figure S8 and Table S16, indicating similar performance of DMRs, likely in the top few hundred, in a classifier. This 10-DMR classifier can be improved upon and tested on a new validation cohort of Golden Retrievers with cBCL to best evaluate the model. Pathway analysis was done on the top 100 DMRs (41 intergenic DMRs were not used in the analysis) and showed enrichment in the glucose homoeostasis and RNA polymerase II transcription factor activity GO processes, PI3K-Akt signalling pathway, and HOX domains, Table S17. Genes in the glucose homoeostasis pathway included IGF1R, CRTC2, SIDT2, and TCF7L2. Genes in the RNA polymerase II transcription factor activity pathways included homeobox genes (NKX2-4, POU4F3, MSX1, and OTX1), ATOH1, MYF6, NEUROG3, TCF7L2, TCF7, and ZMIZ1.

The hypermethylome of cBCL shared characteristics with different types of aggressive human lymphomas

Comparing hypermethylated genes to those found in studies of hDLBCL, 37% were similar to those found in a hDLBCL (all subgroups) [13], and 60% were similar to those found an aggressive hBCL study [18], Figure S11, Table S18, and S19. In Chambwe’s microarray study of hDLBCL, the 140 hDLBCLs were divided into 6 subgroups based on increasing magnitude of differences in methylation compared to normal germinal centre B cells controls. In that study, a higher difference in methylation compared to normal was associated with decreased survival and the clusters with lower methylation differences had a higher frequency of GCB hDLBCLs. The genes with hypermethylated promoters in common with our study and in all hDLBCL’s in the Chambwe study included CH25H, GRIN3A, GSC, MSX1, NEFM, NID1, PDE3A, PTHLH, PUS3, RASGRF2, SIX6, SP9, ZIC5, and ZNF462 [13]. RASGRF2 has also been reported to be hypermethylated in other cBCL studies [17,37]. The protocadherin family contains 15 members in the dog and 6 and 3 of the PCDH genes were hypermethylated in cBCL compared to normal lymphocytes in the current study and in a previous cDLBCL microarray study, respectively, and none were hypomethylated [14]. PCDH1, PCDH7, PCDH10, PCDH12, PCDH17 and PCDH20 were hypermethylated in cBCL in our study and have been reported to be hypermethylated in specific clusters of hDLBCL. Hypermethylation of PCDH17 also has been identified in human urological cancers [38,39], oesophageal cancer [40], colorectal carcinoma [41], and acute lymphoblastic leukaemia [42] and has been associated with a poor prognosis. PCDH10 is a frequent and early methylation target in many cancers including hDLBCL [43,44]. Hypermethylation frequently leads to decreased gene expression of these PCDHs but the functional role of PCDHs in carcinogenesis is unknown. The methylation profile of cBCL in our study did not clearly align with a particular methylation cluster in hDLBCL identified in Chambwe’s study but based on the overlap of specific hypermethylated genes, was most similar to the hDLBCL clusters with higher levels of methylation and worse prognosis [13], Table S18. Comparing enriched pathways, the cBCLs in our study most clearly aligned with a mid-methylation cluster (cluster C), which was characterized by hypermethylation of many important developmental genes, particularly homeobox and forkhead box family genes. Comparing these two hBCL studies to each other, only moderate similarity exists between hypermethylation profiles, suggesting that hBCL is a more heterogeneous disease than cBCL. The decreased epigenetic heterogeneity seen in cBCL may be due to decreased heterogeneity in the disease of cBCL compared to hBCL in general. In a study comparing the incidence of different types of lymphomas in dogs and humans, of the non-Hodgkin B cell lymphomas (excluding cases classified as ‘other’), 73% and 37% of cases were classified as DLBCL, respectively [1]. Genes that were recurrently hypermethylated in multiple studies in cBCL and hDLBCL included FGF5 [16], PITX2, TJP1, WNT2, and NEFL [13,14], Table S18. All of these genes except TJP1 were also hypermethylated in the human pan-cancer study [45] and therefore are likely not specific to a cancer type or species. Previous GEP, exome, and methylation studies in cBCL have failed to stratify cBCL into similar subtypes as hDLBCL and thus our study lends support to the hypothesis that cBCL shares features that are characteristic of the different subtypes of hDLBCL but does not clearly align with a particular subtype.

There was little overlap of the hypomethylome of cBCL with aggressive human lymphomas

We compared hypomethylated genes in our study to hypomethylated genes in a hDLBCL study and an aggressive hBCL study and 14% [13] and 0% [18] of the genes in the human studies were found in our study, Figure S12, Table S19. Sixty-seven per cent of these hypomethylated genes that were in common with our study were also hypermethylated in our study. None of these DMRs were in the promoter region and were located in introns, intergenic, and 3’-UTR regions.

Discussion

Using a genome-wide methylation sequencing approach, the methylation profile of cBCL in Golden Retrievers showed a conserved hypermethylome. Further, there was almost complete overlap of hypermethylated genes with an Italian microarray study of different dog breeds, where only 9.3% of the dogs were Golden Retrievers; in the current North American study all of the dogs were Golden Retrievers [37]. The Italian study also only included dogs diagnosed with DLBCL with histopathology and this current study included Golden Retrievers with BCLs including DLBCL and possibly other clinically aggressive subtypes. We conclude in this study that the hypermethylome in clinically aggressive cBCL is conserved in Golden Retrievers in the United States, but also propose that conservation is seen across breeds, geographical region, and aggressive cBCL subtypes. These results suggest that, despite suspected genetic differences between Golden Retrievers in the American and European lines [46,47,48], aggressive cBCL is a relatively epigenetically homogenous disease. Confirmation of this finding would include future studies that include hierarchical clustering of cBCL from different breeds and countries with data that were generated and processed with similar methods. Expansion of this study to also include T-cell lymphomas (e.g., aggressive peripheral T cell lymphoma and T-zone lymphomas) would explore the hypothesis that the methylome of aggressive cBCL is distinct from other subtypes of canine lymphoma.

This study also describes the hypomethylome but our conclusions are mainly focused on the hypermethylome since the MIRA-seq method is more reliable for distinguishing hypermethylated than hypomethylated regions in cancer compared to normal. Hypomethylated DMRs were less conserved than hypermethylated DMRs across cBCL and hBCL studies and may have a less significant functional consequence in general. In cancers, hypomethylation tends to be global and can affect extensive blocks resulting in chromatin reorganization [49] and chromosomal stability [50] while hypermethylation tends to be more targeted, affecting specific genes. Hypomethylated DMRs have been shown to be more heterogeneous than hypermethylated DMRs in the same patient [51] and thus are likely to be more dynamic across time, breed, region, and lymphoma subtype. The mean fold change of the RPKM values of the hypermethylated and hypomethylated DMRs were 6.92 (range 2.00–324.11) and 2.52 (range 2.00–69.43), respectively (Wilcoxon rank sum test, p < 2.2 e^-16). The lower fold change of hypomethylated DMRs compared to hypermethylated DMRs could contribute to the dynamic nature of the hypomethylome or alternatively would have required more power to find a statistical difference and would have been more subject to variability (e.g., technical, patient, environment, etc.).

A strength of this study is that hypermethylated and hypomethylated DMRs were queried throughout the genome and not restricted to pre-determined loci. The MIRA-seq method identified more DMRs than previous cBCL studies and provided the most comprehensive methylome of cBCL in Golden Retrievers compared to previous studies in cBCL. MIRA-seq has been previously shown to be more sensitive than other methods such as MeDIP-seq and reduced representation bisulphite sequencing (RRBS) microarrays [52]. MIRA-seq is also more specific for the 5-methyl cytosine (5-mC) base modification compared to bisulphite-based techniques that cannot distinguish alternative DNA modifications such as 5-hydroxymethylation of cytosine [53]. DMRs outside of the promoter have also been associated with functional consequences; for example, internal promoter, splice site, and TTS methylation has been linked to isoform switching [54,55]. In a study using the TCGA dataset, it was estimated that isoform switching caused by differential DNA methylation affected 10–21% of genes [54]. Methylation around the TSS was correlated with decreased isoform usage while methylation in the 4^th or further downstream exon bodies increased isoform usage by promoting elongation. It is expected similarly in cBCL that there is a functional consequence to the highly recurrent DMRs. We can speculate that these methylation differences contribute to favouring of alternative transcripts.

The HOX, FOX, ONECUT, TBX, WNT, and NKX genes were recurrently hypermethylated in cBCL [14,15,16,17] and may be candidate biomarkers for early disease detection in cBCL. The HOX, FOX, ONECUT, TBX, and NKX transcription factors are required for differentiation and are reversibly repressed by PRC2, a histone methyltransferase that can cause gene silencing, in human embryonic stem cells [56]. In comparison, cancer cells are more likely to have promoter DNA hypermethylation in these transcription factor families than non-PRC2 targets, supporting a stem cell origin of cancer [57,58]. As a result, the HOX, FOX, ONECUT, TBX, WNT, and NKX genes are commonly hypermethylated across human solid cancers and are candidates for a pan-cancer methylation signature [35, 36, 45, 59,60,61,62,63,64]. As far as the author is aware, the HOX, FOX, TBX, WNT, and NKX genes have been found to be hypermethylated in only subsets of hDLBCLs. HOX and FOX hypermethylation was characteristic in a particular hDLBCL cluster and was seen in the clusters with the highest magnitude of methylation changes compared to normal germinal centre B cells [13]. TBX, NKX, and WNT hypermethylation was also seen in subsets of hDLBCLs with the highest magnitude of methylation changes compared to normal germinal centre B cells but were not characteristic for a particular cluster [13]. ONECUT hypermethylation was not found in any hDLBCL cluster [13]. HOX, FOX, ONECUT, TBX, WNT, and NKX genes were not differentially methylated when all of the hDLBCL samples were compared to normal in this hDLBCL study [13] and thus these genes have not previously been implicated as a key pathway in hDLBCL. With the exception of the HOX genes, these pathways have also not been previously implicated in cBCL. Consequential mutations in the ONECUT, NKX, TBX, WNT, and HOX genes were rare in a whole genome and RNA-seq study of hDLBCL, instead, epigenetic dysregulation may be a prevalent feature in lymphoma [65].

The HOX genes are transcription factors that are responsible for specifying body position in the developing embryo and are also critical for growth and differentiation. The HOX genes have been recurrently implicated in cancer, most commonly as an oncogene, but also as a TSG [66]. As a TSG, the HOX genes have been implicated in breast and lung carcinomas [67]. In acute leukaemias, the HOX genes have been commonly implicated as an oncogene [66] but also less frequently as a TSG [68]. Widespread hypermethylation in the HOX gene clusters are a common finding in human pancreatic and lung carcinomas [60,67] and has been seen in canine acute leukaemias (data not yet published), human acute lymphoblastic leukaemias (data not yet published), B cell lymphomas such as human follicular [69] and mantle cell lymphoma [70] and canine melanoma [71]. Potential mechanisms for the HOX gene hypermethylation seen in this study may be the overexpression of histone methyltransferases EZH2 (H3K27me3) and SETDB1 (H3K9me2/3) or DNMT3A mutation [72]. EZH2 and SETDB1 can interact with DNA methyltransferases to cause hypermethylation. This mechanism has been shown in hDLBCL [73] and hAML [74] and EZH2 overexpression has been shown in cBCL [75]. HOX gene hypermethylation in cBCL may just be an inconsequential product of aberrant hyperactivity of histone methyltranferases, association with other epigenetic features (e.g., CpG density, open chromatin, bivalent domains) [76] or may have biologically significant roles, such as regulation of downstream enhancers [60], activation of naïve stem cell resetting programmes or a more permanent arrest of differentiation programmes [77]. ChIP-seq, transcriptomic, and mutational profiles in cBCL are needed to understand this prevalent finding in aggressive cBCL.

Limitations

The methylation status of particular CpGs cannot be evaluated with MIRA-seq and thus methylation at the base pair level of resolution was not determined in this study. Further resolution can be evaluated with methods such as the EpiTYPER, which was done for selected DMRs in this study. Histology was only done on two patients and the subtype of cBCL could not be confirmed in all cases. The flow cytometry profile inclusion criteria used in this study has been consistent with histologically diagnosed cDLBCL in previous studies [22,78], but this phenotype has not been directly compared to different types of cBCL. This profile also most consistently correlates with histologically diagnosed cDLBCL at the Clinical Immunology Lab at Colorado State University [79]. Since cDLBCL is the most common WHO histologic diagnosis in dogs with high-grade B-cell lymphomas, 93–96%, the dogs included in this study were most likely to have cDLBCL [79,80,81,82]. It is possible that subtypes would have become more evident with increased sample numbers. Restricting the background genomic variability to Golden Retrievers living in the United States was used to overcome sample size limitations in the present study. The methylation profiles could not be correlated with clinical outcome since the patients were not treated with the same drug protocol. Approximately 16–32.7% of patients with aggressive lymphoma (mostly B cell) do not obtain a complete remission with a CHOP-based chemotherapy protocol [83,84,85,86,87]. There is currently no test that can reliably identify these patients. Clinical correlations can answer the question of whether the methylome can identify these patients and if so are there epigenetic mechanisms that can explain chemoresistance. RNA was not specifically collected from most of the patients, and thus gene expression cannot be directly correlated with methylation profiles. High methylation cannot be definitively attributed to hypermethylation since amplifications cause a similar methylation profile. Future experiments should also include whole genome, exome, and RNA sequencing for a complete illustration of the genomic landscape.

Conclusions

Promoter hypermethylation was a prevalent and conserved feature in cBCL in Golden Retrievers. The epigenetic signature of cBCL in Golden Retrievers appeared to be similar to that identified in cDLBCL in different breeds in different countries. Support for subgroups in cBCL was not found in the methylation data. The methylome of cBCL was more homogenous than hDLBCL. The hypermethylome of cBCL did not share many similarities with hDLBCL as a whole but shared modest similarity with methylation subgroups within hDLBCL. Genes with hypermethylated DMRs in the promoter-TSS-first exon often also had multiple hyper and hypomethylated DMRs distributed throughout the gene suggesting multiple layers of epimutations. The prevalence of hypermethylation in transcription factor families in cBCL may represent a fundamental step in oncogenesis. Many of the genes that were hypermethylated in cBCL were also hypermethylated in many other human cancers. A pan-cancer hypermethylation signature that is able to detect a myriad of canine cancers likely exists.

Materials and methods

Patients

Samples were submitted to the Clinical Immunology Laboratory at Colorado State University Veterinary Teaching Hospital (CO), University of Missouri Veterinary Health Center (MO), and Texas A&M (TX). All dogs with lymphoma that were included in the study were client-owned dogs described by their owners as Golden Retrievers. Collection of samples at the University of Missouri was approved by the ACUC under protocol #9726. Samples contributed by the Clinical Immunology Laboratory at Colorado State University were from samples that remained after flow cytometry analysis. All dogs with lymphoma had diffuse peripheral lymphadenopathy and a predicted aggressive clinical course. All dogs had a flow cytometry profile consistent with cBCL, including a CD21+ lymphocytosis of intermediate to large lymphocytes and variable MHCII expression. Full staging diagnostics with abdominal ultrasound, thoracic radiographs, and bone marrow aspirate were not performed on all dogs. Lymph node samples were either biopsy or fine needle aspirate samples that were collected as part of the patient’s diagnostic plan. Dogs were excluded if they were ever treated with chemotherapy for lymphoma prior to sample collection.

The MIRA-seq lymphoma affected group included lymph node samples from 23 Golden Retrievers with BCL (fine needle aspirates). The MIRA-seq normal control dogs included CD21+ cells from 7 dogs (5 from lymph nodes and 2 from blood), CD4 + T cells from 4 dogs (all from lymph nodes), CD8 + T cells from 3 dogs (all from lymph node) and CD5+ cells from 2 dogs (all from blood). Normal lymph nodes were collected from healthy dogs undergoing terminal procedures for continuing education courses or from euthanized client owned dogs. The lymph nodes were dissociated, and single-cell suspensions were incubated with antibodies specific for canine CD5 (clone YKIX322.3, APC conjugated), canine CD21 (clone CA2.1D6, PE conjugated), canine CD4 (clone YKIX302.9, Pacific blue conjugated), and canine CD8 (clone YCATE55.9, FITC conjugated). Genetic variability from contamination with other cell types (e.g., T cells, other white blood cells, and connective tissue cells) was minimized by purifying lymphocyte subsets from normal control dogs by cell sorting. A clinical summary of the dogs included in the study can be found in Table S1.

Flow cytometry and cell sorting

Immunophenotyping by flow cytometry was performed on all lymphoma samples as previously described [79], Table S20. A diagnosis of cBCL was established as previously described [79]; in brief dogs with a CD21 lymphocytosis and forward scatter of >469 were included. CD4+, CD8+, and CD21+ populations, from normal lymph node suspensions and PBMCs, were purified using a MoFlo cell sorter (Beckman Coulter) with four-way sort capability. Purity of all populations was assessed by re-analysis after sorting and was greater than 90%.

PBMC processing for normal dogs

All normal dogs had a CBC with blood smear evaluation by a board-certified clinical pathologist to confirm the absence of atypical lymphocytes or a lymphocytosis at the time of blood PBMC collection. Blood samples were collected in mononuclear cell preparation (CPT) sodium citrate tubes (#362761, BD Biosciences) and processed according to the manufacturer’s protocol to isolate PBMCs.

Magnetic column based manual cell separation for normal dogs

If flow sorting was unavailable, the Miltenyi MACS system was used for separation of B and T cells from canine PBMCs. The manufacturer’s protocol was followed. Specifically, cells were run through a pre-separation filter (#130-041-407, Miltenyi), washed in PBS with 1% BSA, incubated with 1/100 dilution of mouse anti-canine CD21 (#MCA1781 R, Bio-Rad), washed with PBS with 1% BSA, incubated with goat anti-mouse IgG microbeads (130–048-402, Miltenyi), washed with MACS buffer (1:20 dilution of 130–091-376 with #130-091-222, Miltenyi), eluted in MS column (#130-041-201, Miltenyi).

DNA library preparation and sequencing

Genomic DNA was extracted with the DNeasy blood and tissue kit (Qiagen Sciences, Germantown, MD, USA). Spectrophotometric and fluorometric DNA quantification was done with the NanoDrop and Qubit platforms, respectively (NanoDrop products, Wilmington, DE, USA, and Life Technologies, Grand Island, NY, USA). One ug of genomic DNA was sonicated with the Covaris M220 focused ultrasonicator to ~350 bp fragments (Covaris Inc, Woburn, MA, USA). Clean-up was accomplished with 1.6x AMPure XP beads (Beckman Coulter Inc, Brea, CA, USA). Library preparation was performed with the NEBNext Ultra DNA Library Prep Kit for Illumina (#E7370, New England BioLabs Inc, Ipswich, MA, USA). The MethylCollector Ultra kit was used for the MIRA (#55005, Active Motif, Carlsbad, CA, USA). PCR enrichment and indexing was done with the NEBNext Multiplex Oligos for Illumina kit (#E7335, New England BioLabs Inc, Ipswich, MA, USA). PCR was performed on the C1000 Thermal Cycler (Bio-Rad, Hercules, CA, USA) with the following conditions, 98°C for 30s, 11 cycles of 98°C for 10s, 65°C for 30s, 72°C for 30s and a final extension of 72°C for 5 minutes. The conditions were revised to 98°C for 30s, 4 cycles of 98°C for 10s, 65°C for 75s, and a final extension of 72°C for 5 minutes, when the PCR buffer was changed to the Q5 Hot Start HF PCR Master Mix. Libraries were cleaned twice with 1x Agencourt AMPure XP beads (A63881, Beckman Coulter, Indianapolis, IN, USA). Qubit quantification and fragment analysis was done on the resulting library (Advanced Analytical Technologies, Ames, IA, USA). DNA libraries were pooled to 3–4 samples per lane on the HiSeq 2500 sequencer and sequenced as 100 bp single-end reads (Illumina, San Diego, CA, USA).

Quality editing and control

Reads were trimmed 15 bp at the 5ʹ end and filtered for quality (>30% N calls or >10% of the sequence with Phred quality < 20 were omitted) using Cutadapt v1.6 [88]. Trimmomatic v0.32 was used to further select high quality reads [89]. Reads were dropped if the read length <20. Reads were trimmed if the average quality per base < 20 (4-base sliding window). The ends of the reads were trimmed if the quality of the base was <10. Base quality was measured with the simplified log-likelihood approach. FastQC v0.11.2 was used to obtain summary QC statistics for the remaining reads.

Alignment

Remaining reads were then aligned to the Canis familiaris reference CanFam3.1 with BWA-MEM v0.7.10 with default parameters [90].

MEDIPS Analysis of MIRA-seq data

Data was analysed with the Bioconductor MEDIPS package, v1.34.0 [91]. The following parameters were specified: unique reads, reads were extended 300 nucleotides, and each chromosome was divided into 100 nucleotide fragments for binning purposes. Only the autosomes were analysed. Differential methylation was evaluated with the following parameters, Benjamini-Hochberg correction method, q-value of 0.001, EdgeR, ≥10 reads per 100 bp window. DMRs were annotated to Ensembl canfam3.1 version 96 (last modified on 3/13/19, 2:04:00 PM) with HOMER v4.8. Mitochondrial, scaffold, and sex chromosomal DNA were excluded. The promoter was defined as 1000 bp upstream and 100 bp downstream of the TSS (default in HOMER). DMRs were also manually annotated as in a promoter by two additional methods. Firstly, if the start or end of the DMR was within 500 bp upstream or 200 bp downstream of the TSS, the DMR was labelled as belonging to the promoter-TSS region. Secondly, if the DMR was in the region outlined with previous CAGE studies (Cap Analysis of Gene Expression) [92,93], the DMR was labelled as belonging to the promoter-TSS region. The Eukaryotic Promoter Database (EPD) has 7545 experimentally validated promoters, Cf_EPDnew_001 [94]. This database represents the CAGE data that was further analysed (Ensembl version 92) and verified to only contain promoters that are associated with protein coding genes, on full chromosomes and annotated in UCSC with known transcripts. Differential methylation was evaluated in two separate analysis, 1) normal B cells were compared to normal T cells to identify DMRs and 2) normal B cells were compared to cBCL to identify DMRs. Batch correction was done with Combat-seq, sva version 3.38.0, with default parameters [33].

Pathway analysis

Pathway analysis was done with open access software, David v6.8 and ConsensusPathDB v34 with default settings. Overrepresentation analysis based on the hypergeometric test was done with pathway-based sets in ConsensusPathDB with the following parameters (minimum overlap with input list = 2, p < 0.01).

Classifier analysis

To assess the predictive utility of DMRs and identify key predictors, a random forest classifier was constructed using hypermethylated and hypomethylated DMR rpkm values as input and tumour status (normal vs cBCL) as the target with the R randomForest package v4.6.14 [95]. The model was constructed with 100,001 trees and default parameters otherwise. Variable importance was calculated as the Mean Decrease Gini value. Performance was estimated using the built in internal out-of-bag error method and reported as sensitivity, specificity, and overall accuracy from the confusion matrix. Data were visualized with ggplot, varImpPlot, and MDSplot.

DMR Data visualization

Heatmaps were done with gplots in R. Hierarchical clustering was done with the Euclidian average method. The log2 of the reads per kilobase per million mapped reads (RPKM) values were used.

Validation of DNA methylation with the EpiTYPER

Primers for regions in the HOXA13, PITX2, and WT1 genes, Table S21, were designed using the Agena Biosciences EpiDesigner online tool (www.epidesigner.com). DNA was bisulphite converted in duplicate with the EZ DNA Methylation kit (Zymo Research) and pooled prior to methylation analysis. The EZ DNA Methylation kit protocol was followed with the following modifications, an alternative cycling protocol of [1] 95°C for 30s [2], 50°C for 15 min [3], 20 cycles of [1] and [2]. PCR reactions were performed in triplicate with the pooled bisulphite converted DNA with the following mixture for each 5 ul reaction, 2ul of template, 1x buffer (Agena Biosciences), 10 ng dNTP, 0.4 µM primer mix, and 0.2 unit DNA polymerase (Agena Biosciences). The SAP and MassCLEAVE reactions were done according to the manufacturer’s protocol. Methylation percentages were calculated based on the ratio of the unmethylated versus methylated peaks. A low and high methylation control were used for each primer. The low methylation control was produced from a histologically confirmed normal lymph node from a dog without cancer with two cycles of whole-genome amplification with the REPLI-g kit (Qiagen). The high methylation control was produced from a histologically confirmed normal lymph node from a dog without cancer and was purchased from EpigenDx.

Abbreviations

ABC: activated B cells	=
cBCL: canine aggressive B-cell lymphoma	=
cDLBCL: canine diffuse large B-cell lymphoma	=
CGI: CpG island	=
CpG: cytosine-phospho-guanosine dinucleotide	=
DMR: differentially methylated region	=
GCB: germinal center B cells	=
GEP: gene expression profiling	=
GR: Golden Retriever	=
hDLBCL: human diffuse large B-cell lymphoma	=
hBCL: human aggressive B-cell lymphoma	=
NHL: non-Hodgkin lymphoma	=
PBMC: peripheral blood mononuclear cells	=
TSS: transcription stop site	=
TTS: transcription termination site	=
UTR: untranslated region	=

Data availability

Data was deposited into SRA under accession PRJNA808112.

Acknowledgments

The authors would like to acknowledge Hans Rindt, Senthil Kumar, and Maren Fleer at the University of Missouri for technical assistance. The authors would like to acknowledge Obi L Griffith for bioinformatics assistance.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2022.2105033

Additional information

Funding

This work was supported by the Canine Health Foundation under Grant 01918-G.