Genetic characterization of Addison’s disease in Bearded Collies
BMC Genomics volume 21, Article number: 833 (2020)
Primary hypoadrenocorticism (or Addison’s disease, AD) is an autoimmune disease that results in destruction of the adrenal cortex and consequent adrenal insufficiency. The disease has been described in purebred and mixed breed dogs, although some breeds, including the Bearded Collie, are at increased risk for AD. Candidate gene approaches have yielded few associations that appear to be breed-specific. A single other genome-wide association study reported no significant regions of association for AD in Standard Poodles. The present study aimed to identify genomic regions of association for canine AD in Bearded Collies.
Our study consists of the first genome-wide association analysis to identify a genome-wide significant region of association with canine AD (CFA18). Peaks of suggestive association were also noted on chromosomes 11, 16 and 29. Logistic regression analysis supported an additive effect of risk genotypes at these smaller effect loci on the probability of disease associated with carrying a risk genotype on CFA18. Potential candidate genes involved in adrenal steroidogenesis, regulation of immune responses and/or inflammation were identified within the associated regions of chromosomes 11 and 16. The gene-poor regions of chromosomes 18 and 29 may, however, harbor regulatory sequences that can modulate gene expression and contribute to disease susceptibility.
Our findings support the polygenic and complex nature of canine AD and identified a strongly associated locus on CFA18 that, when combined with three other smaller effect loci, was predictive of disease. The results offer progress in the identification of susceptibility loci for canine AD in the Bearded Collie. Further studies are needed to confirm association with the suggested candidate genes and identify actual causative mutations involved with AD susceptibility in this breed.
Primary hypoadrenocorticism or Addison’s disease (AD; OMIA 000520–9615) results from immune-mediated destruction of the adrenal cortex, leading to a life-threatening clinical condition characterized by inadequate secretion of adrenocortical hormones. The disease has been reported in many purebred and mixed breed dogs [1, 2], with prevalence up to 0.5% in the overall canine population [3,4,5,6,7]. However, multiple breeds appear to be predisposed to AD [2, 8,9,10], including Bearded Collies, Portuguese Water Dogs (PWD), Standard Poodles, West Highland White Terriers, Leonbergers, Wheaten Terriers, and Nova Scotia Duck Tolling Retrievers (NSDTRs) [2, 5, 7, 11,12,13], and disease prevalence in these breeds can be as high as 9% [4, 14,15,16].
Similar to humans, studies on canine AD have suggested a female predisposition [7, 12], although males and females appear to be equally affected in Bearded Collies , PWD  and Standard Poodles . Average age of onset is 4 years for all dog breeds [2, 7, 11], and most dogs are diagnosed between 2 and 7 years of age . An exception is the NSDTR, for which average age of onset is 2 to 3 years [7, 13]. However, recent findings suggest that the NSDTRs are affected by two types of AD: a juvenile type often seen as part of a multisystemic autoimmune syndrome and diagnosed at a much earlier age, and an adult-onset AD that corresponds to the AD seen in other dog breeds .
Disease pathogenesis is similar in dogs and humans [2, 8,9,10, 19], with lymphocytic infiltration of adrenal tissue, autoantibody production against adrenocortical antigens and impaired steroidogenesis characterizing AD in both species [19,20,21,22]. Although autoantibodies against adrenal antigens are observed in > 80% of human AD patients , providing important diagnostic and predictive information, in dogs, autoantibodies have only been identified in 24% of AD patients , indicating limited usefulness as a diagnostic tool for canine AD.
Genetic predisposition for disease is evident in both dogs and humans [1, 8, 15, 16, 19, 24, 25], and genes implicated in increased susceptibility to human AD include the major histocompatibility complex (MHC) class II genes (haplotypes DR3-DQ2 and DR4-DQ8), cytotoxic T-lymphocyte-associated protein 4 (CTLA4), protein tyrosine-phosphatase non-receptor type 22 (PTPN22), MHC class II transactivator (CIITA) [19, 24, 26], and the autoimmune regulator (AIRE) gene, which causes an autoimmune polyglandular syndrome that includes AD . Although canine pedigree studies have indicated a recessive autosomal mode of inheritance for AD in the Standard Poodle and PWD [1, 15], a recent genome-wide association study (GWAS) on Standard Poodles failed to demonstrate significant associations with AD, suggesting a more complex mode of inheritance . Candidate gene approaches have implicated many of the same genes identified in human AD with increased susceptibility to canine AD: CTLA4 in the PWD, Cocker Spaniel and Springer Spaniel [10, 14, 28]; PTPN22 in the Cocker Spaniel ; and the canine MHC (or dog leukocyte antigen, DLA) class II genes in the NSDTR, Bearded Collie and Standard Poodle [29,30,31]. Particular single nucleotide polymorphisms (SNPs) in CTLA4 and collagen type IV alpha 4 chain (COL4A4) have been associated with AD in the Springer Spaniel, whereas SNPs in oxysterol binding protein like 9 (OSBPL9) and PTPN22 were implicated in the Cocker Spaniel, and SNPs in syntaxin-binding protein 5 (STXBP5) in the Labrador Retriever. Two different genes (the vitamin D binding protein group-specific component (GC) and the SP110 nuclear body protein (SP110)) were investigated for association with AD in Bearded Collies, though no statistical significance was seen after correcting for multiple testing .
Taken together, these studies indicate the existence of genetic commonalities and breed-specificity in canine AD, underscoring the need for additional research to elucidate the genetic basis of this disease in dogs. A 2019 health survey conducted by the Bearded Collie Foundation for Health showed a relatively high prevalence of autoimmune conditions in the breed (11.1% based on 3072 Bearded Collies included in the survey). Of these, 3.1% had been diagnosed with primary hypoadrenocorticism . Therefore, the present study aimed to identify genetic loci associated with AD in Bearded Collies by interrogating the genome. Given the likely polygenic nature of canine AD, we hypothesized that more than one genomic region would be associated with disease development and that those associated regions would harbor genes related to immune function and regulation.
Genome-wide association study
PLINK and GEMMA analyses on unrelated dogs identified a peak on canine chromosome (CFA) 18 (Fig. 1a and b), with a single SNP (BICF2P1230367 at 28,207,384 bp) reaching genome-wide significance.
To complement the findings in the GWA analyses, 37 additional minimally related dogs (3 male and 11 female AD cases, 11 male and 12 female controls) were genotyped for the top CFA18 SNP and the odds ratios (OR) calculated (Table 1; Additional Table 1); the total number of males and females analyzed across the cases and controls did not significantly differ (p > 0.05). Two genotypes were associated with increased risk for AD on CFA18; dogs homozygous for the alternate genotype (AA) or those with a heterozygous (AG) genotype were at higher risk for disease than those carrying a homozygous CanFam3.1 reference (GG) genotype (OR = 5.8, 95% CI = 2.62–12.92, p = 0.00001, and OR = 4.7, 95% CI = 1.53–14.70, p = 0.00742, respectively).
Recognizing the complex nature of AD, we anticipated smaller effect loci contributing to AD expression, though perhaps not reaching genome-wide significance. A suggestive threshold was therefore set at p ≤ 0.00005 [34, 35]. Peaks of suggestive association were noted on chromosomes 11, 16 and 29. The GEMMA analyses yielded similar observations after accounting for the genetic structure of the sample (Fig. 1c). Specifically, the same top SNPs were observed in both PLINK and GEMMA analyses, although the SNP on CFA11 did not quite reach the suggestive threshold in the latter. Odds ratios were also calculated for these suggestive SNPs in all 140 dogs, which included the additional 37 dogs mentioned above (Table 2).
A homozygous alternate genotype was associated with increased risk for AD at both CFA16 (i.e. GG) and CFA29 (i.e. CC) SNPs, both contributing similar risk for disease (OR = 4.0, 95% CI = 1.87–8.70, p = 0.00038; and OR = 4.2, 95% CI = 1.95–8.88, p = 0.00019, respectively). The highest OR obtained was associated with the homozygous alternate genotype on CFA11 (GG; OR = 21.0, 95% CI = 2.63–167.99, p = 0.00015). However, that risk genotype was observed in a single control dog and was only present in 20% of cases. Interestingly, all dogs carrying a risk genotype on CFA11 (GG) also carried the risk genotype on CFA16 (GG).
When the significant SNP on CFA18 and the suggestive SNPs on chromosomes 11, 16 and 29 were considered, as the number of risk genotypes increased, so did the risk for AD. That is, within these 140 dogs, dogs carrying risk genotypes at any one or more of these four evaluated SNPs were at greater risk for AD than dogs carrying only non-risk genotypes at all four SNPs (Fig. 2). Specifically, dogs carrying one or two risk genotypes were at greater risk for AD than dogs having none (OR = 17.0, 95% CI = 2.21–130.70, p = 0.00028), and dogs carrying risk genotypes at three or four of the SNPs were at greater risk for AD than those carrying one or two risk genotypes (OR = 32.9, 95% CI = 4.23–256.68, p = 9.03 × 10− 7).
To further investigate the observation that risk for AD increased with the number of risk genotypes a dog carried, logistic regression analysis was used to assess the genotypic contribution of each of the four SNPs of interest (the genome-wide significant CFA18 SNP and the suggestive SNPs on chromosomes 11, 16 and 29) on the probability of AD. Specifically, the analysis addressed the question of whether a multi-SNP model was better than the CFA18 single-SNP model in predicting the probability that a Bearded Collie would be affected by AD. Separate 2-SNP, 3-SNP and 4-SNP models were developed using the genotypes of the 55 Addisonian and 85 healthy control Bearded Collies. Lower leave-one-out information criterion (looic) values were noted for the multi-SNP models compared to the single-SNP ones. The lowest looic value was associated with the 4-SNP model, indicating that this model presented the best and most parsimonious fit for AD prediction in the study population (Additional Table 2). Based on the known genotypes and associated phenotypes, logistic regression calculated the probability of disease associated with every possible genotypic combination, even those that were not included in the data set (Additional Fig. 1), but only those associated with the highest risk genotypes are presented in the subsequent tables.
Taking the probability of disease associated with a CFA18 risk genotype as baseline, the addition of one other risk genotype at any of the three suggestive SNPs, in a 2-SNP interaction model, tended to increase the probability for AD, although a significant increase over the CFA18 SNP alone was only noted with the addition of a risk genotype on CFA11 (Table 3).
In contrast with the 2-SNP model, when considering a 3-SNP model, the addition of a risk genotype at any two suggestive SNPs also tended to increase the probability of AD in dogs carrying a risk genotype on CFA18, with significant increases noted for almost all possible 3-SNP risk genotype combinations. For the 4-SNP model, the addition of all four SNPs did not substantively increase the estimate of AD probability over the 3-SNP model likely reflecting that dogs that had a risk genotype on CFA11 also had the risk genotype on CFA16. Conversely, the probability of disease associated with carrying only non-risk genotypes at all four SNPs was extremely low for almost all possible 4-SNP combinations (Table 4).
To refine the regions of association for candidate gene exploration, haplotype analysis was done for the four chromosomes. Eighteen haplotype blocks were significantly associated with AD (seven on CFA11, six on CFA16 and five on CFA18; permuted p-value< 0.05; Additional Table 3), and seven of them overlapped protein-coding genes. The location and distribution of the seven haplotype blocks on CFA11 indicate the existence of two associated regions separated by less than 2 Mb on CFA11 (Fig. 3a): one at the 69 Mb range and the other at 71 Mb. Moreover, one of the haplotype blocks on CFA11 overlapped with a protein-coding gene of utmost importance to the immune system: the toll-like receptor 4 (TLR4) gene.
Similar to CFA11, five haplotype blocks on CFA16 clustered around the 27 Mb range (Fig. 3b), coinciding with the location of the top GWAS SNP (located within haplotype block 4), whereas a sixth haplotype block was located approximately 2 Mb upstream of that region. Moreover, four of the haplotype blocks overlapped protein-coding genes: zinc finger matrin-type 4 (ZMAT4), fibroblast growth factor receptor 1 (FGFR1), leucine zipper and EF-hand containing transmembrane protein 2 (LETM2), eukaryotic translation initiation factor 4E binding protein 1 (EIF4EBP1), adrenoreceptor beta 3 (ADRB3), glutamic-oxaloacetic transaminase 1 like 1 (GOT1L1), RAB11 family interacting protein 1 (RAB11FIP1) and BRF2 RNA polymerase III transcription initiation factor subunit (BRF2).
Unlike with chromosomes 11 and 16, none of the haplotype blocks on CFA18 contained the top GWAS SNP (Fig. 3c), which is located about half-way between haplotype blocks 4 and 5. Four of the haplotype blocks were clustered around the 27 Mb range, spanning a total of 1.2 Mb in a region containing various RNA genes and few protein-coding genes. The fifth block was located almost 11 Mb upstream of that region and overlapped with two protein-coding genes: reelin (RELN) and solute carrier family 26 member 5 (SLC26A5). No haplotype blocks were associated with AD on CFA29 (Fig. 3d).
The present study reports the first genome-wide association analysis to identify a statistically significant region of association with canine AD. Three suggestive smaller effect loci were also identified and shown to contribute to increased risk for AD in Bearded Collies. Previous studies focused on candidate gene approaches yielded few associations, most of which appear to be breed specific [2, 9, 10, 28, 32, 36]. Similar to human AD , MHC (or DLA) class II genes have shown the most consistent associations with canine AD to date, even though different DLA class II haplotypes have been associated with increased risk for the disease in different breeds [9, 29,30,31]. However, DLA class II haplotypes fail to completely explain AD development in both dogs and humans, and only 37.7% of Bearded Collies with AD possess the breed’s identified DLA class II risk haplotype for the disease . This may explain the lack of association with the DLA class II locus on CFA12 in the present study and underscores the importance of identifying additional loci contributing to AD risk.
Although sample size was a limitation of this study, a genome-wide significant region of association was identified on CFA18, revealing two risk genotypes for the top CFA18 SNP (AA and AG). Given the complex nature of canine AD, a single major locus associated with disease risk was not expected. Instead, complex traits are likely governed by multiple loci of smaller effects in addition to environmental factors, and SNPs representing these smaller effect genomic regions are not likely to reach genome-wide significance . Peaks reaching the suggestive threshold for disease association were identified on chromosomes 11, 16, and 29, and a homozygous alternate genotype at the top SNP of each chromosome was associated with increased risk for AD.
Grouping of the dogs based on the number of risk genotypes they carried at all four SNPs indicated that, as the number of risk genotypes increased, so did the risk for AD. Logistic regression analysis further supported a contribution of additional risk genotypes at these smaller effect loci (on chromosomes 11, 16 and 29) to the probability of disease associated with a risk genotype on CFA18, with a 4-SNP model presenting the best fit for predicting the probability of AD in the study population although a 3-SNP model showed similar performance and captured a greater number of AD cases. Notably, despite presenting the best fit, only ~ 10% of AD dogs carried risk genotypes at all 4 SNPs. While one might anticipate a greater number of AD cases to be captured by such a comprehensive model, especially when the study cohort consisted of a single dog breed expected to be genetically similar, complex inherited diseases such as AD may not require the presence of every single risk locus for disease development.
It should also be mentioned that the small sample size meant that many possible genotypic combinations were not represented in the dataset. Nevertheless, logistic regression analysis still allows for prediction of disease risk for all plausible genotype combinations, recognizing that the predictions in these “empty cells” is not very precise and the confidence intervals surrounding the probability estimates will be large. Based upon our assessment of genotype effects at these four chromosomes, there were genotype combinations that could be considered high risk for AD and others showing strong evidence of being associated with dogs not afflicted with AD, as shown in Additional Fig. 1. However, the amount of risk conferred by each SNP, or combinations of SNPs, remains unclear.
The highest OR was associated with a homozygous alternate genotype for the top CFA11 SNP, and all 2-SNP, 3-SNP and 4-SNP genotypic combinations showing the highest probability of disease carried the CFA11 risk genotype. However, this genotype was only seen in 20% (11/55) of cases, which might suggest heterogeneity of our case group despite our strict inclusion criteria. Given that all dogs carrying the risk genotype on CFA11 also carried a risk genotype on CFA16, it is more likely that the very high risk for AD seen in dogs carrying the CFA11 risk genotype is due to an additive effect of the two risk genotypes. This fact may also have confounded the 3-SNP and 4-SNP models, and may be why no significant differences were seen between the probability of disease in dogs carrying three or four risk genotypes.
The genome-wide significant SNP identified in our study lands in a gene-poor region of CFA18. The closest protein-coding gene is LRRC4C (leucine rich repeat containing 4C, also known as netrin-G1 ligand, NGL-1), located approximately 600 kb downstream of the SNP and involved in thalamocortical axon guidance during embryonic development . No other protein-coding genes were identified within 1 Mb upstream or downstream of the SNP. Based upon the literature, LRRC4C has no connection to immune function, inflammation, or adrenal development. However, long noncoding RNA (lncRNA) genes can be found approximately 500Kb upstream and downstream of the associated SNP. Long noncoding RNAs have limited or no protein-coding capacity, but can regulate various biological processes, including immune cell development and function. These molecules can also act as transcriptional regulators of neighboring genes, and recent studies have implicated lncRNAs in autoimmune diseases such as systemic lupus erythematosus (SLE) and rheumatoid arthritis (RA) . Gene-poor regions of the genome have also been shown to harbor regulatory sequences that can bind transcription factors and regulate gene expression. In fact, more than half of the GWAS variants associated with human autoimmune diseases fall in enhancer elements, attributing much of the genetic risk for autoimmune conditions to the modulation of gene expression . Whereas the mapping of regulatory elements to the canine genome is not yet available, research teams affiliated with The Dog Genome Annotation Project (DoGA) have been working to generate such functional annotation to facilitate the study of canine diseases.
Haplotype analysis on CFA18 indicated a second region of association approximately 12 Mb upstream of the genome-wide significant SNP. The haplotype block identified in this region overlaps two genes RELN and SLC26A5. However, based on the literature, neither appear to be potential candidate genes for canine AD. SLC26A5 encodes an anion transporter molecule responsible for electromotility in mammalian auditory sensory hair cells , and mutations in the human gene have been associated with hearing loss (autosomal recessive deafness; OMIM 604943). RELN plays a critical role in brain formation, and mice with Reln deficiency (i.e. reeler mice) present with malformations of the cerebral cortex . Though immune abnormalities have been reported in reeler mice, they were not compatible with autoimmune diseases and characterized by reduced cytokine production, T cell proliferation and B cell responses .
The top CFA11 SNP was located within a lncRNA gene adjacent to TLR4. As mentioned above, lncRNAs may be involved in autoimmunity through the regulation of neighboring genes . Moreover, five haplotype blocks were identified in that region, one of which overlapped with TLR4. TLR4 is a membrane-bound pattern recognition receptor expressed in macrophages, dendritic cells and T lymphocytes, and binding of this receptor results in the activation of signaling pathways that lead to the production of pro-inflammatory cytokines, chemokines and costimulatory molecules that activate the adaptive immune response . Dysregulation of TLR4 signaling may result in overt activation of proinflammatory signals, thus increasing the risk for autoreactive lymphocyte activation. Autoimmune diseases that have been associated with TLR4 include RA, SLE and type 1 diabetes .
Two haplotype blocks associated with AD were also identified approximately 1.5 Mb upstream of the top CFA11 GWAS SNP. Two immune genes can be found within 50Kb upstream of the haplotype blocks: TNF superfamily member 15 (TNFSF15) and TNF superfamily member 8 (TNFSF8). TNFSF15 is a member of the tumor necrosis factor family of ligands that can influence the proliferation, activation and differentiation of immune cells, including T cells. TNFSF15 expression can be induced in monocytes, macrophages and dendritic cells through TLR signaling, for instance, and in T lymphocytes in response to T-cell receptor (TCR) binding under inflammatory conditions . Polymorphisms in this gene have been associated with several autoimmune diseases, such as RA, inflammatory bowel disease, Crohn’s disease , SLE, and Grave’s disease – a polygenic organ-specific autoimmune thyroid condition . Similarly, TNFSF8 (also known as CD30 ligand) also plays a role in modulating T cell responses, where binding of these ligands on T lymphocytes results in augmented cell proliferation and cytokine production . Despite its recognized potential role in autoimmunity, to the best of our knowledge, no associations between polymorphisms in TNFSF8 have been reported for any specific autoimmune disorder.
The top SNP on CFA16 was intergenic and located in a gene-rich region of the chromosome. Similar to chromosomes 11 and 18, haplotype blocks associated with AD suggested two regions of association approximately 1 Mb apart on CFA16. Four of the six haplotype blocks overlapped with protein-coding genes: ZMAT4, FGFR1, LETM2, EIF4EBP1, GOT1L1, RAB11FIP1, BRF2 and beta-3 adrenergic receptor (ADRB3). However, literature review showed no evidence of their direct connection to immune function or adrenal development, making them unlikely candidates for susceptibility to canine AD. While EIF4EBP1 and FGFR1 are involved in the PI3K-Akt signaling pathway, which does play a role in regulating immune cell effector functions through various mechanisms , a more likely candidate gene for AD susceptibility would be indoleamine 2,3-dioxygenase 2 (IDO2), one of two neighboring genes involved in tryptophan catabolism and modulation of immune responses. The better characterized gene, IDO1, is known to modulate the immune system by altering regulatory T cell populations. Although IDO2 is less efficient at catabolizing tryptophan, it is expressed in antigen presenting cells and has been recently implicated in driving B-cell mediated and T cell-dependent autoimmune responses [50, 51].
It should also be noted that a gene involved in steroidogenesis can be found 200Kb downstream of the top CFA16 SNP: the steroidogenic acute regulatory protein (STAR) gene. This gene encodes a mitochondrial transporter protein that transfers cholesterol from the outer to the inner mitochondrial membrane where steroidogenesis is initiated. Various inactivating mutations in STAR result in lipoid congenital adrenal hyperplasia, a severe (and usually fatal) condition that manifests within the first weeks to months of life. Severe STAR deficiency also causes failed testosterone production in utero resulting in males being born with female genitalia . While inactivating mutations in STAR are extremely unlikely to be involved in canine AD, normal steroidogenesis in the adrenal glands is thought to suppress antigen presentation within the adrenal cortex . It is, therefore, possible that mutations in STAR causing even a mild reduction in steroidogenesis could promote a breakdown of tolerance to adrenocortical antigens, thus facilitating the onset of autoimmunity.
No haplotype blocks associated with AD were identified on CFA29. The two top SNPs on that chromosome marked a 1 Mb region of the chromosome that contained three protein-coding genes, two of which were uncharacterized. No potential candidate genes for canine AD were identified in the associated region of this chromosome.
Susceptibility to complex human autoimmune diseases, including AD, is commonly associated with an accumulation of small effect variants in multiple immune-related genes that are usually insufficient to cause disease on their own . Moreover, complex traits such as AD are influenced by environmental factors that may act as triggers for disease expression in some but not all carriers of risk genotypes. Though based upon a relatively small sample size, our findings are consistent with these observations, and suggest a combination of risk genotypes in different genomic loci contributing to AD development in Bearded Collies. While a larger sample size may help further elucidate the contributions of each associated region, it is unlikely that we would be able to expand the sample to any great degree given the relatively small actual and effective population size of Bearded Collies . Incomplete penetrance and potential genetic heterogeneity of canine AD may also limit our ability to determine the amount of risk each region confers. Low prevalence, complex nature, incomplete penetrance and potential heterogeneity makes it challenging to investigate the genetic basis of human AD , and similar challenges exist with canine AD, even within predisposed dog breeds. Nonetheless, our study offers progress in the identification of susceptibility loci for canine AD in Bearded Collies.
Despite the small sample size, the present study identified a genome-wide significant region of association with canine AD on CFA18, in addition to suggestive regions of association on chromosomes 11, 16 and 29. Logistic regression supported an additive contribution of risk genotypes at each of the associated loci to the probability of developing AD in Bearded Collies. Associated loci overlap with potential candidate genes for AD susceptibility that have roles regulating adrenal steroidogenesis, immune responses and/or inflammation. These findings offer progress in the identification of susceptibility loci for canine AD Bearded Collies, though further studies are needed to confirm association of the identified candidate genes and determine actual causative mutations involved with AD susceptibility in this breed.
Blood or buccal swab samples from privately owned Addisonian and healthy Bearded Collies were submitted to our laboratory by the dog’s owner or veterinarian. Addisonian dogs (AD cases) included in this study were selected based on a veterinary diagnosis of primary hypoadrenocorticism that included an ACTH stimulation test (serum cortisol levels < 55 nmol/L before and after ACTH administration) and evidence of electrolyte imbalance (sodium to potassium ratio < 27:1) and had no other concurrent illnesses (diagnosed or suspected). The AD cases had an average age of AD diagnosis of 50 months (range 5 to 120 months). Dogs included as healthy controls were at least 10 years old (average 12.5 years), with no history of AD or other autoimmune diseases (diagnosed or suspected). Samples were subjected to DNA extraction as previously described , and DNA aliquots quantified using a Nanodrop® spectrophotometer and stored at − 20 °C until further analysis.
Genome-wide association study
One hundred and forty-three DNA samples, consisting of 59 AD cases (22 males, 37 females) and 84 controls (34 males, 50 females), were genotyped at GeneSeek (Lincoln, NE) using the Illumina CanineHD BeadChip (San Diego, CA) containing 173,662 SNPs. Reference alleles and SNP locations were based on the CanFam3.1 reference genome; alleles differing from the reference allele at any given location were referred to as alternate alleles. Individuals and SNPs with less than 95% call rates as well as SNPs with minor allele frequency less than 1% were removed from the dataset using PLINK version 1.9 . Autosomal SNPs with a minor allele frequency greater than 1% were then used to generate a genetic relationship matrix for the remaining dogs using GCTA ; closely related individuals were excluded from GWAS analysis based on a relatedness cutoff of 0.3. Of the 143 dogs that underwent GWA analyses, 36 (17 AD cases, 19 controls) were removed from the analysis based on a relatedness cutoff of 0.3, and four dogs (one AD case and three controls) were excluded from analysis due to poor genotyping. Association analyses using the 98,540 SNPs remaining after quality control was done on the 103 unrelated dogs (41 AD cases and 62 controls) using PLINK version 1.9  following exclusion of SNPs that deviated from Hardy-Weinberg equilibrium among controls (p > 0.0001) and/or with a minor allele frequency less than 5%. The univariate linear mixed model approach implemented by GEMMA  was also used to account for any remaining population structure. As previously noted , most of the SNPs used in a canine GWAS are in linkage disequilibrium with multiple other SNPs, thus the number of SNPs used for analysis does not accurately represent the number of independent tests. To avoid overcorrecting for multiple testing, genome-wide significance was established using a Bonferroni-corrected alpha value of 0.05 divided by the effective number of independent SNPs as determined by GEC (Genetic type 1 error calculator) . Based on published studies addressing the genetics of complex traits in dogs, a threshold for suggestive association was set at p ≤ 0.00005 [34, 35].
To narrow down regions of interest for candidate gene exploration, haplotype analysis of the GWAS data was done using Haploview version 4.2  and 25,000 permutations on the chromosomes containing SNPs with suggestive and/or genome-wide significant association with AD. Genes overlapping with haplotype blocks were explored for potential involvement in AD development.
Genotyping of the four GWAS SNPs of interest on additional dogs
Polymerase chain reaction (PCR) was used to genotype significant and suggestive SNPs of interest (CFA11: BICF2G630307993, CFA16: BICF2G630109748, CFA18: BICF2P1230367 and CFA29: BICF2P977298) on additional dogs selected to be minimally related to those used for GWAS based on pedigree analysis, that is sharing no more than one grandparent. Primers flanking each of the SNPs were designed using Primer3  (CFA11-F: TGCCACATTCCCTTCCTTCT, CFA11-R: AGCATGACAAGACAGGACGA; CFA16-F: TTGTTTTGTTGTGTTTCCCTCTT, CFA16-R: TATCACAGTAAGGGGCCATAGG; CFA18-F: TTGGTTGTGCGTAGTCCTCT, CFA18-R: TCATCGCTCTAGTCACTGGG; CFA29-F: GCTAAGTACGCCTTGCAACC, CFA29-R: CCTGTGCTTGGATGTGATTG). A standard 31 cycle PCR protocol was used with an annealing temperature of 60 °C for the CFA11, CFA16 and CFA18 primer sets and 62 °C for the CFA29 primers. A 25 μL reaction and Promega GoTaq® Flexi DNA Polymerase (Promega, WI, USA) were used. Amplicon size was confirmed by running 6 μL of the PCR product on a 1% agarose gel. Exosap-IT™ express (Thermo Fisher Scientific, Waltham, MA, USA) was used to purify the PCR products according to the manufacturer’s recommendations, and these were then sequenced by capillary electrophoresis on an ABI 3730 DNA analyzer (Applied Biosystems, Foster City, CA, USA). Sequences were read using FinchTV v.1.4 (Geospiza, USA) and genotypes obtained by aligning the sequences against the CanFam3.1 reference sequence using CLC viewer v.7.8.1 (QIAGEN Bioinformatics, Denmark).
Odds ratio (OR) and two-tailed Fisher’s exact p-values were calculated using a 2 × 2 contingency table in VassarStats (http://vassarstats.net/odds2x2.html); OR calculations were based on the number of AD cases and healthy controls carrying a particular genotype compared to the number of cases and controls not carrying that genotype. Additionally, OR calculations included genotypic pairwise comparisons (i.e. homozygous alternate versus heterozygous; homozygous alternate versus homozygous reference; and heterozygous versus homozygous reference) to identify which genotypes conferred risk for AD. Statistical significance was set at p < 0.05.
Logistic regression was used to model the risk of disease as a function of the observed genotypes at each of the four SNPs of interest (CFA11: BICF2G630307993, CFA16: BICF2G630109748, CFA18: BICF2P1230367 and CFA29: BICF2P977298) across the unrelated dogs. Disease probability was defined by methods previously described . Briefly, the logit of disease probability was modeled as a function of the SNP genotypes representing 3 classes (e.g., AA, AG and GG with 2 degrees of freedom in the analysis). Genotypic effects and predictions of disease risk were determined using the Bayesian statistical package Stan  executed in R . A series of models were built using various combinations of the potentially important SNPs and the model with the best predictive accuracy for disease was identified using the leave-one-out cross-validation process , which calculates a statistic referred to as the leave-one-out information criterion (looic). Similar to that observed with the Akaike Information Criterion (AIC) for frequentist methods, the model with the smallest looic value is considered the best fit for predicting disease status. Looic calculation was performed for each of the models using the R package loo . The receiver operating characteristic curve, fitted with the R package pROC , was then used to assess the value of the models in disease prediction and define the optimal discrimination threshold for assigning a healthy or AD status to dogs based on each individual’s known SNP genotypes.
Availability of data and materials
All data generated or analyzed during this study are included in this published article and its supplementary information files.
Genome wide association study
Online Mendelian inheritance in animals
Portuguese Water Dog
Nova Scotia Duck Tolling Retriever
Major histocompatibility complex
Cytotoxic T-lymphocyte-associated protein 4
Protein tyrosine-phosphatase non-receptor type 22
MHC class II transactivator
Single nucleotide polymorphism
Oxysterol binding protein like 9 gene
Dog leukocyte antigen
Collagen type IV alpha 4 chain
Syntaxin-binding protein 5
Vitamin D binding protein group-specific component
SP110 nuclear body protein
The toll-like receptor 4
Zinc finger matrin-type 4
Fibroblast growth factor receptor 1
Leucine zipper and EF-hand containing transmembrane protein 2
Eukaryotic translation initiation factor 4E binding protein 1
Adrenoreceptor beta 3
Glutamic-oxaloacetic transaminase 1 like 1
RAB11 family interacting protein 1
BRF2 RNA polymerase III transcription initiation factor subunit
Solute carrier family 26 member 5
Leucine rich repeat containing 4C
Systemic lupus erythematosus
The Dog Genome Annotation Project
Online Mendelian inheritance in man
Long noncoding ribonucleic acid
TNF superfamily member 15
TNF superfamily member 8
Beta-3 adrenergic receptor
Indoleamine 2,3-dioxygenase 2
Indoleamine 2,3-dioxygenase 1
Steroidogenic acute regulatory protein
Polymerase chain reaction
Akaike information criterion
Oberbauer A, Bell J, Belanger J, Famula T. Genetic evaluation of Addison's disease in the Portuguese water dog. BMC Vet Res. 2006;2(1):15.
Boag AM, Catchpole B. A review of the genetics of Hypoadrenocorticism. Top Companion Anim Med. 2014;29(4):96–101.
Bellumori TP, Famula TR, Bannasch DL, Belanger JM, Oberbauer AM. Prevalence of inherited disorders among mixed-breed and purebred dogs: 27,254 cases (1995-2010). J Am Vet Med Assoc. 2013;242(11):1549–55.
Decôme M, Blais M-C. Prevalence and clinical features of hypoadrenocorticism in great Pyrenees dogs in a referred population: 11 cases. Can Vet J= La revue veterinaire canadienne. 2017;58(10):1093–9.
Hanson JM, Tengvall K, Bonnett BN, Hedhammar Å. Naturally occurring adrenocortical insufficiency – an epidemiological study based on a Swedish-insured dog population of 525,028 dogs. J Vet Intern Med. 2016;30(1):76–84.
Wiles BM, Llewellyn-Zaidi AM, Evans KM, O'Neill DG, Lewis TW. Large-scale survey to estimate the prevalence of disorders for 192 kennel Club registered breeds. Canine Genet Epidemiol. 2017;4:8.
Klein SC, Peterson ME. Canine hypoadrenocorticism: part I, vol. 51; 2010.
Charmandari E, Nicolaides NC, Chrousos GP. Adrenal insufficiency. Lancet. 2014;383(9935):2152–67.
Massey J, Boag A, Short AD, Scholey RA, Henthorn PS, Littman MP, et al. MHC class II association study in eight breeds of dog with hypoadrenocorticism. Immunogenetics. 2013;65(4):291–7.
Short AD, Boag A, Catchpole B, Kennedy LJ, Massey J, Rothwell S, et al. A candidate gene analysis of canine Hypoadrenocorticism in 3 dog breeds. J Hered. 2013;104(6):807–20.
Burkitt JM. Chapter 76 - Hypoadrenocorticism. In: Silverstein DC, Hopper K, editors. Small animal critical care medicine. Saint Louis: W.B. Saunders; 2009. p. 321–4.
Van Lanen K, Sande A. Canine Hypoadrenocorticism: pathogenesis, diagnosis, and treatment. Top Companion Anim Med. 2014;29(4):88–95.
Hughes AM, Nelson RW, Famula TR, Bannasch DL. Clinical features and heritability of hypoadrenocorticism in Nova Scotia duck tolling retrievers: 25 cases (1994-2006). J Am Vet Med Assoc. 2007;231(3):407–12.
Chase K, Sargan D, Miller K, Ostrander EA, Lark KG. Understanding the genetics of autoimmune disease: two loci that regulate late onset Addison's disease in Portuguese water dogs. Int J Immunogenet. 2006;33(3):179–84.
Famula TR, Belanger JM, Oberbauer AM. Heritability and complex segregation analysis of hypoadrenocorticism in the standard poodle. J Small Anim Pract. 2003;44(1):8–12.
Oberbauer AM, Benemann KS, Belanger JM, Wagner DR, Ward JH, Famula TR. Inheritance of hypoadrenocorticism in bearded collies. Am J Vet Res. 2002;63(5):643–7.
Peterson ME, Kintzer PP, Kass PH. Pretreatment clinical and laboratory findings in dogs with hypoadrenocorticism: 225 cases (1979-1993). J Am Vet Med Assoc. 1996;208(1):85–91.
Brown E, Young A, Wolf Z, Wade CM, Hughes A, Foreman O, et al. Juvenile Hypoadrenocorticism in the Nova Scotia duck tolling retriever: a recessive monogenic autoimmune disease. Denver: American College of Veterinary Internal Medicine Forum: 2016; 2016. p. 1407–519. Journal of Veterinary Internal Medicine.
Mitchell AL, Pearce SHS. Autoimmune Addison disease: pathophysiology and genetic complexity. Nat Rev Endocrinol. 2012;8(5):306–16.
Boag AM, Christie MR, McLaughlin KA, Syme HM, Graham P, Catchpole B. Autoantibodies against cytochrome P450 side-chain cleavage enzyme in dogs (Canis lupus familiaris) affected with Hypoadrenocorticism (Addison’s disease). PLoS One. 2015;10(11):e0143458.
Soderbergh A, Winqvist O, Norheim I, Rorsman F, Husebye ES, Dolva O, et al. Adrenal autoantibodies and organ-specific autoimmunity in patients with Addison's disease. Clin Endocrinol (Oxf). 1996;45(4):453–60.
Friedenberg SG, Brown DL, Meurs KM, Law JM. Lymphocyte subsets in the adrenal glands of dogs with primary Hypoadrenocorticism. Vet Pathol. 2018;55(1):177–81.
Boag AM, McLaughlin K, Christie M, Graham P, Syme H, Catchpole B. Analysis of P450 side-chain cleavage enzyme autoantibodies in dogs affected with hypoadrenocorticism. (2014 ACVIM forum research abstracts program). J Vet Intern Med. 2014;28(3):976–1134.
Betterle C, Morlin L. Autoimmune Addison's disease. Endocr Dev. 2011;20:161–72.
Løvås K, Husebye ES. Addison's disease. Lancet. 2005;365(9476):2058–61.
Skinningsrud B, Husebye ES, Pearce SH, McDonald DO, Brandal K, Wolff AB, et al. Polymorphisms in CLEC16A and CIITA at 16p13 are associated with primary adrenal insufficiency. J Clin Endocrinol Metab. 2008;93(9):3310–7.
Friedenberg SG, Lunn KF, Meurs KM. Evaluation of the genetic basis of primary hypoadrenocorticism in standard poodles using SNP array genotyping and whole-genome sequencing. Mamm Genome. 2017;28(1–2):56–65.
Boag AM, Short A, Kennedy LJ, Syme H, Graham PA, Catchpole B. Polymorphisms in the CTLA4 promoter sequence are associated with canine hypoadrenocorticism. Canine Med Genet. 2020;7(1):2.
Gershony LC, Belanger JM, Short AD, Le M, Hytönen MK, Lohi H, et al. DLA class II risk haplotypes for autoimmune diseases in the bearded collie offer insight to autoimmunity signatures across dog breeds. Canine Genet Epidemiol. 2019;6(1):2.
Hughes AM, Jokinen P, Bannasch DL, Lohi H, Oberbauer AM. Association of a dog leukocyte antigen class II haplotype with hypoadrenocorticism in Nova Scotia duck tolling retrievers. Tissue Antigens. 2010;75(6):684–90.
Treeful A, Rendahl A, Friedenberg S. Sex-specific DLA class II haplotypes are associated with primary hypoadrenocorticism in standard poodle dogs. J Immunol. 2019;202(1 Supplement):73.79.
Short AD, Catchpole B, Boag AM, Kennedy LJ, Massey J, Rothwell S, Henthorn PS, Littman MP, Husebye E, Ollier B. Putative candidate genes for canine hypoadrenocorticism (Addison's disease) in multiple dog breeds Veterinary Record. 2014;175:430.
Bearded Collie Foundation for Health. BeaCon open health registry report; 2019.
Karlsson EK, Sigurdsson S, Ivansson E, Thomas R, Elvers I, Wright J, et al. Genome-wide analyses implicate 33 loci in heritable dog osteosarcoma, including regulatory variants near CDKN2A/B. Genome Biol. 2013;14(12):R132.
Baker LA, Kirkpatrick B, Rosa GJM, Gianola D, Valente B, Sumner JP, et al. Genome-wide association analysis in dogs implicates 99 loci as risk variants for anterior cruciate ligament rupture. PLoS One. 2017;12(4):e0173810.
Hughes AM, Bannasch DL, Kellett K, Oberbauer AM. Examination of candidate genes for hypoadrenocorticism in Nova Scotia duck tolling retrievers. Vet J. 2011;187(2):212–6.
Lin JC, Ho W-H, Gurney A, Rosenthal A. The netrin-G1 ligand NGL-1 promotes the outgrowth of thalamocortical axons. Nat Neurosci. 2003;6(12):1270–6.
Gao Y, Li S, Zhang Z, Yu X, Zheng J. The role of long non-coding RNAs in the pathogenesis of RA, SLE, and SS. Front Med (Lausanne). 2018;5:193.
Fu Y, Tessneer KL, Li C, Gaffney PM. From association to mechanism in complex disease genetics: the role of the 3D genome. Arthritis Res Ther. 2018;20(1):216.
Schaechinger TJ, Oliver D. Nonmammalian orthologs of prestin (SLC26A5) are electrogenic divalent/chloride anion exchangers. Proc Natl Acad Sci. 2007;104(18):7693–8.
Tissir F, Goffinet AM. Reelin and brain development. Nat Rev Neurosci. 2003;4(6):496–505.
Greenjohnson JM, Zalcman S, Vriend CY, Nance DM, Greenberg AH. Suppressed T cell and macrophage function in the "Reeler" (rl/rl) mutant, a murine strain with elevated cerebellar norepinephrine concentration. Brain Behav Immun. 1995;9(1):47–60.
Mohammad Hosseini A, Majidi J, Baradaran B, Yousefi M. Toll-like receptors in the pathogenesis of autoimmune diseases. Adv Pharm Bull. 2015;5(Suppl 1):605–14.
Liu Y, Yin H, Zhao M, Lu Q. TLR2 and TLR4 in autoimmune diseases: a comprehensive review. Clin Rev Allergy Immunol. 2014;47(2):136–47.
Bayry J. TL1A in the inflammatory network in autoimmune diseases. Nat Rev Rheumatol. 2010;6(2):67–8.
Aiba, Yoshihiro, and Minoru Nakamura. The role of TL1A and DR3 in autoimmune and inflammatory diseases. Mediators Inflamm. 2013;2013(258164):9. https://doi.org/10.1155/2013/258164.
Zhang M, Liu S, Xu J, Lv S, Fan Y, Zhang Y, Zhang Y, Wu Y, Su Y, Yu H, Song S, He J, Li H. TNFSF15 polymorphisms are associated with graves' disease and graves' ophthalmopathy in a Han Chinese population. Curr Eye Res. 2020;45(7):888–95. https://doi.org/10.1080/02713683.2019.1705494. Epub 2019 Dec 27.
Mackay F, Kalled SL. TNF ligands and receptors in autoimmunity: an update. Curr Opin Immunol. 2002;14(6):783–90.
Patel RK, Mohan C. PI3K/AKT signaling and systemic autoimmunity. Immunol Res. 2005;31(1):47–55.
Merlo LMF, DuHadaway JB, Grabler S, Prendergast GC, Muller AJ, Mandik-Nayak L. IDO2 modulates T cell–dependent autoimmune responses through a B cell–intrinsic mechanism. J Immunol. 2016;196(11):4487.
Merlo LMF, Pigott E, DuHadaway JB, Grabler S, Metz R, Prendergast GC, et al. IDO2 is a critical mediator of autoantibody production and inflammatory pathogenesis in a mouse model of autoimmune arthritis. J Immunol. 2014;192(5):2082.
King SR, Bhangoo A, Stocco DM. Functional and physiological consequences of StAR deficiency: role in lipoid congenital adrenal hyperplasia. Pediatr Adrenal Dis. 2011;20Karger Publishers:47–53.
Theofilopoulos AN, Kono DH, Baccala R. The multiple pathways to autoimmunity. Nat Immunol. 2017;18:716.
Bell JS. AKC Bearded Collie Stud Book & Genetic Diversity Analysis. In., vol. 2017http://www.beaconforhealth.org/Final_Stud_Book_Analysis%20with%20edit%20restrictions.pdf; 2017.
Rincon G, Tengvall K, Belanger JM, Lagoutte L, Medrano JF, Andre C, et al. Comparison of buccal and blood-derived canine DNA, either native or whole genome amplified, for array-based genome-wide association studies. BMC Res Notes. 2011;4:226.
Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience. 2015;4(1).
Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88(1):76–82.
Zhou X, Stephens M. Genome-wide efficient mixed-model analysis for association studies. Nat Genet. 2012;44(7):821–4.
Li MX, Yeung JM, Cherny SS, Sham PC. Evaluating the effective numbers of independent tests and significant p-value thresholds in commercial genotyping arrays and public imputation reference datasets. Hum Genet. 2012;131(5):747–56.
Barrett JC, Fry B, Maller J, Daly MJ. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics. 2005;21(2):263–5.
Koressaar T, Remm M. Enhancements and modifications of primer design program Primer3. Bioinformatics. 2007;23(10):1289–91.
Belanger JM, Famula TR, Gershony LC, Palij MK, Oberbauer AM. Genome-wide association analysis of idiopathic epilepsy in the Belgian shepherd. Canine Med Genet. 2020;7(1):12.
Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, et al. Stan: a probabilistic programming language. J Stat Softw. 2017;76(1).
R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020. https://www.R-project.org/.
Vehtari A, Gelman A, Gabry J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput. 2017;27(5):1413–32.
Vehtari A, Gabry J, Magnusson M, Yao Y, Bürkner P, Paananen T, Gelman A. loo: Efficient leave-one-out cross-validation and WAIC for Bayesian models. R package version 2.3.1. 2020. https://mc-stan.org/loo.
Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez J-C, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011;12(1):77.
The authors would like to gratefully acknowledge the infrastructure support of the Department of Animal Science, College of Agricultural and Environmental Sciences, and the California Agricultural Experiment Station of the University of California, Davis. The authors would also like to thank Kaisu Hiltunen for technical help and Elsa Sell and Linda Aronson for their invaluable help in recruiting Bearded Collie samples, as well as all the owners who have actively participated in our research projects, and every dog that contributed to the study.
This study was funded by the American Kennel Club Canine Health Foundation (#02488, #02187-MOU, and #1236-A) and The Bearded Collie Foundation for Health (BeaCon). The authors would also like to thank the Academy of Finland, Helsinki Institute of Life Science and Jane and Aatos Erkko Foundation, the Brazilian National Council for Scientific and Technological Development (CNPq) and the Ann T. Bowling Fellowship for financial support. The Bearded Collie Foundation for Health (BeaCon) contributed with sample recruitment. Funding agencies were not involved in study design, sample collection, data analysis and interpretation or writing of the manuscript.
Ethics approval and consent to participate
The study followed all applicable guidelines for the care and use of animals, and all procedures were in accordance with the ethical standards of the University of California, Davis (IACUC #20402) and University of Helsinki, Finland (permit ESAVI/6054/04.10.03/ 2012). Written informed consent was obtained from all dog owners participating in the study.
Consent for publication
The authors declare that they have no conflicts of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The original online version of this article was revised: the caption for Figure 1 was erroneously put in the main text in the original article.
Figure 1. Logistic regression analysis for the four SNPs of interest. Probability of AD associated with the 81 plausible 4-SNP genotypic combinations based on known genotypes and phenotypes of 140 Bearded Collies (55 AD cases, 85 controls). Logistic regression calculates probabilities for genotypic combinations that are not included in the data set, so one can still predict the disease risk for all 81 plausible genotypes recognizing that the predictions in these “empty cells” is not very precise and the confidence intervals surrounding the probability estimates will be large. Based upon the observed data, probability estimates in excess of 0.42 were associated with AD risk whereas genetic combinations that had a probability below 0.42 were associated with reduced AD risk. Thus, in this figure, 0.42 was set as the baseline zero value. Probability estimates in excess of that figure are depicted above the line and those lower are depicted as below the line. For example, the highest probability calculated was 0.974, which would be represented as 0.55 above the threshold value of 0.42.
Table 1. Genotypes of the Bearded Collies analyzed at the four SNPs of interest. Genotypes of the 140 Bearded Collies (55 AD cases, 85 healthy controls) analyzed at each of the four SNPs of interest: BICF2G630307993 on canine chromosome (CFA) 11, BICF2G630109748 on CFA16, BICF2P1230367 on CFA18 and BICF2P977298 on CFA29.
Table 2. Model assessment. Values for the leave-one-out information criterion (looic) and standard error associated with each of the predictive SNP models tested. The model with the smallest looic value is considered the best, most parsimonious, fit.
Table 3. Haplotype blocks associated with Addison’s disease in Bearded Collies. Haplotype blocks associated with Addison’s disease on canine chromosomes 11, 16 and 18 in 103 unrelated Bearded Collies (41 cases, 62 controls). Locations represent the position of the first and last single nucleotide polymorphism (SNP) in each block and are based on the CanFam3.1 reference genome.
About this article
Cite this article
Gershony, L.C., Belanger, J.M., Hytönen, M.K. et al. Genetic characterization of Addison’s disease in Bearded Collies. BMC Genomics 21, 833 (2020). https://doi.org/10.1186/s12864-020-07243-0