Understanding mechanisms of vitiligo development in Smyth line of chickens by transcriptomic microarray analysis of evolving autoimmune lesions

Background The Smyth line (SL) of chicken is an excellent avian model for human autoimmune vitiligo. The etiology of vitiligo is complicated and far from clear. In order to better understand critical components leading to vitiligo development, cDNA microarray technology was used to compare gene expression profiles in the target tissue (the growing feather) of SL chickens at different vitiligo (SLV) states. Results Compared to the reference sample, which was from Brown line chickens (the parental control), 395, 522, 524 and 526 out of the 44 k genes were differentially expressed (DE) (P ≤ 0.05) in feather samples collected from SL chickens that never developed SLV (NV), from SLV chickens prior to SLV onset (EV), during active loss of pigmentation (AV), and after complete loss of melanocytes (CV). Comparisons of gene expression levels within SL samples (NV, EV, AV and CV) revealed 206 DE genes, which could be categorized into immune system-, melanocyte-, stress-, and apoptosis-related genes based on the biological functions of their corresponding proteins. The autoimmune nature of SLV was supported by predominant presence of immune system related DE genes and their remarkably elevated expression in AV samples compared to NV, EV and/or CV samples. Melanocyte loss was confirmed by decreased expression of genes for melanocyte related proteins in AV and CV samples compared to NV and EV samples. In addition, SLV development was also accompanied by altered expression of genes associated with disturbed redox status and apoptosis. Ingenuity Pathway Analysis of DE genes provided functional interpretations involving but not limited to innate and adaptive immune response, oxidative stress and cell death. Conclusions The microarray results provided comprehensive information at the transcriptome level supporting the multifactorial etiology of vitiligo, where together with apparent inflammatory/innate immune activity and oxidative stress, the adaptive immune response plays a predominant role in melanocyte loss.


Background
Vitiligo is a postnatal hypopigmentary disease characterized by appearance of white patches due to loss of functioning pigment producing cells (melanocytes) in affected skin. Vitiligo renders patients more susceptible to various forms of skin cancer and has negative psycho-social effects on life quality of vitiligo patients, especially when the depigmented skin is present in sun-exposed areas [1,2]. The etiology of vitiligo is multifaceted, involving genetic, immune, metabolic and environmental factors and the loss of melanocytes in affected tissues is generally considered to result from autoimmune destruction of melanocytes [3].
The Smyth line (SL) of chicken is an excellent animal model for human autoimmune vitiligo, for which many phenotypical and etiological similarities to human vitiligo have been documented [4][5][6]. The genetically vitiligosusceptible SL chickens [7][8][9] develop postnatal and spontaneous autoimmune loss of feather pigmentation (SLV), with 80-95% of hatch-mates expressing the disorder during early adulthood [5,8]. Etiopathologically, SLV development and progression is associated with mononuclear leukocyte infiltration into the target tissue [6], melanocytespecific humoral [10][11][12] and cell-mediated immunity [13][14][15][16][17] and elevated levels of oxidative stress [18]. In addition, routine vaccination at hatch with live herpesvirus of turkey against lymphoma-causing Marek's disease virus was identified as a dependable environmental factor related to SLV expression [19]. Taken together, like in humans, SLV is a multifactorial disease, which underlines the suitability of the SL chicken as the avian model for human vitiligo.
Several vitiligo-susceptibility genes have been identified for human vitiligo using the candidate gene approach, the genome-wide linkage approach and the gene expression approach [20]. However, identification of susceptibility genes in humans was complicated by the tremendous difference in genetic background among humans [21]. By contrast, SL chickens are highly inbred [22] and share the same MHC-haplotype and background genes as the Brown line (BL) of chicken from which the SL was derived. Together the vitiligo-expressing SL chicken and its parental BL control provide an excellent model to focus research efforts on vitiligo-related aspects and activities. Additional attributes of the SL chicken model that make it highly suitable for vitiligo study are the predictably high incidence of vitiligo expression (80-95%) [5] and the fact that melanocytes are located in growing feathers [8]. Collection of growing feathers is a minimally invasive procedure and growing feathers regenerate, providing opportunity to repeatedly sample and study the spontaneously evolving lesions in the same individual throughout the natural process of vitiligo development.
Using the SL chicken vitiligo model, the purpose of this study was to examine gene-expression associated with vitiligo development by comparing transcriptomic profiles in samples collected from SL chickens that never developed SLV (NV) and those collected from the same SLV chickens before vitiligo onset (EV), during active (AV) and after complete loss of pigmentation (CV). The results from the current study will further shed light on the critical events related to vitiligo development in SL chickens and in humans.

Animals
Eighteen SL and four MHC (B 101/101 )-and age-matched BL chicks were randomly selected from 15 SL and 4 BL families. The BL chicken is the parental control from which the SL chicken was derived. This line continues to have a rare incidence of vitiligo (< 2%). One-day old chicks were vaccinated with live herpesvirus of turkey (Fort Dodge Animal Health, Fort Dodge, Iowa) to protect against Marek's disease. All chicks were kept in one floor pen at the Arkansas Experiment Station Poultry Farm in Fayetteville, AR, USA with feed and water provided ad libitum. The animal use was approved by the University of Arkansas Institutional Animal Care and Use Committee.

Feather collection
The living part of growing feathers (feather tip) [6] was collected twice a week starting from 4 weeks of age until the chickens were 20 weeks old (young adulthood; natural endpoint when feather growth is complete). For each chicken and collection, three feather tips were snap frozen with liquid nitrogen in Tissue-Tek ® O.C.T freezing medium (Sakura Finetek Inc., Torrance, CA). All samples were kept at -80°C until use. For each chicken, vitiligo development was evaluated at each feather collection based on pigmentation loss in the growing feathers using a scale of 1 to 5 following previously described criteria [16].

Feather sample selection for the microarray study
By the end of the 20-week feather collection, 12 of the SL chickens had developed SLV (SLV chickens) while the other 6 SL (non-vitiliginous SL chickens) and the 4 BL chickens had not developed vitiligo. To represent various stages of SLV development (early, active, and complete vitiligo) feather tips from each of the 12 SLV chickens were selected as follows for RNA isolation and down-stream analysis: two collections before (within 1 week) visible vitiligo onset (early vitiligo; EV sample; vitiligo score 1), two collections during active vitiligo (AV sample; vitiligo score 2-4; partially depigmented regenerating feathers were selected), and two collections at least one week after complete depigmentation (CV sample; vitiligo score 5). For the 6 non-vitiliginous SL (NV sample, vitiligo score 1) and the 4 non-vitiliginous BL control chickens (BL sample), feather tips obtained at the same age when the selected SLV samples were collected (age-matched) were chosen for RNA isolation.

RNA isolation and quantification
Feather tips were taken out of O. C. T. blocks and homogenized by Tissue Tearor™ (Biospec Products, Inc. Model: 985370-395) at -20°C in lysis buffer provided by Qiagen RNeasy ® Mini kit (Qiagen Inc., Valencia, CA). Total RNA was isolated from homogenates following manufacturer's instruction with on-column DNA digestion (Qiagen Inc., Valencia, CA). RNA was eluted in 30 μl RNase-, DNase-free water and stored at -80°C until use. RNA integrity and concentration were determined with an Experion™ automated electrophoresis system and Experion™ RNA StdSens analysis kit (Bio-Rad laboratories, Hercules, CA) following manufacturer's instruction.

Making RNA pools for microarray analysis
For microarray analysis, separate pools of RNA were prepared for NV, EV, AV, CV and BL samples by mixing equal amounts of RNA from the pre-selected feather tips of each bird in the respective samples. For this, the 12 vitiliginous and 6 non-vitiliginous SL chickens were evenly divided into 3 groups of 4 and 2 birds each, representing three SLV-and three NV-replications, respectively. For each SLV group, RNA pools of their individual EV, AV and CV samples were prepared. Using RNA isolated from the age-matched feather samples, one RNA pool was made for each NV group and one RNA pool was made for the BL control group. The final concentrations of RNA pools were determined by the Thermo Scientific Nano-Drop™ 1000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA).

Two-color microarray
A two-color labeling microarray system together with chicken 4 × 44 k Agilent gene expression microarray slides (array IDs: 251506810338, 251506810339 and 251506810340; Agilent Technologies, Inc., Santa Clara, CA) was employed to compare global transcriptomic profiles between NV, EV, AV and CV samples by using the BL sample as the reference (reference design). RNA pools from SL chickens (NV, EV, AV and CV) were labeled with Cy5 (fluoresces red) and the BL pool was labeled with Cy3 (fluoresces green) by incorporating cyanine 5-and cyanine 3-CTP, respectively, during cRNA generation using T7 RNA polymerase. The two-color microarray workflow from sample amplification and labeling to hybridization was carried out as previously described [23] Microarray data collection and analysis Raw data were collected and analyzed with a GenePix 4000B scanner and GenePix Pro 6.1 Software (Molecular Devices, Inc., Sunnyvale, CA), respectively. Backgroundcorrected red and green intensities for each spot were used in the subsequent analysis. Locally weighted linear regression (LOWESS) was employed to normalize the fluorescent intensities to retain only biological differences by removing undesirable systemic variations in microarray experiments. Moderated t-statistic based on the empirical Bayes method [24] was applied to each gene to identify differentially expressed genes in NV, EV, AV and CV samples relative to BL controls. Comparisons between NV, EV, AV and CV samples (within SL comparisons) were performed by One-way ANOVA test using JMP Genomics 4 (multiple means comparison) (SAS Institute Inc, Cary, NC). Genes with P ≤ 0.05 were considered differentially expressed (DE). Log2 data obtained by gene-expression analyses were transformed and reported as linear fold change values in the result section.

Microarray validation with qRT-PCR
Eleven genes were selected for validation by Taqman ® qRT-PCR using the same RNA samples as in the microarray. Primers and probes were designed by Primer Express 3.0 (Applied Biosystems, Foster City, CA) and all primer and probe sets are available upon request. RNA (2 μg) from NV, EV, AV, CV, and BL pools was reversetranscribed into cDNA in a 30 μl reaction volume using a High Capacity cDNA reverse transcription kit without RNase inhibitor according to the manufacturer's protocol (Applied Biosystems, Foster City, CA). The resulting cDNA was diluted with RNase-and DNase-free water at a ratio of 1:3. One μl of this product was subjected to PCR using Taqman 2X Universal PCR master mix with glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as the endogenous control gene and the cDNA from the BL RNA pool as the calibrator. The following PCR conditions were used: 1 cycle at 50°C for 2 minutes, 95°C for 10 minutes, and 45 cycles of denaturation at 95°C for 15 seconds, and annealing/extension at 60°C for 60 seconds. The relative gene expression was determined by the delta delta Ct (ΔΔCt) method [25].

Results
Differentially expressed (DE, P ≤ 0.05) genes in SL samples (NV, EV, AV and CV) relative to BL control samples Relative to the BL sample, 395, 522, 524 and 526 DE genes were identified for NV, EV, AV and CV samples, respectively, which gave rise to 151, 213, 212 and 212 DE genes for Ingenuity Pathways Analysis (IPA) after deletion of genes without information, genes from cDNA libraries and clones and genes related to hypothetical proteins. The top 10% of up-and down-regulated DE genes for individual SLV states are summarized in Table 1.
The magnitude of the fold change of the top 10% upregulated genes was lowest in NV samples, intermediate in EV and CV and highest in AV samples. The three upregulated DE genes with the highest expression were immunoglobulin J chain (IGJ), acidic chitinase (CHIA) and T cell receptor (TCR) delta in NV; chemokine ligand 13 (CXCL13), guanylate binding protein (GBP) and lipopolysaccharide-induced TNF factor (LITAF) in EV; sodium channel alpha subunit, IGJ and CXCL13 in AV; and, IGJ, POU class 2 associating factor 1 (POU2AF1) and myeloid antimicrobial peptide 27 in CV. The magnitude of the fold change of down-regulated DE genes was similar in NV, EV, and AV samples, but higher in CV samples due to marked decreased expression of solute carrier family 24 member 5 (SLC24A5), tyrosinase (TYR), and matrix metalloproteinase 115 (MMP115). Some of the top 10% DE genes were shared by more than two SL samples, although their expression-levels varied. Up-regulated DE genes shared by EV, AV and CV samples (SLV samples) included LITAF, tumor necrosis factor superfamily (TNFSF) 13B, interleukin 21 receptor (IL21R), complement component 3a receptor 1 (C3AR1) and sodium channel alpha subunit, with higher expression levels observed in AV than in EV and CV samples. Down-regulation was noticed in lipoprotein The 88 DE genes could be roughly divided into 4 functional groups which were immunity-, melanocyte-, stressrelated and others (Table 2) based on functions of their corresponding proteins. More than half of the DE genes were immune system related, including molecules of both innate and adaptive branches of the immune system with more DE genes relating to innate than adaptive immune response activities ( Table 2). The group of DE genes related to stress contained DE genes associated with cellular stress and apoptosis. Compared to SLV samples, the majority of DE genes in NV samples had the lowest expression level except for all DE genes related to melanocyte functions (MMP115, TRP1, TYR and SLC24A5, SLC24A2, GRP143, V-ATPase C2 subunit and Shikimate 5-dehydrogenase). Additionally, DE genes related to stress (NPY and CRYAB) and intra/ inter-cellular transport [gap junction protein alpha 5 (GJA5), transmembrane protein 9 (TMEM9) and synaptotagmin 12 (SYT12)], exhibited the highest expression in NV samples (Table 2). For most DE genes, the expression levels in EV and CV samples were intermediate to their expression in NV and AV samples. All genes related to immune system activities had the highest expression levels in AV compared to NV, EV and CV samples ( Table 2). The expression levels of DE genes in CV samples tended to be comparable to EV and/or NV samples. Expression of melanocyte-related DE genes (i.e. MMP115, TYR, TRP1, SLC24A2, SLC24A5, V-ATPase C2 subunit), however, was decreased significantly in AV and reached the lowest expression in CV samples compared to those from NV and EV samples. DE genes in the others category, e.g. intra/inter-cellular transport (GJA5 and TMEM9), also exhibited the lowest expression levels in CV compared to NV, EV and AV samples ( Table 2).

Ingenuity Pathway Analysis (IPA) for SL samples
The three major output components provided by Ingenuity Pathway Analysis (IPA) are function, network and pathway interpretations (IPA, Version 9.0; Ingenuity Systems ® ; http://www.ingenuity.com). The functional analysis identifies biological functions and/or diseases that are most closely related to the data set. Network interpretations are algorithmically generated by overlying genes from a data set onto a global molecular network developed from information contained in Ingenuity's Knowledge Base. Canonical pathway analysis identifies the pathways from the Ingenuity Pathways Analysis library of canonical pathways that most significantly fit the data set (IPA, Version 9.0; Ingenuity Systems ® ; http://www.ingenuity.com)

Functional analysis
IPA analysis provided similar composition and ranking of functions for all SL samples ( Figure 1). Common functions not only included those related to normal biological functions, such as cellular-function and -development and lipid metabolism, but also those associated with a variety of diseases, e.g. immunological disease, neurological disease, dermatological disease/condition and cancer. Among them, the inflammatory response was the most significant (lowest P values which were due to chance alone) functional interpretation and the first function identified for all SL groups. However, for any particular common function identified in the IPA analysis of SL samples, the associations were made with higher confidence (lower P values) for DE genes in SLV (EV, AV and CV) samples than in NV samples. Unique functional interpretations of NV DE genes included digestive system development and function as well as organismal survival. Based on function analysis of DE genes identified by within SL comparisons, amino acid/carbohydrate metabolism, energy production, protein folding and developmental disorder were uniquely demonstrated, however with very low confidence.

Network generation
Network output provided number wise more networks for DE genes in SLV samples than for DE genes in NV samples or DE genes identified by within SL comparison. Functional analysis of networks demonstrated that functions for the NV samples and the within SL comparison group were part of the functional networks identified for DE genes in SLV samples. In addition to shared functions between all groups, humoral immune response and immunological disease were only present in EV samples, and cell-mediated immune response was unique to CV samples. The 1 st network for NV samples was associated with functions of cellular movement, immune cell trafficking  *: NV indicates that growing feather samples were from SL chickens that never developed vitiligo; EV, AV, and CV indicate that growing feather samples were collected from SL chickens with vitiligo (SLV) within 1 week before SLV onset, during active depigmentation and at least one week after complete loss of pigmentation, respectively. Normalized microarray data were filtered to include only fluorescence values above 100 and analyzed by JMP genomics 4 to determine differences in NV, EV, AV and CV gene expression levels. Values in the body of the table were mean expression (n = 3) relative to BL samples, where numbers > 1 and < 1 represents up-regulated and down-regulated expression, respectively. Means without a common letter were significantly different at P ≤ 0.05 and cell morphology; for EV samples, functions of inflammatory response, antimicrobial response and humoral immune response; and for AV and CV samples, inflammatory response, cellular growth and proliferation and hematological system development and function, and cellular development. These number one networks are illustrated in Figure 2.

Canonical pathway analysis
Pathway composition was similar for DE genes in SL samples, with higher levels of similarity between DE genes for SLV samples. For a specific common pathway higher significance [suggested by a higher ratio (number of molecules from the data set divided by the total number of molecules in the pathway obtained from the Ingenuity's Knowledge Base) and a lower P value] was exhibited for SLV samples than for NV samples. Pathway analysis of the DE data-set generated by within SL comparison had numerically the most pathways which included all pathway components for all SL samples. The pathway of communication between innate and adaptive immune cells was most significant, constituting the number-one pathway for all SL samples and within SL comparison.

Microarray validation by qRT-PCR
Microarray data were validated by qRT-PCR based on comparable gene expression levels and similar expression trends of selected targets throughout SLV development (Table 3). Fold changes in complement component receptor 2 (CR2) and POU2AF1 from qRT-PCR, however, were noticeable higher than those from microarray.

Discussion
The Smyth line (SL) of chicken is an excellent avian model for human autoimmune vitiligo [4][5][6], which is a multifactorial disease of complex etiology resulting in loss of melanocytes in the skin. In the current study, we took advantage of the unique attributes of the SL chicken model (outlined in the Introduction) and the power of gene expression analysis by microarray to investigate the pathomechanism of vitiligo in SL chickens (SLV). To our best knowledge, this is the first paper using microarray analysis to study transcriptomic expression profiles in evolving vitiligo lesions. Expression levels of genes from this microarray study were validated based on their comparability to those from qRT-PCR (Table 3), which was suggested to have higher detection sensitivity than microarray [26]. While higher expression was observed for CR2 and POU2AF1 by qRT-PCR than with microarray, their expression trends were the same with both methodologies, supporting the validity of the microarray study. The expression difference may be due to the upper fluorescent detection limit of the microarray scanner. Differentially expressed (DE) genes identified by microarray analysis and the interpretations by Ingenuity Pathway Analysis http://www.ingenuity.com support the multifactorial nature of vitiligo. For the sake of simplicity, the discussion will focus on observations regarding the involvement of the immune system, melanocyte biology, redox status, neurology, apoptosis and other physiopathological factors identified to play a role in the spontaneous development of SLV.

Immune system
The autoinflammatory/autoimmune nature of SLV was supported by expression profiling (Table 1 and 2) and IPA functional interpretations of DE genes (e.g. inflammatory response and immunological disease) (Figure 1  and 2). DE genes of both T and B cell markers (Table 1 and 2) and functions of humoral and cell-mediated immune response (Figure 1 and 2) in the current study are in line with previous studies demonstrating involvement of both types of immune responses in melanocyte loss in the etiology of SLV [6,10,12,16]. Furthermore, it appears that the humoral immune response was particularly important in early SLV development, since DE genes in EV samples were significantly (1 st network) associated with this function (Figure 2). In addition, the humoral immune response also appeared important in SLV progression (AV) and complete depigmentation (CV) as supported by significantly up-regulated DE genes related to B cell development and/or activation such as IGJ (immunoglobulin J chain), POU2AF1 (a protein essential for the B cell response and germinal center formation), BTK (Bruton's tyrosine kinase), mucolipin 2 (a target for BTK in B cells), SPI1 (a lymphoid specific enhancer involved in B cell differentiation and activation), and PIK3AP1 (B cell activation) in these samples (Table 1 and 2). On the other hand, DE genes associated with the cell-mediated immune response were identified only in CV samples and with lower significance than DE genes associated with the humoral immune response in EV samples. A more relevant role of the humoral immune response in the pathomechanism of SLV is contradictory to the demonstrated more direct role of the cell-mediated immune response in this disorder [6,16]. Confirmation of this finding will definitely impact future direction of vitiligo research. Aside from involvement of adaptive immunity in SLV, it is apparent that innate immunity also plays an important role in SLV development since the majority of immunerelated DE genes were derived from this branch of the immune system (Table 1 and 2). Innate immunity has not received much attention in vitiligo studies in either the SL animal model or human patients. Elevated expression of genes for chemokines (e.g. CXCL13, CCL19), pattern recognition receptors (e.g. MMR, TLR15, MD-2 and PRAM-1), TNF and type 1 and 2 interferon cytokines (e.g. LITAF, TNFSF13B, ISG12-2, ITIF5, IRF1 and IRF8) ( Table  1 and 2) suggested both anti-bacterial and anti-viral activity especially in early and active stages of SLV (EV and AV). The complement system also appears to participate in SLV development and progression as supported by significantly increased expression of genes for C3, C3AR1 and CR2 in AV samples (Table 1 and 2). Functional analysis of networks by IPA revealed that DE genes identified in EV and AV samples were significantly associated (1 st and 3 rd network) with the functions of anti-microbial response and infectious disease, respectively. These observations suggest a potential role of the innate immune response as a precipitating and effector factor in the development of the melanocyte-specific adaptive immune response in SLV. Initiation of antimicrobial activity in SLV is in agreement with the strong association of live herpesvirus of turkey administration at hatch (routine vaccination to protect SL and BL chickens from Marek's disease virus) as a reliable environmental trigger of SLV expression in vitiligosusceptible SL chickens [19]. Additionally, as innate immunity is also initiated in response to cellular stress, the observed alterations in melanocyte biology and redox status (see below) may also be responsible for activation of innate immune responses in SLV [27]. Figure 1 Graphical demonstration of associated functions from Ingenuity Pathway Analysis (IPA) of differentially expressed (P ≤ 0.05; DE) genes for NV, EV, AV, CV relative to BL samples. The functional analysis identifies biological functions and/or diseases that are most significant to the data set. Molecules from the DE gene dataset that are associated with biological functions and/or diseases in Ingenuity's Knowledge Base were considered for the analysis. Right-tailed Fisher's exact test was used to calculate a P-value determining the probability that each biological function and/or disease assigned to that data set is due to chance alone. The y-axis displays the functional categories that are identified in analyses. The x-axis demonstrates the significance which is the value of -log (P). Functions are listed from most significant to least and the orange vertical line denotes the cutoff for significance (P-value of 0.05). For each analysis, only the top 18 functional categories are displayed due to large size of the data files. Figure 2 Network #1 obtained by Ingenuity Pathway Analysis (IPA) of differentially expressed (P ≤ 0.05; DE) genes for NV, EV, AV, CV relative to BL samples. DE genes from the data sets were overlaid onto a global molecular network developed from information contained in Ingenuity's Knowledge Base to algorithmically generate networks that graphically represent relationships between molecules. Molecules are represented as nodes, and the biological relationship between two nodes is represented as an edge (line). The intensity of the node color indicates the degree of up-(red) or down-(green) regulation. Nodes are displayed using various shapes that represent the functional class of the gene product. Edges are displayed with various labels that describe the nature of the relationship between the nodes.

Melanocyte biology
It was previously shown that intrinsic abnormalities of melanocytes, the target cell in vitiligo, were present in SL chickens compared to BL and LBL (vitiligo-resistant, similar plumage color pattern control) controls [7,28]. In the current study, melanocyte loss appeared to be preceded by disturbances in melanosome structure, intracellular molecule transport and tyrosine synthesis as demonstrated by decreased expression of MMP115, SLC24A2, SLC24A5, GPR143, and Shikimate 5-dehydrogenase in EV samples (Table 1). Dysfunction and loss of melanogenesis became obvious in AV and CV samples, as suggested by significantly decreased gene expression in TYR and TRP1 in these samples compared to NV and EV samples ( Table 2). This finding is consistent with findings in human studies, where down-regulation of TYR and TRP1 was revealed in lesional skin compared to non-lesional skin of vitiligo patients and skin from healthy donors [29][30][31].

Redox status
Oxidative stress was previously demonstrated in vitiligo patients [32][33][34][35] and in SLV chickens [18]. In humans, oxidative stress involvement was further supported by the repigmenting effect of antioxidant substances in the treatment of vitiligo lesions [36][37][38]. In line with findings from these studies, there were DE genes for all SL samples and those generated from within SL comparison identified by IPA to be associated with the function of free radical scavenging (including functional annotation of synthesis/production of oxidative species). NV samples had similar DE gene expression for both pro-and anti-oxidative proteins relative to the BL control samples whose expression was considered as the value of 1 ( Table 2). Up-regulation of these DE genes in SLV samples, especially AV samples (Table 2) suggested association of SLV development with disturbances in redox status. Among these, remarkably elevated expression was noticed for Gp91-phox and cytochrome b ( Table 2), suggesting an important role of these proteins in changing the redox status in SLV chickens. Gp91-phox is a subunit of NADPH oxidase containing the heme binding site important in reactive oxygen species (ROS) production [39]. The source of NADPH-dependent ROS production in vitiligo development is not clear, however, after activation, phagocytes and B cells are capable of NAPDH-dependent release of ROS [40,41]. With more DE genes related to B cell than phagocyte development and/or activation in the current study (Table 1 and 2) and the greater B cell than macrophage infiltration in feather tissues at all vitiligo states reported by Shi and Erf [6], it is likely that B cells could be an important source of ROS in SLV development. As revealed by IPA, decreased transmembrane potential of mitochondria and the mitochondrial membrane was one of the many functional interpretations of DE genes in EV, AV and CV samples suggesting mitochondrial involvement in ROS production. Elevated expression of cytochrome b and acetyl-CoA carboxylase 2, which is important in fatty acid oxidation, in EV, AV and CV compared to NV samples points to mitochondrial production of ROS probably due to premature electron leakage of oxygen at the location of cytochrome b in the electron transport chain [42]. a, SL samples included NV, EV, AV and CV samples which were from SL chickens that never developed vitiligo; from vitiliginous SL chickens within 1 week before SLV onset, during active depigmentation and at least one week after complete loss of pigmentation, respectively b, Gene expression in NV, EV, AV and CV was presented by fold changes relative to expression in BL samples in microarray and in qRT-PCR, where fold-change was determined by the delta delta Ct method. The same RNA pools for NV, EV, AV, CV and BL samples used in the microarray study were reverse transcribed to cDNA, which was subjected to qPCR with GAPDH as the endogenous control gene and cDNA from BL sample as the calibrator. All values are fold change means of three replications analyzed by JMP Genomics 4 for microarray and SYSTAT for qRT-PCR data Neurology Neural involvement in vitiligo was supported by appearance of localized vitiligo lesions following nerve damage in affected skin and by morphological changes in terminal cutaneous nerves in vitiligo patients [43]. In line with observations in humans, neurological disease was presented as one of the functional interpretations by IPA for DE genes in SL samples and DE genes identified based on within SL comparisons in the current study. In addition, expression of TAC1 (gene for neurotransmitter substance P, neurokinin A, neuropeptide K and neuropeptide gamma), KCNE1L (a protein capable of regulating neurotransmitter release and neuronal excitability), ALS2 (plays a role in axon and dendrite development) and NPY (neuropeptide Y) were decreased in EV, AV, and CV samples compared to NV samples (Table 1 and 2), supporting neural involvement in SLV development. NPY expression was also investigated in human studies; the results, however, were not consistent [44,45] and differed from the decreased expression with SLV progression in the current study. Other studies revealed an increased number of nerve growth factor receptor-, calcitonin gene-related peptide-and protein gene product 9.5-positive nerve fibers in involved skin from vitiligo patients compared to healthy controls [45,46]. Moreover, increased serum and urinary levels of catecholamines and their relative metabolites were found in vitiligo patients, especially in early and active stages of vitiligo, than in controls [47][48][49]. It was speculated that a role of catecholamines in vitiligo onset was due to their potential to undergo oxidation, resulting in formation of quinones, semiquinone radicals and oxyradicals [47,49].

Apoptosis
While the mechanism by which melanocytes are destroyed in vitiligo is still in debate, it has been suggested that they disappear via apoptosis (presumably initiated by cytotoxic T cells) rather than necrosis [50].
In the current study, significant up-regulation of apoptosis-related DE genes, especially in AV compared to NV, EV and/or CV samples (Table 2), indicated a close relationship of apoptosis with SLV development. These apoptosis-related up-regulated DE genes included GZMA, LITAF (induced by p53 and has been implied in p53 induced apoptosis), caspase-associated recruitment domain, member 11 (CARD11, a positive regulator of cell apoptosis), gasdermin 1 (involved in apoptosis induction), TNF receptor associated factor 5 (TRAF5) and programmed cell death 1 (PDCD1) ( Table 2). In addition, participation of apoptosis in SLV etiology was also supported by the IPA association of functions such as cell death ( Figure 1) and pathways of apoptosis signaling and cytotoxic T cell-mediated apoptosis of target cells with DE genes in SL samples and DE genes identified by within SL comparison. These findings are consistent with elevated levels of apoptotic (TUNEL+) cells in the melanocyte region of growing feathers from vitiliginous SL chickens (particularly during active depigmentation) compared to non-vitiliginous SL, BL and LBL controls [17]. The mechanism of apoptosis was also indicated by decreased levels of anti-apoptotic protein Bcl2 in lesional skin compared to nonlesional skin of vitiligo patients [51,52].

Others
The decreased expression of the gene for keratin especially in EV and AV compared to NV samples (Table 1) may suggest disturbance in keratinocyte function in early stages of SLV development. Since keratinocytes play important roles in supporting melanocyte growth and regulating melanogenesis in melanocytes [53], keratinocyte malfunction can have detrimental effects on melanocyte functions. Constant and concurrent down-regulation in matrix metalloproteinase MMP9 and MMP13 in SL chickens could be due to aberrant keratinocyte function, based on the observations in human studies that keratinocyte-derived MMP9 was important in repigmentation by enhancing melanocyte migration [54,55]. Three DE genes TMEM9, TMEM22 and sorting Nexin 10 may suggest abnormal inter/intracellular molecule/organelle trafficking in SLV chickens. The reason for noticeably decreased expression in lipoprotein VSAF in SL chickens compared to BL is not clear due to lack of functional information on this protein. Significant down-regulation of parathyroid hormone (PTH) in EV, AV and CV compared to NV samples may suggest a possible role of PTH in SLV expression as its biological function is to increase serum calcium levels and vitamin D synthesis, both of which are important in melanogenesis.

Conclusions
The multifactorial etiology of autoimmune SLV is supported by the DE gene profiles from the microarray study and by their interpretations by IPA. Basically, in susceptible SL chickens SLV may occur primarily through melanocyte-specific immune destruction of aberrant melanocytes under certain pathobiological circumstances where pre-existing "danger" signals and oxidative stress may act as essential triggers.