Introduction

As noted by Darwin, chickens have the greatest phenotypic diversity among birds. Compared with their wild ancestor, the Red Junglefowl, domestic chickens harbor many typical domestication characteristics, such as larger body size, tame behavior, and higher egg productivity. Identifying the genetic changes underlying these phenotypic changes should provide new insights into the successful domestication of animals and further inform future breeding programs. A pioneering study by Rubin et al.1 used re-sequencing of genomes of commercial chickens as well as populations of Red Junglefowl to identify several selective sweeps associated with reproduction and growth. However, their approach had limitations in comprehensively elucidating the genetic differences that account for the phenotypic differences between domestic chickens and Red Junglefowl. This was due to, for example, the use of commercial chickens that evolved from Red Junglefowl under strong selection by humans for specific traits associated with growth (e.g., meat production) and reproduction (e.g., egg production) as wells as the use of pooled re-sequenced genomes. Despite the great leap in genetic technologies witnessed in the recent past, the genetic mechanisms underlying many interesting evolutionary differences between domestic chickens and their wild ancestors, for example, reduced vision prowess in chickens, remain unclear.

Vision is one of the most crucial abilities affecting the survival of animal species, as it influences an array of core behavioral traits associated with mating, foraging, and predator avoidance2. However, domesticated animals including dogs, horses, and chickens, which are bred and protected by humans, seem to exhibit markedly weaker visual acuity compared to their wild counterparts3, 4,5,6,7. Domestic chickens (except game fowl) share with these other domestic animals a reduced visual ability as compared to the Red Junglefowl, assumedly the cost of domestication due to the relaxation of selective constraints. The living condition changes that accompany domestication are a reflection of human interventions such as offering protection, provision of food and promotion of the animal breeding. These interventions reduce the pressure for sharp vision in the domestic chickens3,4. However, unlike other phenotypic differences like growth and reproduction1, no reports exist for the genetic mechanisms behind the change in visual acuity observed in domestic chickens.

To promote a better understanding of the underlying genes/variants responsible for phenotypic changes that occurred in the domestic chickens, specifically reduced visual capability, we generated genome sequences with high coverage (average 18.9×) from Red Junglefowl and indigenous chickens bred in villages of Yunnan province, Southwest China, one of the regions where domestic chickens originated from8. Based on these genomes, we characterized potential mechanisms underlying the phenotypic evolution of domestic chickens by comparing genome sequences and integrating transcriptome data from different sub-regions of the brain and the retina.

Results

Positively selected genes in domestic chickens compared with the Red Junglefowl

Variants/genes underlying phenotypic changes in the domestic chicken likely evolved rapidly after domestication. Based on the genomic variation data obtained in this study, we identified candidate genes that potentially underwent rapid evolution in the domestic chickens, and thus might have contributed to the phenotypic differences between domestic chickens and the Red Junglefowl. We first employed an outlier approach to identify candidate selected genes based on the empirical data. This method is free from assumptions of demographic history and is considered to be robust to the confounding effect of population demographic history. Typically, regions or loci evolving rapidly that have experienced selection would show specific signatures of variation, including high population differentiation, significantly reduced nucleotide diversity levels and long-range haplotype homozygosity9. Based on these principles, we examined four different parameters to identify footprints of artificial selection associated with the evolution and domestication of chicken (Figure 1, Supplementary information, Figure S1): FST10, nucleotide diversity (Pi), cross-population extended haplotype homozygosity (XP-EHH)11 and the cross-population composite likelihood ratio (XP-CLR)12. The first percentile rank was used as a threshold to identify candidates throughout the analysis.

Figure 1
figure 1

Analysis of the signatures of positive selection in the genome of domestic chickens. Genomic landscape of the XP-CLR values (A) and XP-EHH values (B) in the genome of domestic chickens. (C) Numbers of PSGs identified by the four methods listed in each of the Venn diagram components.

A number of protein-coding genes with significant higher FST (184 genes), XP-EHH (212 genes), XP-CLR (335 genes), and a lower value for nucleotide diversity (Pi, 284 genes) were identified, and considered to be potential candidate genes that experienced selection during domestication. Of these genes, only 36 were identified by all four methods (Figure 1C, Supplementary information, Table S12). It is not surprising that the overlap among the positively selected genes (PSGs) detected by the different statistical approaches was underwhelming. First, the distinct statistics applied to scan the genome for positive selection are based on different signatures of population variation9. Another reason might be attributable to a pitfall in the outlier approach, whereby an outlier value in one analysis, falling in the 99th percentile of the empirical distribution, may not be classified as an outlier in another study where it may fall within the 98th percentile of the study's empirical distribution13. This pattern is very common in studies on positive natural selection in humans. For example, only 14.1% loci were identified in two or more studies with large-scale genome-wide scans of PSGs13. In the current study, possibly due to the above stated factors, no significant enrichment for any functional category associated with vision was found among the candidates identified by the FST test.

Gene ontology (GO) analysis for the genes showing signatures of positive selection pointed out a significant enrichment for categories associated with the development of the nervous system (Supplementary information, Tables S8-S11). This finding is not a surprise and is consistent with the rapid evolution of genes in the nervous systems of other domestic animals, which has collectively been attributable to behavioral shifts that accompany domestication14,15,16,17,18. In addition, functional enrichment analysis of genes in regions with the top 1% XP-EHH scores revealed that many candidate genes related to the nervous system were over-represented in the following Human Phenotype Ontology (HPO) categories: abnormality of the nervous system (31 genes), abnormality of the central nervous system (30 genes), and cognitive impairment (25 genes) (P < 0.05) (Supplementary information, Table S8). After retrieving regions with the top 1% FST values, we found that 15 genes significantly over-represented in the HPO categories are related to brain disorders — abnormality of the hindbrain, abnormality of the metencephalon, and abnormality of the cerebellum (Supplementary information, Table S9). Our analysis also picked out these functional categories as regions with low nucleotide diversity (Supplementary information, Table S10). Taken together, these findings suggest that genes involved in both the function and development of the central nervous system showed strong directional selection, which may underlie some of the behavioral alterations occurring during chicken domestication.

For the top 10 SNPs with the highest XP-EHH values, only a single protein-coding gene, ATP13A4, was identified. The XP-CLR test also showed that this gene had the most significant XP-CLR values at 10-kb, 5-kb and 2-kb grid sizes (Figure 1; Supplementary information, Figure S2). ATP13A4 plays vital roles in early neuronal development19. It is highly expressed in specific areas of the brain, including the lateral inferior frontal cortex and the temporoparietal cortex, and is associated with benign epilepsy with centrotemporal spikes20. It is also associated with speech apraxia among human children21,22. These results collectively suggest that ATP13A4 likely plays a key role in behavioral evolution during chicken domestication. Similarly, the gene UNC80, which has a function in the nervous system23,24, also experienced positive selection and had the highest identified XP-EHH score by the sliding window analysis (Figure 1B).

We also identified several other PSGs that are involved in body size development and growth rate, indicative of their role in evolution of the comparatively larger body size among domestic chickens compared to the Red Junglefowl. For example, using different analyses, we identified two genes, TSHR and TBC1D1, both of which were previously reported to regulate seasonal reproduction and growth among domestic chickens1,25. IGF2BP3 was identified by its high score in the 10-kb grid size XP-CLR test (Figure 1), and further supported by the 5-kb and 2-kb grid size analyses (Supplementary information, Figure S2). Three other methods that we used in the present study further detected selection signals for this gene (Figure 1, Supplementary information, Figure S1). This is not entirely surprising, as IGF2BP3, an important member of the insulin-like growth factor binding protein family, has pro-growth functions via regulating IGFs (IGF1 and IGF2) to modulate the insulin signaling pathway26,27,28. Moreover, IGF2BP3 is located within a QTL region associated with chicken body size (http://www.animalgenome.org/cgi-bin/QTLdb/GG/index) and has differential expression in the muscles of dwarf and normal-sized chickens29.

To further confirm the findings in this study, we used the Sanger sequencing method to re-sequence and genotype 67 SNPs that cover seven candidate PSGs, TSHR, ATP13A4, RPGRIP1L, COL18A1, VIT, OPA1 and IGF2BP3, in larger chicken populations. We found that the SNPs genotyped in the larger populations still showed higher FST values, which is strongly correlated with the genomic analysis (P < 0.00001) (Supplementary information, Figure S3 and Table S5). When we compared our candidate PSGs with those identified in the study of Rubin et al.1, only a limited number of genes were shared (Supplementary information, Figure S9). Two important genes, TSHR and TBC1D1, previously reported by Rubin et al. to be involved in regulation of seasonal reproduction and growth in domestic chickens, were again identified to undergo positive selection in our study by the four methods applied.

Positive selection rather than relaxation of purifying selection of vision-related genes in domestic chickens

An additional group of genes possessing signatures of positive selection among the tested domestic chickens was found to be significantly enriched for categories associated with the development and maintenance of the eye as well as vision. Particularly, the PSGs detected by the XP-CLR and Pi tests included abnormality of the posterior segment of the eye (19 genes, HP:0004329), nystagmus (18 genes, HP:0000639), abnormality of the fundus (18 genes, HP:0001098), aplasia/hypoplasia affecting the eye (10 genes, HP:0008056), and visual loss (5 genes, HP:0000572) (full details in Table 1). Specific examples in this group are ZNF469 and RD3, with selection signals supported by both XP-CLR and Pi, which have vital roles in vision. ZNF469 encodes a zinc-finger protein and mutations in this gene are associated with brittle cornea syndrome and keratoconus, a common inherited ocular disorder resulting in progressive corneal thinning30,31. RD3 on the other hand encodes retinal degeneration 3 protein, and mutations in this gene cause Leber congenital amaurosis type 12, a disease that results in photoreceptor degeneration32.

Table 1 Gene functional enrichment categories involved in the vision function found in positively selected genes detected by the methods Pi, XP-CLR and XP-EHH, and differentially expressed genes detected by Cuffdiff in retina

Positive selection on vision-related genes in the domestic chicken was unexpected since domestic chickens, as previously reported, have a significantly weakened vision compared to their ancestor, the Red Junglefowl3,4,5. Several other studies have proposed that the weakened visual ability observed in domestic chickens is attributed to relaxation of functional constraints during domestication3,4.

Functional enrichment analysis of the 36 genes that overlapped among the candidate gene lists identified by the four different methods showed that several genes are associated with vision-related functional categories — abnormality of the eye/vision (OPA1, COL18A1, RPGRIP1L and NBAS), and visual impairment (OPA1, COL18A1 and RPGRIP1L) (Supplementary information, Table S13). COL18A1 encodes the alpha subunit of type XVIII collagen and mutations in this gene are associated with Knobloch syndrome33,34, which features high myopia, vitreoretinal degeneration with retinal detachment, and macular abnormalities in humans. Col18a1−/− knockout mice display abnormalities in their iris and ciliary bodies, and retina vessel regression. These mice also show abnormal deposits between the basal infoldings of the retinal pigment epithelium, which could result in deterioration of retinal pigment epithelium function and attenuation of visual function as well as pathological electroretinograms (ERG)35,36,37,38. OPA1, encoding a dynamin-like mitochondrial GTPase, is reported to be involved in an autosomal dominant optic atrophy, which is known as Kjer's optic atrophy. It affects retinal ganglion cells and axons forming the optic nerve and leads to progressive visual loss39,40.

Demographic history, relaxation of purifying selection and hitchhiking could confound the detection of positive selection. The outlier approach to identify candidate genes used above is based on empirical data that is free from any assumption concerning demographic history, and is considered to be robust to the confounding effects of population demographic history, where genome-wide forces affect the patterns of variation at all loci in a genome in a similar manner, whereas directional selection specifically acts on certain loci9. Notwithstanding the outlier approach already applied, we further performed coalescent simulation analysis based on simulated demographic histories (Supplementary information, Figure S4), and found that the candidate PSGs identified by the outlier approach still harbored statistically significant signals compared with simulated data (Figure 2A-2C). To exclude the possibility that hitchhiking possibly contributed to the signatures of selection on the vision-related genes, we retrieved the vision-related genes potentially undergoing positive selection and compared the pattern of variation in regions upstream and downstream of the genes to that of variation within these vision-related genes. As expected, these vision-related genes harbored higher levels of population differentiation, and higher values of XP-EHH values, than adjacent genes and regions (Figure 2D-2F).

Figure 2
figure 2

Positive selection on genes involved in vision. (A-C) Coalescent simulations (ms) of positive selection detected by FST (A), Pi (B) and XP-EHH (C). (D-F) FST(D), Pi (E) and XP-EHH (F) values for the visual-related candidate genes compared with the values for their upstream/downstream sequences with the same length as the gene (Group A), adjacent genes (Group B) and upstream/downstream 50 kb sequences (Group C). (G) The ratio of the numbers of non-synonymous SNPs to the numbers of synonymous SNPs (pN/pS). (H) the alternative allele frequency (AAF) of the visual-related candidate genes in the village chickens and Red Junglefowl.

The ratio of non-synonymous to synonymous substitution rates (pN/pS) in the coding sequences was used to assess the relaxation of purifying selection, in which case pN/pS should increase41,42. In this analysis, we retrieved SNPs in vision-related genes identified under potential positive selection by the four methods, and counted the numbers of SNPs unique to the domestic chicken or Red Junglefowl. We found lower pN/pS in domestic chicken than in the Red Junglefowl (Figure 2G). Compared with the Red Junglefowl, domestic chickens harbored a lower frequency of alternative alleles (Figure 2H). In addition, we employed different methods based on different population genetic theories to detect positive selection, which could possibly reduce the false discovery rate due to the relaxation of purifying selection and/or demographic history. Our analyses indicated a strong unlikelihood that relaxation of purifying selection occurred in the genes associated with vision.

Differential expression of PSGs in the retina between domestic chickens and Red Junglefowl

Since our results on positive selection in vision-related genes were somewhat surprising, given previous viewpoints3,4, we sought to better understand the evolution of the domestic chicken by sequencing (via RNA-seq) transcriptomes from five tissues (cerebral cortex, corpus striatum, optic lobe, cerebellar vermis, and retina) from the Red Junglefowl and domestic chickens. We calculated differences in expression levels for each gene between Red Junglefowl and domestic chickens (see Materials and Methods). This analysis showed that genes with signatures of positive selection presented significantly higher levels of expressional difference in retinal tissue between the domestic chicken and Red Junglefowl than other genes (Figure 3A-3D, FST, P = 0.075; Pi, P = 2.99 × 10−7; XP-CLR, P = 0.015; and XP-EHH, P = 5.50 × 104, Mann-Whitney U-test). We sought to verify and replicate this result by profiling another 10 transcriptomes (6 transcriptomes of 3 retina from 3 Red Junglefowls and 4 transcriptomes of 2 retina from 2 village chickens) with two replications per sample using a different sequencing platform (see Materials and Methods). The results remained synonymous to the earlier analysis when the data from two platforms were analyzed separately (Supplementary information, Figure S5A, FST, P = 0.104; Pi, P = 0.018; XP-CLR, P = 0.178; and XP-EHH, P = 0.005, Mann-Whitney U test), or jointly (Supplementary information, Figure S5B, FST, P = 0.040; Pi, P = 2.23 × 10−4; XP-CLR, P = 0.049; and XP-EHH, P = 8.52 × 10−4, Mann-Whitney U test). However, the four brain regions did not show significantly differential expression patterns except for the corpus striatum by the Pi method (Figure 3B, P = 0.036, Mann-Whitney U test). It is possible that some PSGs are upregulated while others are downregulated in the brains of the village chickens, resulting in the non-significant variations noted when all PSGs were analyzed together. In addition, a greater number of PSGs exhibited expression level changes between the Red Junglefowl and village chicken in the eye than any of the four brain regions (Supplementary information, Figure S6). This pattern supports the hypothesis that positive selection drove changes in the regulation of gene expression that attenuated vision in the domestic chickens.

Figure 3
figure 3

Gene expression for PSGs in five tissues (cerebral cortex (CC), corpus striatum (CS), cerebellar vermis (CV), optic lobe (OL), and retina) comparing domestic chickens and Red Junglefowl. (A-D) PSGs as identified by FST (A), Pi (B), XP-CLR (C) and XP-EHH (D). (E) Difference in expression for six genes between the village chicken and Red Junglefowl in the retina verified by qPCR. (F) Correlation of the gene expression measured by qPCR versus RNA-seq. Fold change, implying (FPKM + 1)VC/(FPKM + 1)RJF for each gene. Statistical significance is indicated by *, where P < 0.05, and **, where P < 0.01.

We then used Cuffdiff43 to retrieve a series of genes that exhibited significant differential expression between the Red Junglefowl and domestic chicken in the cerebral cortex (137 genes), corpus striatum (24 genes), cerebellar vermis (9 genes), optic lobe (19 genes), and the retina (57 genes). Functional enrichment analysis of the 57 differentially expressed genes in the retina between the Red Junglefowl and domestic chickens showed that four genes (RPGR, GUCA1A, TRPM1 and PDE6B) were involved in several HPO categories including decreased central vision, congenital stationary night blindness, and abnormal ERG (Table 1). We verified the differential expression of these four genes, as well as two additional genes that are important for vision, RHO44 and NR2E345,46,47, via quantitative real-time PCR (qPCR) for retina mRNA (Figure 3E). Differences in gene expression measured by qPCR and RNA-seq were significantly correlated (Figure 3F). The significant changes in the mRNA expression in the retina between Red Junglefowl and the domestic chickens potentially have functional consequences in vision due to their respective functions. For example, NR2E3 (photoreceptor cell-specific nuclear receptor, group E, member 3) is a vital transcriptional regulator essential for photoreceptor development and differentiation due to its capability as an activator of rod cell development and a repressor of cone development45,46,47. In domestic chickens, downregulation of NR2E3 expression is associated with weaker visual ability and lower optical sensitivity3. Similarly, RHO (rhodopsin), GUCA1A (guanylate cyclase activator 1A) and PDE6B (phosphodiesterase 6B) are all vital genes in the phototransduction pathway44,48,49. Among domestic chickens, expression of GUCA1A, PDE6B and RHO were downregulated in the retina, suggesting a depression in the phototransduction pathway in the domestic chicken. In addition, RPRG and TRPM1 exhibited upregulation in domestic chickens. Both genes play important roles in vision50,51, for instance, RPGR encodes a retinitis pigmentosa GTPase regulator and overexpression of a RPGR transcript leads to severe photoreceptor degeneration in mice52.

Functional analysis of VIT, a PSG presumably involved in vision

PSGs presented significantly higher levels of differential gene expression in retinal tissue between the domestic chicken and Red Junglefowl and differentially expressed genes in retinal tissue were also involved in visual ability. However, only two genes, RPGR and VIT, were identified to undergo positive selection in domestic chicken (Figure 1B; Supplementary information, Figure S1B; Figures 3E, 4A-4C). RPGR was supported with signal of positive selection by both Pi and XP-EHH tests (the first percentile rank as threshold) (Figure 1B; Supplementary information, Figure S1B) and had an increased level of expression in village chickens compared to Red Junglefowl (Figure 3E). Mutations in RPGR are associated with X-linked retinitis pigmentosa (XLRP) and result in severe and progressive loss of vision in dogs and humans50. The second gene, VIT, showed a higher level of population differentiation (FST) and lower level of nucleotide diversity (Pi) in village chickens (Figure 4A), and was found to show significant differential expression in the retinal tissue in domestic chicken compared to Red Junglefowl (Figure 4B). Results from qPCR further confirmed the downregulation of VIT in the village chicken compared with Red Junglefowl (Figure 4C). Until now, the function of VIT remains unclear. However, it demonstrated specific expression in the eye in the mice (Supplementary information, Figure S7), leading us to reason that VIT is likely involved in vision and that downregulation of VIT expression might play a role in the evolution of vision in domestic chicken.

Figure 4
figure 4

Signal of positive selection and knockdown array of VIT in zebrafish. (A) FST and Pi values showed selection signal for VIT gene. (B) RNA-seq analysis presented differential expression of VIT. (C) qPCR confirmation of the downregulation of VIT in retinal tissue between the domestic chicken and Red Junglefowl. (D) Expression analysis of VIT in mice. Expression profiling in retinal pigment epithelium (RPE) of mice exposed to 10 000 lux of cool white fluorescent light for 18 h was based on expression data of retina from Hmx1 loss-of-function mice (GSE37773) and wild-type mice (GSE47002), and expression data from photoreceptor cells of Nrl−/− mice and Crx::Nr2e3/Nrl−/− mice (GSE5338). (E-F) Lateral view of a zebrafish with measurements of body length (E) and eyeball diameter (F). Scale bar, 0.5 mm. (G-H) Body lengths (G) and eyeball diameters (H) of embryos injected with VIT-e3i3-MO (E3I3-MO) or wild-type control at 32 h post-fertilization (hpf) (VIT morphants and control embryos, n = 10 each; mean ± SEM). Body length is significantly decreased in VIT morphants; *P < 0.0001; NS, not significant.

To obtain a better understanding of the function of VIT, we first examined the VIT expression data available for the mouse. Expression of VIT was significantly upregulated in the retinal pigment epithelium of mice after exposure to 10 000 lux of cool white fluorescent light for 18 h (Figure 4D). Expression of VIT also changed after the knockout of several genes that are crucial for the development of vision. For example, VIT was upregulated in the retina of Hmx1-knockout mice (Figure 4D). Hmx1 (H6 family homeobox 1) is a transcription factor crucial for the development of the eye. Other transcriptional regulators like NRL, CRX and NR2E3, are also important for photoreceptor differentiation46. From a past study, endogenous NR2E3 transcript or protein was not detectable in the retina of Nrl−/− mice46. However, when exogenous NR2E3 was expressed in post-mitotic photoreceptor precursors using the Crx promoter in Nrl−/− mice, expression of VIT was significantly upregulated (Figure 4D). Consistently, the expression levels of both NR2E3 and VIT were lower in village chicken than in Red Junglefowl. These expression data suggest a potential function of VIT in vision.

To further study the function of VIT, we performed a knockdown array of VIT in zebrafish (Daniorerio) embryos by targeting its expression with a specific antisense morpholino (MO) oligonucleotide preventing the proper splicing of exon 3 (E3I3-MO) (Figure 4E-4H). While the body length of VIT-knockdown zebrafish is significantly decreased compared with controls (Figure 4G), no significant change in eye size was observed (Figure 4H), suggesting that the phenotypic effects of VIT knockdown are in tissues other than the eye in the zebrafish. However, since expression of VIT is highly eye-specific in mice (Supplementary information, Figure S7), downregulation of VIT might have greater influence on the eye than other tissues in species such as chicken. We then examined the effect of VIT knockdown on the general morphology of the zebrafish retina by using Tg (Ath5-gal4:UAS-GCaMP1.6) larvae, where most photoreceptors, retinal ganglion cells and some amacrine cells are labeled with GCaMP1.6 (Figure 5A). Compared with control larvae, the general structure of the retinal layers and the projection of the optic nerves were not changed in the VIT morphant (Figure 5A-5D), consistent with the above observation that there is no significant change in eye size between control and VIT morphant (Figure 4H).

Figure 5
figure 5

Electrophysiological recordings and behavioral tests showed that VIT is necessary for normal physiological function of zebrafish retina. (A, B) Confocal images showing GCaMP1.6 expression in photoreceptors (PhR), retinal ganglion cells (RGC) and amacrine cells (AC) of a 3-dpf Tg (Ath5-gal4:UAS-GCaMP1.6) larva (A) and a VIT morphant (B). The two red dash circles indicate the position of inner and outer plexiform layers. Scale bar, 100 μm. (C, D) Confocal images showing the projection of optic nerves from two eyes of a 3-dpf wild-type larva (C) and a VIT morphant (D). Scale bar, 100 μm. (E, G) ERG recording averaged from ten responses to a 2-s light ON stimulus in 3-dpf (E) or 5-dpf (G) wild-type larvae and VIT morphant. (F, H) Summary of ERG recording data at 3 (F) or 5 dpf (H). (I, J) Examples showing a looming stimulus-induced escape behavior in a 5-dpf wild-type larva (I), but not in a 5-dpf VIT morphant (J). The red arrows indicate the body axis of the larvae before and after looming-induced C-startle escape. Scale bar, 1 cm. (K) Summary of data showing the probability of looming-induced escape behaviors in 5-dpf wild-type larvae and VIT morphant. The number in the inset of F, H and K represents the number of larvae examined. **P < 0.01, ***P < 0.001; two-tailed unpaired Student's t-test for the data in F, H and K. Data are represented as mean ± SEM.

To determine whether the knockdown of VIT gene leads to defective retinal function, we performed electrophysiological recordings and behavioral tests. Application of a 2-s whole-field flash near the eyes of zebrafish elicited a clear ERG response from the eye, where b and d waves occurred at the onset and the offset of the flash, respectively (Figure 5E). Knockdown of VIT dramatically decreased the amplitude of b and d waves in both 3-dpf and 5-dpf larvae (Figure 5E-5H). Since b and d waves are generated mainly by ON and OFF bipolar cells, respectively53, this result indicates that VIT might be necessary for the normal function of bipolar cells. Furthermore, we examined the effect of VIT knockdown on visual-induced escape behaviors. A looming stimulus, an expanding black disc, was applied to mimic the approaching of a predator in the natural environment. Wild-type larva exhibited a fast C-shape escape response to an approaching looming stimulus (Figure 5I), but this visual behavior was largely impaired in the VIT morphant (Figure 5J and 5K).

Apoptosis plays an important role in embryogenesis and tissue homeostasis54. The early steps in eye development involve extensive cell death associated with morphogenesis, and the suppression of apoptosis in cells of the lens lineage by fibroblast growth factors is an important component of lens morphogenesis55. We therefore examined the role of VIT in apoptosis in developing eyes by staining zebrafish embryos with acridine orange (AO). We detected few apoptotic cells in wild-type control zebrafish. In contrast, significantly increased levels of staining were observed throughout the eyes of zebrafish injected with 4 ng of VIT-e3i3-MO (the antisense MO oligonucleotide), suggesting increased levels of apoptosis (Figure 6). In addition, we detected increased levels of apoptosis in other tissues including brain and heart (Figure 6, Supplementary information, Figure S8).

Figure 6
figure 6

Morpholino knockdown of VIT induces potent apoptosis in the eye and brain. Wild-type control embryos and embryos injected with VIT-e3i3-MO (E3I3-MO) were stained with acridine orange (AO) at 32 hpf. Apoptotic cells are visible as bright green spots or black spots, and less bright homogenous green or black staining, an unspecific background staining. (A-D) Uninjected wild-type control zebrafish exhibited few or no apoptotic cells in whole organism. In contrast, significantly increased staining was observed throughout the brain and eye in zebrafish injected with 4 ng VIT-e3i3-MO (E-H, red arrows). Lateral view, anterior, left (A-H). (I) Quantification of apoptosis particle number in brain shows a 10-fold increase in VIT morphants (n = 10) at 32 hpf. (J) Quantification of apoptosis particle number in eye shows a 9-fold increase in VIT morphants (n = 10) at 32 hpf.

All these analyses collectively indicate that the VIT gene is necessary for normal physiological function of zebrafish retina, and the evolutionary changes in VIT might play a vital role in the evolution of vision in the domestic chicken.

Discussion

By an extensive analysis of the genome sequences of Red Junglefowl and village chickens, we found that a series of vision-related genes have undergone positive selection rather than relaxation of purifying selection in domestic chickens. In particular, VIT likely has a function in vision and the evolution of VIT might play a role in the evolution of vision in domestic chicken. A caveat of our approach is that the individuals studied possibly do not represent the entirety of genetic variation seen among chickens.

Functional enrichment analysis showed that a group of genes identified by several approaches based on the genomic and transcriptomic data were involved in the vision-related categories. We compared our result with a previous report by Rubin et al.1. Only limited number of genes overlapped (Supplementary information, Figure S9), among which were two important domestication-related genes, TSHR and TBC1D1. The limited overlap between our report and study by Rubin et al.1 might be due to differences in sequencing and statistical approaches or the use of different chicken breeds. The complicated history of domestic chicken could also contribute to the observed differences56. However, functional enrichment analyses of PSGs reported by Rubin and colleagues showed that many of the genes identified in all domestic lines (AD), commercial broiler lines (CB) as well as in layer lines were functionally involved in vision-related categories (Supplementary information, Table S14-S16), which is synonymous with our results. Since the precise nature of the variations that were selected in our study remain unclear, additional functional, biotechnological and physiological experiments will be needed to identify the mechanisms that explain the attenuated vision in domestic chickens.

Using the parsimony principle, the comparatively weaker eyesight in domestic chickens has been explained as a consequence of relaxation of the functional constraints, as domestic chickens are far less dependent on vision for their behaviors like foraging, avoiding predators and mate choice. They live in a predictable farm environment compared to the Red Junglefowl, which does not require acute vision for survival. However, our present study identified many candidate genes involved in the development of vision with signatures of positive selection, including lower nucleotide diversity and long range of haplotype homozygosity, rather than signatures of relaxation of purifying selection. Moreover, the PSGs show significant changes in their expression in the retina between village chickens and Red Junglefowl. Weakened visual ability in domestic chickens as a by-product of artificial selection on other traits, such as large body size, could be another possibility. Larger body size, for instance, could induce larger eyes that may lose optical sensitivity3. Nevertheless, it is implausible that the large number of genes involved in vision with signatures of positive selection identified in our study (Table 1) could be due to hitchhiking with other genes involved in other selected traits, such as body size. First, vision-related genes are distributed on multiple chromosomes, and thus each gene would have to hitchhike. Second, the levels of population differentiation and XP-EHH values within each of the vision-related genes were higher than those for the adjacent regions/genes, showing that hitchhiking could not generate the signatures of selection seen in these vision-related genes (Figure 2).

It is evident that hitchhiking, relaxation of selection or demographic history cannot offer a concrete explanation to our findings. However, nearly 150 years ago, Darwin mentioned a kind of selection that may be unconscious57. Given the results of this study, we propose that during domestication, this unconscious selection may have driven the weakening of vision among domestic chickens. Birds gather visual information through a scanning behavior (also known as vigilance) to make decisions relevant for survival (e.g., detecting predators, making mating choices and finding food)58. Presumably, the progenitors of domestic chickens that harbored weaker vision may also have possessed a reduced fear response hence vigilance, allowing them to be more easily caught, reared, and domestic by humans. Chickens with higher constant vigilance would not be easily kept in a cage, particularly in the evening, and thus might not have been selected for domestication. In addition, growth and egg production improves with blindness in chickens. Blind chickens tend to thrive better under more crowded conditions than their sighted counterparts59. It is therefore plausible to reason that selection on behavior — which is well documented for many domestic species17 — unknowingly promoted the selection for chickens with weakened vision.

Dogs also harbor weaker visual acuity compared to their wild ancestors, however, in contrast to chickens, there was no evidence for positive selection in vision-related genes in previous studies14,15,16. Vision is important in birds for survival58, but other sensory abilities including hearing and smell, might play more important roles in mammals. In dogs, many hypotheses have been proposed to explain the tameness necessary for successful domestication17. In chickens, besides behavioral evolution, the evolution of vision might have also been important for successful domestication. Clearly, more targeted studies should clarify this hypothesis, but at the moment this explanation offers a compelling model that may shift some of our perspectives on the genetic history of human-driven domestication of key animals, and raise the appreciation of unconscious selection in the characterization of the evolution of domestic animals60.

Materials and Methods

This study was approved by the Ethics and Experimental Animal Committee of Kunming Institute of Zoology, Chinese Academy of Science, China.

Genomic and RNA-seq data collection

Genome sequences of domestic chickens from the area of origin might be more suitable for revealing the phenotypic changes after early domestication from the Red Junglefowls (RJF), than those from commercial layer and broiler chickens bred during the 20th century in Europe. To identify potential underlying gene/variants responsible for phenotypic changes in domestic chickens, we used genome re-sequencing data with high-depth from 5 RJFs and 8 village chicken (VCs) bred as a food source in the villages of Yunnan Province, Southwest China, an area believed to be one of the origins for domestic chickens8 (Supplementary information, Table S1). We also obtained RNA-seq data for the retina and four regions of the brain (cerebral cortex, corpus striatum, optic lobe, and cerebellar vermis) from RJF and VC. RNA-seq libraries with insert size around 300 bp were prepared using the Illumina standard RNA-seq library preparation pipeline and sequenced on the Illumina Hiseq 2000 platform (Supplementary information, Table S2). In total, transcriptomes of three retinas, two cerebral cortex, two cerebellar vermis, one corpus striatum and two optic lobes for VCs and two retinas, one cerebral cortex, one cerebellar vermis, one corpus striatum and one optic lobe for RJFs were used to perform comparative analysis. In addition, another ten transcriptomes were re-sequenced to verify the results, including six transcriptomes of three retinas from three RJF, and four transcriptomes of two retinas from two VCs with two experimental repeats per sample on the Hiseq 4000 platform (Supplementary information, Table S2). Sequencing data from this study have been submitted to the NCBI Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra) under accession number SRP040477.

Genomic sequence alignment and genotyping

Sequence reads were filtered by removing adaptors and low-quality bases using cutadaptor and Btrim61 (http://graphics.med.yale.edu/trim/). Only qualified paired-end reads were aligned onto the chicken reference genome (Galgal4) using BWA-MEM with default settings, except “-t 8 -M” options (https://github.com/lh3/bwa). A series of post-processes were then carried out with the alignment bam format file, including sorting, duplicate marking, local realignment and base quality recalibration. All of them were carried out using the relevant tools from the Picards 1.56 (http://picard.sourceforge.net) and Genome Analysis Toolkit 2.662 (GATK). SNPs and indels were called and filtered using UnifiedGenotyper and VariantFiltration command in GATK. Loci with more than 2 alleles were removed. All SNPs were assigned to specific genomic regions and genes using ANNOVAR63 based on the Ensembl chicken annotations.

A total of 17 115 375 SNPs were identified, around half (51.6%) of which were mapped to intergenic regions in the genome (Supplementary information, Table S3). Among the SNPs that locate to protein-coding regions, 70 990 were non-synonymous and 174 748 were synonymous, with 445 genes (including 446 transcripts) having SNPs that cause the gain or loss of a stop codon (Supplementary information, Table S4).

Genome-wide selective sweep scans

To identify genomic regions harboring footprints of positive selection in native chickens, we employed multiple tests to investigate selection in both populations and in a single population. Here, we used the FST, Pi, XP-EHH and XP-CLR statistical methods. FST values for each SNP were calculated between RJF and VC as previously described10. Nucleotide diversities (Δπ or ΔPi) = πRJFπVC were calculated using a sliding window analysis with a window size of 50 kb and a step size of 25 kb. XP-EHH value for each variant was calculated according to11. XP-CLR test12 was performed with scripts available at http://genetics.med.harvard.edu/reich/Reich_Lab/Software.html using the following parameters: sliding window size 0.1 cM, grid size 10 k, maximum number of SNPs within a window 300, correlation value for 2 SNPs weighted with a cutoff of 0.99. Here, VC was taken to be the object population, and RJF was chosen as the reference population. Haplotypes for each chromosome were deduced by shapeit.v2.r72764 (http://www.shapeit.fr/). General genetic map for chicken used here was 2.8 cM/Mb for chr1-9 and 6.4 cM/Mb for chr10-28 and chr3265.

SNP validation in larger population

We choose 67 SNPs that were mapped to the TSHR, ATP13A4, RPGRIP1L, COL18A1, VIT, OPA1, IGF2BP3 and MAP3K5 genes (see Supplementary information, Table S5) to validate in larger population, including 25 RJFs and 20 VCs by Sanger sequencing. Primer pairs used for these genes are listed in Supplementary information, Table S6.

Demographic history and coalescent simulations

Demographic histories for RJF and VC were inferred using MSMC66 based on haplotypes from multiple individuals in each population. We performed our analysis using 2 individuals (4 haplotypes), 3 individuals (6 haplotypes), 4 individuals (8 haplotypes) and 5 individuals (10 haplotypes) for each group (VC and RJF) independently. To assess the significance of the selective sweep signals identified by the above-mentioned methods, we stimulated two groups of genome sequences (one for VC and one for RJF) under a neutral evolutionary model taking into account the inferred demographic history. Coalescent simulations were performed using the ms program67. In total, 26 sequences of 50 kb were simulated 1 000 times and applied to the FST and Pi test. In addition, 26 sequences of 1 Mb were simulated 1 000 times and the mean XP-EHH value of the middle 50 kb (475-525 kb) in each simulated sequence was calculated as simulated data. Population sizes for the common ancestors (> 8 kya and 2 kya) of RJF and VC were obtained from our MSMC analysis. Since MSMC has a low power to estimate population size at relatively recent times, the effective population sizes for present day RJF and VC populations were taken from elsewhere68 (1.6 × 105 and 4 × 105, respectively). Generation time (g) and mutation rate per year (u) for chicken used here is 1 year and 1.91 × 10−9, respectively69.

Scripts used for simulation were as follows:

For FST and Pi:

ms 26 1000 -T -I 2 10 16 -t 191 -en 0.003 1 0.16 -en 0.003 2 0.22 -en 0.001 1 0.8 -eg 0.001 1 0 -en 0.001 2 0.32 -eg 0.001 2 0 -ej 0.004 1 2 | tail -n +4 | grep -v // >fst-pi-treefile

seq-gen -mHKY -l 50000 -s .017 <fst-pi-treefile> fst-pi-seqfile.ms

For XP-EHH:

ms 26 1000 -I 2 10 16 -t 3820 -en 0.003 1 0.16 -en 0.003 2 0.22 -en 0.001 1 0.8 -eg 0.001 1 0 -en 0.001 2 0.32 -eg 0.001 2 0 -ej 0.004 1 2 -p 6 >XP-EHH.ms

Test for the effect of relaxation of constraints and hitchhiking

Relaxation of purifying selection should increase pN/pS (the ratio of the number of non-synonymous substitutions per non-synonymous site (N) to the number of synonymous substitutions per synonymous site (S)) as selection against non-synonymous substitutions is lowered41,42. Here, we retrieved SNPs in vision-related genes identified under potential positive selection by the four methods, according to functional category enrichment, and calculated the ratio of the number of pN SNPs to the number of synonymous SNPs (N/S) in the VC and RJF populations. To exclude the possibility of hitchhiking, we compared the FST, Pi and XP-EHH values for the visual-related candidate genes with their related upstream/downstream sequences of the same length as the gene, sequences at distance of 50 kb and in adjacent genes.

Comparison of gene expression

The RNA-seq data were trimmed to filter out adaptor sequences and low-quality bases using cutadaptor and Btrim (http://graphics.med.yale.edu/trim/) and were mapped onto the chicken reference using Tophat43,70. Cufflinks and cuffcompare43 were then used to assemble transcripts and compare the assembled transcripts with the annotated reference transcripts, generating a newmerged GTF annotation file. New transcripts assembled here were defined as: having at least 2 exons and length ≥ 200 bp with class code j, c, o, e, and new loci with all class codes being u, or i, or x. Cufflinks was run using the GTF file with parameter “-G” to retrieve expression profile for each gene. Cuffdiff43 was used to detect the significance of the gene expression differences from four brain regions (cerebral cortex, corpus striatum, optic lobe, and cerebellar vermis) and retina between the RJF and VC. P-value was corrected by FDR (default by cuffdiff program).

To assess whether the expression changes in retina, cerebral cortex, corpus striatum, optic lobe, and cerebellar vermis between RJF and VC may have been driven by positive selection at local regulatory sites during domestication, a series of statistic tests were performed. The expression level (FPKM) for each gene in each tissue was retrieved and transformed according to log2(FPKM + 1)68. The difference of expression level for each gene between VC and RJF was calculated using log2((FPKM + 1)VC/(FPKM + 1)RJF). We then compared the expression levels of PSGs identified by FST, Pi, XP-EHH and XP-CLR with those of other genes (all genes in whole genome excluding PSGs) by Mann-Whitney U test and Kolmogorov-Smirnov test in each tissue. P-values were corrected by FDR. Since the two different batches of data employed different sequencing platforms and different runs, we performed the comparative analysis on data generated on the Hiseq 2000 and the Hiseq 4000 platforms independently, and also analyzed the combined data generated on both platforms.

Selective sweep region annotation and gene functional enrichment analysis

Candidates' selective sweeps detected by the above-mentioned methods were annotated using the Variant Effect Predictor available at http://asia.ensembl.org/info/docs/tools/index.html. Functional enrichments of the protein-coding genes including GO categories, KEGG pathway and HPO were analyzed using g: Profiler71.

VIT expression data in mouse

Expression data for VIT in different tissues of the mouse was retrieved from Biogps (http://biogps.org/). Expression profiling by array in retinal pigment epithelium (RPE) of mouse exposed to 10 000 lux of cool white fluorescent light for 18 h was from a previous study72 (GSE37773). The expression data from retina after Hmx1 loss and in wild-type mice were from another study73 (GSE47002). Expression data from photoreceptor cells of Nrl−/− mice and Crx::Nr2e3/Nrl−/− mice were also from a previous study46 (GSE5338).

RNA extraction and real-time quantitative PCR assay

Total RNA was isolated from chicken retina tissues using TRNzol-A+ (Tiangen Biotech, Beijing, China) and purified using RNeasy Micro kit (Qiagen, Germany). Concentration and integrity of the RNA was measured using electrophoresis and NanoDrop spectrophotometer 2000. Total RNA (2 μg) was used to synthesize single-strand cDNA using the PrimeScript RT-PCR Kit (TaKaRa, Japan) in a final volume of 25 μL reaction mixture according to the manufacturer's instructions. The relative mRNA expression levels of RPGR (retinitis pigmentosa GTPase regulator), GUCAIA (guanylate cyclase activator 1A), TRPM1 (transient receptor potential cation channel, subfamily M, member 1), PDE6B (phosphodiesterase 6B), NR2E3 (nuclear receptor subfamily 2, group E, member 3), VIT (vitrin) and RHO (rhodopsin) genes were measured using real-time quantitative PCR (qPCR) and the relative standard curve method, with normalization to the house-keeping gene GAPDH. qPCR was performed on the platform of the iQ2 system (BioRad Laboratories, Hercules, CA, USA) with SYBR Premix Ex Taq II kit (TaKaRa, DRR081A). Student's t-test was used to analyze the differences and measure the statistical significance. Primer pairs used for these genes are listed in Supplementary information, Table S7.

Zebrafish care and maintenance

Adult wild-type AB strain zebrafish were maintained at 28.5 °C on a 14 h light/10 h dark cycle74. Five to six pairs of zebrafish were set up for each natural mating each time. On average, 200-300 embryos were generated. Embryos were maintained at 28.5 °C in fish water (0.2% Instant Ocean Salt in deionized water). The embryos were washed and staged according to published guidelines75. The zebrafish facility at Shanghai Research Center for Model Organisms is accredited by the Association for Assessment and Accreditation of Laboratory Animal Care (AAALAC) International.

Zebrafish microinjections

Gene Tools, LLC (http://www.gene-tools.com/) designed the MO oligonucleotide. Antisense MO (GeneTools) was microinjected into fertilized one-cell stage embryos according to standard protocols76. The sequence of the exon 3-intron 3 splice VIT MO (VIT-e3i3-MO) was 5′-CCTGAATAGTCTACAGTACCTGAGA-3′. For the VIT knockdown experiment, 4 ng of VIT-e3i3-MO was used per injection. Total RNA was extracted from 80 to 100 embryos per group in Trizol (Invitrogen) according to the manufacturer's instructions. RNA was reverse transcribed using the PrimeScript RT reagent Kit with gDNA Eraser (Takara). Primers spanning VIT exon 2 (forward primer: 5′-GTTCGCCTCCATATCCAGCATCTGC-3′) and exon 9 (reverse primer: 5′-TTCCCCTGGACGCTCTGAGACTGTC-3′) were used for RT-PCR analysis for confirmation of the efficacy of the E3I3-MO.

AO staining for apoptosis

Wild-type control embryos and embryos injected with 4 ng VIT-e3i3-MO were immersed in 5 μg/ml AO (acridinium chloride hemi-[zinc chloride], Sigma-Aldrich) in fish water for 60 min at 32-hpf and 80-hpf. Next, zebrafish were rinsed thoroughly in fish water three times (5 min/wash) and anaesthetized with 0.016% MS-222 (Tricaine methanesulfonate, Sigma-Aldrich, St. Louis, MO). Zebrafish were then oriented on their lateral side and mounted with methylcellulose in a depression slide for observation by fluorescence microscopy.

In vivo ERG recording

Zebrafish larvae were first paralyzed with α-bungarotoxin (100 μg/ml, Sigma), and then embedded in 1% low melting-point agarose (Sigma) with one eye upward in a custom-made chamber as described previously77. The cornea and lens of the upward eye were removed to expose retinal surface by using a glass micropipette with a tip opening of 1 μm. After dissection, the larvae were transferred to a recording setup, and perfused with extracellular solution, which consists of (in mM) 134 NaCl, 2.9 KCl, 4 CaCl2, 10 HEPES and 10 glucose (290 mOsmol/l, pH 7.8). ERG responses were recorded with an EPC-10 amplifier (Heka, Germany) through inserting the pipette of 3-μm tip opening into the interface of photoreceptors and bipolar cells. The ERG signals were amplified at 1 000 total gain and low-pass filtered at 100 Hz. Whole-field light stimuli were given by a white LED controlled by the Master 8 stimulator (A.M.P.I., Israel).

Visual escape behavioral test

Larval behavior was monitored by an infrared-sensitive high-speed camera at the acquisition rate of 250 Hz (Redlake Motionscope M3, US). In one experiment, the behaviors of 12 larvae were simultaneously recorded. Each larva was put in an individual 3.5-cm Petri dish under the light background for 30 min to adapt. During each test, 10 trials were carried out for calculating the probability of escape behavior with a 5-min interval. A successful escape was scored when a C-shape movement finished within 20 ms after movement onset78. Visual stimuli were generated by a self-written Matlab program and delivered through a mini-projector. The dark shadow would gradually cover the fish's body when an expanding black disc was projected on the Petri dish at the expanding speed of 5.4 cm/s.

Image acquisition

Embryos and larvae were analyzed with Nikon SMZ 1500 fluorescence microscope and subsequently photographed with digital cameras. A subset of images was adjusted for levels, brightness, contrast, hue and saturation with Adobe Photoshop 7.0 software (Adobe, San Jose, CA, USA) to optimally visualize the expression patterns. Quantitative image analyses processed using image-based morphometric analysis (NIS-Elements D3.1, Japan) and ImageJ software (National Institutes of Health, Bethesda, MD, USA; http://rsbweb.nih.gov/ij/). Inverted fluorescent images were used for processing. Positive signals were defined by particle number using ImageJ. In total, 10 animals for each treatment were analyzed and the total signals per animal were averaged. All data were presented as mean ± SEM. Statistical analysis and graphical representation of the data were performed using GraphPad Prism 5.0 (GraphPad Software, San Diego, CA, USA). Statistical significance was performed using a Student's t-test as appropriate. Statistical significance is indicated by *, where P < 0.05, and ***, where P < 0.0001.

Author Contributions

YPZ and DDW led the project and designed the study. MSW, RWZ, LYS and JLD performed experiments. MSW, DDW, LYS, YL, MSP, RWZ, JLD, HQL and LZ performed the analyses. MSW, HQL and LZ sampled chickens and tissues. MSW, DDW, DMI and YGY drafted the manuscript. All authors read and improved the manuscript.

Competing Financial Interests

The authors declare no competing financial interests.