All the individuals (n = 1,411) studied in this research were from the Canadian Center for Fur Animal Research (CCFAR) at Dalhousie University, Faculty of Agriculture (Truro, Nova Scotia, Canada), from 2013 to 2021. In 2013, an outbreak of AD occurred at CCFAR, resulting in most of the mink being dead or culled in the barn. The exact source of the outbreak was not determined definitively, but it was suspected that AMDV-contaminated feed and contact with wild animals carrying AMDV were the most likely causes. Thus, within 3 years of the disease outbreak, approximately 150 mink (120 dam and 30 sires) from six AD-positive farms in Nova Scotia (Canada), which were believed to be resilient to AD and have been phenotypically selected for AD-resilient mink for many years, were introduced and used as breeders at CCFAR. AD was first identified in the province of Nova Scotia in Canada in 1941 (Agriculture_and_Marketing_of_Nova_Scotia_Government, 1998), and many mink farms in the province started selecting AD-resilient mink based on specific AD-resilient traits (Farid and Ferns, 2017). The studied mink included five different color types: black (n = 177), demi (n = 542), mahogany (n = 527), pastel (n = 152), and stardust (n = 13). The color type was identified by experienced technicians at CCFAR at weaning.

Counterimmunoelectrophoresis (CIEP) and antigen-based enzyme-linked immunosorbent assay tests (ELISA-G) were used to measure the immune response of the studied mink to AMDV exposure. The tests were conducted using established protocols described by Hu et al. (2021). In brief, blood samples of the studied mink were collected in mid-November before selecting breeders and in mid-February before mating from 2013 to 2021. The blood samples were sent to the Animal Health Laboratory at the University of Guelph (Ontario, Canada) and Middleton Veterinary Services (Nova Scotia, Canada) for CIEP and ELISA-G tests, respectively. The CIEP tests were used to detect the existence of anti-AMDV antibodies, and the results were recorded as 0 (negative: none or extremely low antibody level detected) or 1 (positive: detectable antibody level). The ELISA-G tests were applied to measure amounts of antibodies against AMDV using optical density, and the test results included eight categories with 1-point increments from 0 (none or an extremely low level of antibodies) to 7 (extremely high antibody level).

The Kleiber ratio (KR) and feed conversion ratio (FCR) were used to measure the growth and feed efficiency of the studied mink, respectively. The body weights of the studied individuals were collected using the established protocols described by Do et al. (2021). In brief, the body weight (BW) of the mink was measured at both birth and weaning (around 6 weeks after birth) and every 3 weeks from 13 to 28 weeks after birth. The average daily gain (ADG) and mid-test metabolic BW (BW^0.75) were calculated by the following equations, respectively: A D G = F i n a l   B W − I n i t i a l   B W Number   of   days   on   the   test , B W 0.75 = I n i t i a l   B W + F i n a l   B W 2 0.75 , where final BW was the BW on the last day of the feeding trial and the initial BW was the BW at 13 weeks of age. The Kleiber ratio (KR) was calculated using the following equation: K R = A D G B W 0.75 .

The feed intake data on the studied mink were collected using the established protocols described by Davoudi et al. (2022). In brief, mink were raised individually in separate cages, and feed was distributed daily to cages. The daily feed intake (DFI) of each mink was measured by calculating the difference between the amount of feed left over and the feed provided. The individual DFI records obtained during the experiment were averaged to obtain the individual average daily feed intake (ADFI). The FCR was calculated using the following equation: F C R = A D F I A D G   .

The pelt quality tests were conducted using the same method described by Hu et al. (2021). The live grading of overall pelt quality (QUA) was performed to measure the qualities of the mink pelt when they were alive. The grading was conducted based on the North American Fur Auctions live animal grading procedure by a skilled technician. The grading focused on checking the color consistency, fur roughness, and overall gloss. The QUA was scored into three categories from 1 (poor) to 3 (best).

The female reproductive performance was measured using the same approach described by Hu et al. (2021). Female reproductive performance was measured and recorded by the technicians in CCFAR during each annual reproduction cycle from 2006 to 2021. In this study, the number of newborn kits that survived 24 h after birth was used to quantify the reproduction performance of dams under AMDV exposure.

Studied individuals were grouped into pairwise subgroups based on their immune response, general resilience, and female reproductive performance. Studied individuals with both CIEP and ELISA-G test results were grouped into low- or high-immune response subgroups based on their CIEP and ELISA-G test results. Individuals with zero ELISA-G scores and negative CIEP results were grouped into low-immune response subgroups, and individuals with 5–7 ELISA-G scores and positive CIEP results were grouped into the high-immune response subgroup (Table 1). In this study, we not only grouped the entire populations of individuals into low- or high-immune response subgroups but also the individuals within the same color type, which included black, demi, mahogany, and pastel color types (Table 1). For resilience indicator traits, two methods were used to group CIEP-positive individuals into pairwise groups (Table 2). Studied individuals with positive CIEP results that had BW, feed intake, and pelt quality records were grouped into resilient or susceptible subgroups. The CIEP-positive individuals, which had bottom 20% FCR (14.38–22.49), top 20% KR (7.14–9.17), and score 3 (high pelt quality) for QUA were grouped into the resilient subgroup, and the CIEP-positive individuals, which had top 20% FCR (41.46–72.28), bottom 20% KR (2.19–4.41), and score 1 (low pelt quality) for QUA, were grouped into the susceptible subgroup (Table 2). Meanwhile, studied CIEP-positive dams, which had records for the number of newborn kits that survived 24 h after birth, were grouped into low- or high-female reproductive performance subgroups. The CIEP-positive dams with less than four newborn kits that survived 24 h after birth were grouped into the low-reproductive performance subgroup, and the CIEP-positive dams that had more than nine newborn kits that survived 24 h after birth were grouped into the high-reproductive performance subgroup (Table 2).

Tongue tissues from studied individuals were collected before pelting. The DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany) was used to extract the DNA from the tongue tissue based on the manufacturer’s instructions. The NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE) was applied to measure the quantity and quality of the extracted DNA samples. The 260/280-nm readings for all samples ranged from 1.7 to 2.0. All samples had a final concentration of 20 ng and were finally genotyped by an Axiom Affymetrix Mink 70K SNP panel (Neogen, Lincoln, Nebraska, USA) (Do et al., 2024).

PLINK (Purcell et al., 2007) was used to conduct animal and SNP data quality control. SNPs which had a minor allele frequency lower than 1%, call rate lower than 90%, excess of heterozygosity higher than 15%, Mendelian error frequency larger than 5%, and were out of Hardy–Weinberg equilibrium with a very low probability (1 × 10⁻⁵) were excluded from the analyses. Meanwhile, mink, which had a call rate lower than 90%, were also removed from the dataset. After quality control, 26,406 SNPs and 1,411 animals remained for further analyses.

Three methods, including the pairwise fixation index (Fst) (Weir and Cockerham, 1984), nucleotide diversity (θπ) (Nei and Li, 1979), and cross-population extended haplotype homozygosity (XP-EHH) (Sabeti et al., 2007), were performed to detect the selection signatures. The Fst and θπ methods directly utilize the SNP genotype, while the XP-EHH method uses phased data. The Fst analysis was conducted for each SNP based on the method proposed by Weir and Cockerham (1984) using VCFtools software (Danecek et al., 2011). The Z-transformation was performed using the scale function in the R program (Team, 2022) to normalize the Fst values. All negative Fst values were set to zero. The Fst values of all SNPs were ranked, and the SNPs with the top 5% Fst values were considered candidate selection signatures. The θπ analysis was conducted for each SNP based on the method proposed by Nei and Li (1979) using VCFtools software (Danecek et al., 2011). The θπ ratios were computed as θπ (subgroup1)/θπ (subgroup2) for all pairs of groups and were then log2-transformed (log2 (θπ ratios)). The SNPs with the top 5% θπ ratio values were considered candidate selection signatures. The XP-EHH approach was calculated for each SNP using selscan software (Torres et al., 2018). The missing genotypes were removed using VCFtools software (Danecek et al., 2011), and the genotypes were phased using Beagle software (Browning et al., 2018) because selscan software can only handle the phased genotypes without missing genotypes. The original obtained XP-EHH values were normalized using the norm function within the selscan software program. Then, the pnorm function in the R program was applied to calculate the p-values of the normalized XP-EHH values. The p-values of the normalized XP-EHH values were log-transformed, and the SNPs with -log (p-value) more than two were considered candidate selection signatures. Only the SNPs detected as candidate selection signatures by at least two methods were used for genome annotation, Gene Ontology, and functional analysis.

The potential selection regions were defined by extending 350 kb both upstream and downstream of the candidate selection signatures. The regions were defined based on the previous study that suggested that linkage disequilibrium (r ² < 0.2) in the current studied American mink population did not exceed 350 kb (Hu et al., 2023). Genome annotation was conducted using bedtools software (Quinlan, 2014) referring to the genome assembly of Neogale vison (Karimi et al., 2022). The Gene Ontology (GO) terms, including biological process (GO:BP), cellular component (GO:CC), and molecular function (GO:MF), were assigned to annotated genes using PANTHER 14.1 (Thomas et al., 2003). The overrepresentation tests of annotated genes were conducted using Fisher’s exact test and adjusted by the false discovery rate (FDR) correction, and the terms with the FDR-adjusted p-value (q-value) <0.05 were considered the overrepresented terms. Meanwhile, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted using the clusterProfiler package (Yu et al., 2012) in the R program with FDR control.

The genome-wide distribution of selection signatures associated with the immune response trait is presented in Figure 1. Additionally, Figure 2 illustrates the selection signatures that overlap across the three methods. Supplementary Material S1 presents the candidate selection signatures (chromosome number and location on the chromosome) detected by each method. When considering the entire population of individuals, a total of 619 SNPs were detected as candidate selection signatures by at least two methods. Notably, 33 SNPs were detected by all three methods and were considered strongly selected candidates (Figure 2). Furthermore, when analyzing individuals within specific color types, 444, 512, 385, and 335 were detected as candidate selection signatures for black, demi, mahogany, and pastel color-type mink, respectively. In addition, 57, 31, 32, and 45 SNPs were detected by all three methods for black, demi, mahogany, and pastel color-type mink, respectively, highlighting strong selection signature candidates specific to the immune response trait within each color type (Figure 2).

The candidate genes annotated from the selection signatures in the whole population and different color types are listed in Supplementary Material S1. Figure 3 shows the overlapped annotated candidate genes among the whole population and different color types. A total of 1,611 candidate genes were annotated from the selection signatures detected from the whole population (Figure 3 and Supplementary Material S2A). For black, demi, mahogany, and pastel color-type individuals, 1,355, 1,645, 1,436, and 1,042 candidate genes were annotated, respectively (Figure 3). Among the candidate genes annotated from the whole population and different color types, many genes were found to be associated with the AD-characterized phenotypes, including the immune system process, growth, pigmentation (except for the black color-type), reproduction, and response to stimulus (Supplementary Material S2, Figures 3, 4). Table 3 provides a list of genes for all color types. However, no common gene was detected among all color types for immune response traits (Figure 3).

Figure 4 presents the functional classifications of candidate genes related to the immune response trait. The genes annotated from the whole population, black, demi, mahogany, and pastel color-type individuals were classified into 17, 18, 17, 20, and 18 GO:BP categories, respectively, where the top four biological processes were the cellular process, metabolic process, biological regulation, and response to stimulus in all cases (Figure 4). The cellular anatomical entity and protein-containing complex were the two cellular components detected from the whole population and all color types. Regarding the GO:MF classifications, both analyses for whole population individuals and black color-type individuals detected 11 GO:MF, while analyses for demi, mahogany, and pastel color-type individuals detected 12 GO:MF. The top four GO:MF categories for the whole population and all color types were binding activity, catalytic activity, transcription regulator activity, and molecular transducer activity (Figure 4).

Table 4 (GO:BP) and Table 5 (GO:CC and GO:MF) present the overrepresentations of candidate genes. A total of 27, 18, 50, 26, and 18 significant (q-value<0.05) overrepresented GO enrichment terms were detected from the whole population, black, demi, mahogany, and pastel color-type individuals, respectively (Tables 4 and 5). Among all detected significant (q-value<0.05) overrepresented GO enrichment terms, some of them were commonly detected in the whole population and all color types (Tables 4 and 5), including two GO:BP (detection of chemical stimulus involved in sensory perception of smell (GO:0050911) and sensory perception of smell (GO:0007608)), three GO:CC (cytoplasm (GO:0005737), intracellular anatomical structure (GO:0005622), and organelle (GO:0043226)), and one molecular function (olfactory receptor activity (GO:0004984)). However, some were only detected in a certain color type of mink or the entire population. For example, the biological processes of adaptive immune response (GO:0002250) and system development (GO:0048731) were only detected for the whole population. Several metabolism-related GO:BP (e.g., heterocycle metabolic process (GO:0046483), macromolecule metabolic process (GO:0043170), and nitrogen compound metabolic process (GO:0006807)) and several GO:CC (e.g., envelope (GO:0031975), protein-containing complex (GO:0032991), and synapse (GO:0045202)) were only detected in demi color-type mink. Two unique GO:MF, namely, catalytic activity (GO:0003824) and hydrolase activity (GO:0016787), were only detected in mahogany color-type mink.

The genome-wide distribution of selection signatures for general resilience and female reproductive performance traits are presented in Figure 5, and the overlapped selection signatures are presented in Figure 6. Supplementary Material S3 presents the SNPs detected as candidate selection signatures by each method and the SNPs detected as candidate selection signatures by at least two methods. For general resilience traits, 569 SNPs were detected as candidate selection signatures by at least two methods, and 57 SNPs were detected as candidate selection signatures by all three methods (Figure 6). For female reproductive performance traits, 526 SNPs were detected as candidate selection signatures by at least two methods, and 16 SNPs were detected as candidate selection signatures by all three methods (Figure 6).

The candidate genes annotated from the candidate selection are listed in Supplementary Material S3, and the functional classifications of candidate genes are shown in Figure 7. A total of 1,933 genes were annotated from the selection signatures for the general resilience trait (Supplementary Material S2B). Several annotated genes were related to AD resilient traits, including growth, immune system process, pigmentation, and reproduction. The functional classifications of the annotated genes resulted in 18 GO:BP, 2 GO:CC, and 11 GO:MF. The annotation of selection signatures related to the female reproductive performance trait resulted in a total of 1,538 genes (Supplementary Material S2B). Except for 10 genes related to reproduction, several other genes were related to some AD-resilience indicator traits, including growth, immune system process, and pigmentation. The annotated genes were classified into 18 GO:BP, 2 GO:CC, and 11 GO:MF.

The overrepresentations of candidate genes are presented in Table 6. For general resilience traits, nine significant (q-value<0.05) overrepresented GO:BP were detected, and they were involved primarily in development, cellular process, and sensory perception of smell. Meanwhile, nine GO:CC (mostly related to organelle) and three GO:MF (related to olfactory receptor activity or binding) were significant (q-value<0.05) in the overrepresentation tests of candidate genes related to the general resilience. The overrepresentation tests of candidate genes from female reproductive performance traits resulted in 10 significant (q-value<0.05) overrepresented GO:BP (mostly related to development and detection of stimulus), 6 significant overrepresented GO:CC (mostly related organelle), and 4 significant overrepresented GO:MF (related to olfactory receptor activity or binding).

The overlapped genes among immune response, general resilience, and female reproductive performance traits are presented in Figure 8. In brief, 1,347, 1,680, and 1,277 unique genes were annotated from the selection signatures related to immune response, general resilience, and female reproductive performance traits, respectively. Sixteen genes, namely, ARHGAP19 (chr2: 209,800,177–209,812,754 bp), COL14A1 (chr4: 19,536,084–19,755,397 bp), DEPTOR (chr4: 19,803,040–19,933,557 bp), EXOSC1 (chr2: 209,917,603–209,928,479 bp), FAM135B (chr4: 4,922,384–4,961,329 bp), FRAT1 (chr2: 209,829,018–209,831,565 bp), FRAT2 (chr2: 209,839,939–209,842,216 bp), LOC122905718 (chr4: 5,215,903–5,216,010 bp), MMS19 (chr2: 209,938,058-209,971,536 bp), MRPL13 (chr4: 19,501,767–19,518,965 bp), PGAM1 (chr2: 209,910,223–209,917,518 bp), PTCHD4 (chr1: 94,047,590–94,230,349 bp), RRP12 (chr2: 209,866,295–209,895,940 bp), TBX18 (chr1: 44,912,005–44,922,875 bp), UBTD1 (chr2: 209,972,082–210,024,677 bp), and ZDHHC16 (chr2: 209,928,631–209,937,825 bp), were detected from all three studied traits.

The significant (q < 0.05) KEGG pathways of candidate genes from immune response, general resilience, and female reproductive performance traits are listed in Table 7. For immune response traits, only one significant (q < 0.05) pathway, the longevity regulating pathway, was detected. For female reproductive performance traits, one significant (q < 0.05) pathway, the mitogen-activated protein kinase (MAPK) signaling pathway, was detected. No significant (q < 0.05) KEGG pathway was detected for the general resilience trait.

A total of 1,611 annotated genes were related to the immune response trait by considering the whole population (Supplementary Material S2A). Some of these annotated genes were related to the immune system process (Figure 4; Tables 4 and 5). A total of 23 genes, namely, CCL26, CD28, CGAS, DEF6, EPOR, FAS, FYB1, GGT1, HIC2, IL16, JAK3, MEIS1, MFAP3, PATZ1, RUNX2, SHFL, SIGLEC15, THEMIS, TNFRSF21, TOX, TREM2, YES1, and ZNF572, were related to the immune system process, which might play important roles in immune-mediated responses to AMDV infection. Three genes, namely, TNFRSF21 (chr1: 93,535,621-93,549,929 bp), CCL26 (chr4: 10,918,021-10,922,028 bp), and TREM2 (chr1: 87,946,316-87,950,960 bp), are related to inflammatory processes (Stubbs et al., 2010; Santer et al., 2012; Liu et al., 2020). This may be due to several inflammations, which include interstitial nephritis, myocarditis, hepatitis, splenitis, meningoencephalitis, pneumonia, glomerulonephritis, and arteritis, caused by AD infection (Jepsen et al., 2009). Four genes, SIGLEC15 (chr3:151,209,875–151,221,497 bp), JAK3 (chr6: 213,007,629–213,025,411 bp), DEF6 (ch1: 117,504,375-117,525,807 bp), and FAS (chr2: 164,464,382-164,489,036 bp), were related to autoimmune disorders in humans (Hsu et al., 2012; García-Bermúdez et al., 2015; Serwas et al., 2019; Läubli and Varki, 2020), and AD is defined as an immune complex-mediated disorder disease in mink (Bloom et al., 1988). Three genes, IL16 (chr13: 150,353,853-150,448,017 bp), THEMIS (chr1: 73,061,569-73,141,431 bp), and CD28 (chr3: 15,437,373-15,465,940 bp), were related to T-cell proliferation (June et al., 1987; Wilson et al., 2004; Fu et al., 2009), and the TOX (chr4: 75,470,800-75,728,160 bp) gene was detected to be a crucial transcription factor involved in the exhaustion of CD8⁺ T cells (Seo et al., 2019). The detection of those genes might be related to the proliferation of CD8⁺ T cells after AD infection as CD8⁺ T cells were found to double in numbers during the development of AD (Aasted, 1989). The EPOR (chr6: 216,150,503–216,156,409 bp) gene was discovered to be associated with the production of red blood cells, and severe anemia was observed in AD-infected mink few months after infection (McGuire et al., 1979). The FYB1 (chr1: 286,138,312–286,169,159 bp) gene was found to be related to thrombocytopenia (Levin et al., 2015), which is one of the typical symptoms of AD infection (Gordon et al., 1967). The CGAS (chr1: 114,736,330–114,760,223 bp) gene was related to the production of the type I interferons and activation of inflammasomes (Wang et al., 2017; Decout et al., 2021), and the increase in the number of type I interferons was observed in the host during AD infection (Jensen et al., 2003). Meanwhile, the overrepresentation tests on the annotated genes detected one significant (q < 0.05) GO:BP, adaptive immune response (GO:0002250), related to immune response, where eight genes (IL12B, TNFRSF21, TAP1, JAK3, TAP2, C7, THEMIS, and C6) were involved.

The immune-response-related genes detected in this study were different from the genes detected by a previous study (Karimi et al., 2021). Seven genes, namely, TRAF3IP2, WDR7, SWAP70, TNFRSF11A, CBFB, IGF2R, and GPR65, were detected and related to the immune system process, as explained by Karimi et al. (2021), and none of these genes were detected in the current study. Several potential reasons could lead to these discrepancies: 1) the use of different types of genomic data (GBS in their study vs genotypes in this study); 2) the use of different grouping methods, where kidney lesion levels and virus loads were also considered in grouping animals in their study in addition to antibody titer; 3) the different ways the animal contracted AMDV (intranasal inoculation in their study vs natural exposure in this study); and 4) the color types of studied mink (only black in their study vs multiple colors in this study).

A total of 20, 26, 19, and 9 genes were related to the immune response in black, demi, mahogany, and pastel color types, respectively (Table 3). Most of the genes detected from a single-color type were unique from the rest of the color types. For black, demi, and mahogany color types, there were few genes in common between the two color types, but no common gene was detected among all of them. For pastel, eight of nine detected genes (only TNFRSF1B gene was common with the black color type) were unique from the rest of the three darker color types, which might indicate that pastel color-type mink has different immune responses to AD infection compared with the other three darker color types of mink. This could be a potential reason to explain the previous finding by Ellis (1996), where the mink with lighter color types were observed to be more susceptible to the AMDV than darker mink.

The KEGG pathway analysis of annotated genes from the whole population or different color types only detected one significant KEGG pathway, the longevity regulation pathway, in mahogany color-type mink. The relationship between longevity and immune response is complex. A strong and well-functioning immune system is crucial for protecting an organism from infections and other threats, and therefore, may contribute to increased longevity (Xia et al., 2019). Furthermore, chronic inflammation and overactivation of the immune system have been linked to aging and age-related diseases, which can shorten lifespan (Rea et al., 2018). Aleutian disease is defined as an immune complex disease and can cause persistent and chronic infection in mink (Porter and Cho, 1980; Stolze and Kaaden, 1987). Thus, the detection of the longevity regulation pathway may be related to chronic infection and autoimmune disorders caused by AD.

Since the general response trait used in this study was a combination of three AD resilient traits, which include growth, feeding efficiency, and pelt quality, we focused on genomic regions containing genes related to these traits. A total of 1,933 genes were related to the general resilience trait (Supplementary Material S3B). From these, several annotated genes were related to body growth. For example, PRKAG3 (chr3: 29,115,382–29,120,999 bp), a regulatory subunit of the AMP-activated protein kinase, was detected in this study and found to be related to body growth in several livestock species including swine (Ryan et al., 2012), sheep (Ibrahim, 2015), and beef cattle (Li et al., 2012). The PLAG1 (chr1: 58,588,578–58,696,784 bp) gene was also detected in this study. This gene is a positive regulator of insulin-like growth factor 2 (Voz et al., 2000; Van Dyck et al., 2007) that is known to affect body weight in both livestock (e.g., swine (Van Laere et al., 2003) and beef cattle (Huang et al., 2013)) and humans (Sandhu et al., 2003). The TMEM18 (chr1: 8,940,691–9,011,524 bp) gene detected in this study has been reported to be associated with growth traits and obesity in rats (Rask-Andersen et al., 2012), cattle (Ma et al., 2012), and humans (Almén et al., 2010; Haupt et al., 2010). Meanwhile, three genes, TPRA1 (chr6: 165,580,535-165,594,267 bp), MCM2 (chr6: 165,600,097-165,621,492 bp), and Tbx18 (chr1: 44,894,084-44,922,875 bp), which were all found to be related to embryo development in mice (Aki et al., 2008; Wehn and Chapman, 2010; Xu et al., 2022), were also detected in this study indicating that AD may influence the early stages of mink development, and therefore, influence growth. Meanwhile, several genes related to feed efficiency were also detected, for example, MRAP2 (chr1: 45,464,087-45,513,205 bp) and GLP1R (chr1: 85,852,947-85,884,385 bp). MRAP2 (Berruien and Smith, 2020) and GLP1R (Dailey and Moran, 2013) were found to play important roles in regulating appetite, and AD has been reported to cause adverse influences on the appetite of infected mink (Jensen et al., 2016). The annotated gene, HCRTR2 (chr1: 101,039,715-101,147,858 bp), is an orexin receptor and plays an important role in feeding behavior and balance of energy metabolism (Spinazzi et al., 2006; Belkina and Denis, 2012). In addition, several annotated genes, including ESRRG (chr10: 14,353,012-14,953,383 bp), LZTFL1 (chr6: 208,512,050-208,525,767 bp), and ELOVL4 (chr1: 49,176,236-49,208,238 bp), were reported to play key roles in regulating metabolism processes (Zhang et al., 2001; Alaynick et al., 2007; Wei et al., 2018). In addition to the genes related to growth and feed efficiency, the DCT gene (chr5: 151,495,518-151,529,776 bp), related to pigmentation (Guyonneau et al., 2004), was also detected in this study, and this might be related to the hair depigmentation, which causes single white hairs in the fur (sprinklers) impacting the pelt quality of AD infected mink (Farid and Ferns, 2011).

A total of 1,538 genes were detected from female reproductive performance traits (Figure 7; Table 6). From these, several genes, including SLX4 (chr14: 18,318,274-18,338,014 bp), TDRD6 (chr1: 92,984,166-92,997,290 bp), TACR3 (chr11: 104,943,216-105,013,153 bp), SHOC1 (chr9: 21,365,597-21,443,498 bp), FBXW11 (chr1: 255,864,225-255,951,567 bp), EPC2 (chr3: 82,648,215- 82,734,211 bp), GSC (chr13: 10,120,630-10,122,701 bp), and DICER1 (chr3: 9,803,946-9,876,750 bp), were related to reproduction. SLX4 (Hamer and de Rooij, 2018), TDRD6 (Vasileva et al., 2009), SHOC1 (Zhang Q. et al., 2019), and FBXW11 (O’Doherty et al., 2018) play important roles in the development of germ cells. The gene TACR3 plays a key role in reproductive functions, and loss-of-function mutations in this gene can lead to hypogonadotropic hypogonadism and infertility in humans (Guran et al., 2009; Topaloglu et al., 2009; Young et al., 2010). EPC2, GSC, and DICER1 genes are all important for the development of early embryos in animals and have been related to the reproduction traits (e.g., litter size) in swine and cattle (Kaczmarek et al., 2020; Chen J. et al., 2022; Chen Z. et al., 2022; Wang et al., 2022). The reproduction-related genes detected in this study were different from the genes detected in the signature selection study for response to Aleutian disease by Karimi et al. (2021). In that study, the genes FBXO5, CATSPER4, GOT2, and CatSperβ were annotated from the candidate selection regions and related to reproductive performance, which were not detected in our study. The different genomic data, grouping methods, and population structures could be the potential reasons that lead to the differences between these studies.

The KEGG pathway analyses of annotated genes detected only one significant (q < 0.05) pathway, the MAPK signaling pathway. The MAPK signaling pathway is involved in female reproductive performance by regulating the proliferation, differentiation, and apoptosis of granulosa cells in the follicle, ultimately affecting folliculogenesis and oocyte maturation (Zhang and Liu, 2002; Sun et al., 2016; Huang et al., 2022). The MAPK pathway also plays a role in regulating luteinizing hormone secretion, which stimulates ovulation and formation of the corpus luteum (Kahnamouyi et al., 2018). Additionally, MAPK signaling has been implicated in regulating the menstrual cycle and endometrial function (Zhou et al., 2010; Makieva et al., 2018). Meanwhile, previous studies found that abnormal MAPK signaling can cause reproductive disorders (e.g., infertility and embryonic death) in swine (Prochazka et al., 2012; Prochazka and Nemcova, 2019) and cattle (Sigdel et al., 2021; Tahir et al., 2021). The detection of the MAPK signaling pathway in this study may indicate that AD infection may lead to the disorder of the MAPK signaling pathway, therefore influencing female reproductive performance.

A total of 16 genes, ARHGAP19, COL14A1, DEPTOR, EXOSC1, FAM135B, FRAT1, FRAT2, LOC122905718, MMS19, MRPL13, PGAM1, PTCHD4, RRP12, TBX18, UBTD1, and ZDHHC16, were common to all three studied traits. From these, five genes were related to growth in livestock species in previous studies. For example, the gene ARHGAP19 was found to be related to body weight in yak (Jiang et al., 2022); the gene FAM135B was related to body weight growth in cattle (Serão et al., 2013; Seabury et al., 2017), the genes COL14A1 (Cardoso et al., 2018) and PTCHD4 (Doyle et al., 2020) were found to play important roles in muscle development in cattle and swine, the gene EXOSC1 (Dall'Olio et al., 2020; Ropka-Molik et al., 2018) has been related to muscle growth, and PGAM1 was found to relate to the development of adipose tissue (Xing et al., 2019). Meanwhile, two genes, FAM135B (Serão et al., 2013; Seabury et al., 2017) and COL14A1 (de Lima et al., 2020), were also found to be related to feed efficiency in cattle. In addition, several genes were found to be related to reproduction in previous studies. UBTD1 (Kongmanas et al., 2015), ZDHHC16 (Uzbekova et al., 2021; Caetano et al., 2023), RRP12 (Tiensuu et al., 2019), MMS19 (Tsai et al., 2017), and PGAM1 (Zhang et al., 2015) genes were found to be associated with the development of germ cells. The genes FAM135B and FRAT1 were detected to be associated with the reproductive performance in swine (Zhang Z. et al., 2019) and cows (Melo et al., 2017), respectively. The Tbx18 gene was related to mouse embryonic development (Wehn and Chapman, 2010).

Two GO terms, olfactory receptor activity (GO:0004984) and detection of chemical stimulus involved in sensory perception of smell (GO:0050911), were significant (q < 0.05) among all three studied traits. This may indicate that AD may influence the sense of smell in mink, although the relationship between AD and smell has not been reported in the literature. However, reduced appetite of infected mink seems to corroborate the loss of their sense of smell because the smell is vital for mink feeding behavior (Saunders, 1988). In American mink, Adney et al. (2022) speculated that individuals experimentally infected with SARS-CoV-2 may have altered sense of smell because they observed neutrophilic infiltrate in the olfactory epithelium. Thus, future studies could assess the condition of the olfactory epithelium of AD-infected mink to determine if infection could influence their sense of smell.