Two Proximally Close Priority Candidate Genes for diplopodia-1, an Autosomal Inherited Craniofacial-Limb Syndrome in the Chicken: MRE11 and GPR83.

Next-generation sequencing (NGS) and expression technologies were utilized to investigate the genes and sequence elements in a 586 kb region of chicken chromosome 1 associated with the autosomal recessive diplopodia-1 (dp-1) mutation. This mutation shows a syndromic phenotype similar to known human developmental abnormalities (e.g., cleft palate, polydactyly, omphalocele [exposed viscera]). Toward our goal to ascertain the variant responsible, the entire 586 kb region was sequenced following utilization of a specifically designed capture array and to confirm/validate fine-mapping results. Bioinformatic analyses identified a total of 6142 sequence variants, which included SNPs, indels, and gaps. Of these, 778 SNPs, 146 micro-indels, and 581 gaps were unique to the UCD-Dp-1.003 inbred congenic line; those found within exons and splice sites were studied for contribution to the mutant phenotype. Upon further validation with additional mutant samples, a smaller subset (of variants [51]) remains linked to the mutation. Additionally, utilization of specific samples in the NGS technology was advantageous in that fine-mapping methodologies eliminated an additional 326 kb of sequence information on chromosome 1. Predicted and confirmed protein-coding genes within the smaller 260 kb region were assessed for their developmental expression patterns over several stages of early embryogenesis in regions/tissues of interest (e.g., digits, craniofacial region). Based on these results and known function in other vertebrates, 2 genes within 5 kb of each other, MRE11 and GPR83, are proposed as high-priority candidates for the dp-1 mutation. © The American Genetic Association 2018.

Introduction

Animal models contribute enormously to our understanding of the genetic and molecular basis of human diseases and disorders. The chicken embryo offers a unique advantage to the study of vertebrate developmental biology, due to the many naturally occurring mutants that exhibit developmental defects (Pisenti et al. 1999; Delany 2004; Robb et al. 2011; Robb and Delany 2012a) as well as the benefit of in ovo embryogenesis allowing for easy access, assessment, and manipulation. The many benefits of utilizing the chick embryo as a model system for human health-related research has been described in both primary research articles and reviews (Schock et al. 2016; Davey et al. 2018; Stern 2018; The Chick Embryo Model System 2018). One area of biomedical research that makes excellent use of the chick embryo model system is the study of inherited congenital malformations (Burt 2007; Davey and Tickle 2007; Schock et al. 2016; Davey et al. 2018). The inherited mutations in the chicken contribute new knowledge regarding the etiology and pathways underpinning complex syndromes involving craniofacial, limb and organ malformations (Davey et al. 2006; Schock et al. 2016). Diplopodia is a congenital malformation occurring in tetrapods involving duplication of elements of the foot and limb as well as other structural abnormalities. Human and chicken diplopodia (Figure 1) share phenotypic features (Karchinov 1973; Jones et al. 1978; Narang et al. 1982; Hamanishi et al. 1985; Brower et al. 2003; Khan et al. 2008; Kadir et al. 2011) including extreme polydactyly (more than one additional digit), micromelia (truncated limbs), and visceral abnormalities. Although similar to polydactyly in that supernumerary digits are duplicated, diplopodia differs from common polydactyly in that affected individuals have additional metatarsal and tarsal bones in addition to the extra digits. Diplopodia, derived from the Greek roots diplo, meaning “double” and pod meaning “foot,” refers to the pre-axial polydactyly observed in the mutant embryo (Figure 1), resulting in the doubling of the structures of the foot (i.e., metatarsals) (Taylor and Gunns 1947).

Figure 1.

Developmental phenotype variation of diplopodia-1 mutant embryos. All individuals shown are at 10 days of embryogenesis. A. Normal (+/+) UCD-Dp-1.003 chicken embryo. B. UCD-Dp-1.003 mutant embryo (−/−) displaying truncation of the limbs, mild cleft-palate, mild dwarfism (shorter stature), exposed viscera (omphalocele), and pre-axial polydactyly (5 digits on both wings and legs). C, D. Wing digit number variation is observed in UCD-Dp-1.003 mutants (5 digits shown in C; 4 digits shown in D; normal digit number is 3). E, F. Hind limb (leg) digit number variation is observed in UCD-Dp-1.003 mutants. Note the inter-individual variation of digit length and individual morphology (i.e., webbing versus no webbing at same stage of development) in the wings and legs. G–J. Variation of visceral exposure is observed in UCD-Dp-1.003 mutant embryos. K–N. Variability in clefting (mild through severe) is also observed in UCD-Dp-1.003. Note that dp-1 mutant embryos displaying severe clefting, a phenotype not as common but yet observed over the course of this study, were not available on the day of photography and it was deemed inappropriate/unnecessary to set additional eggs for observations when the severe phenotype representation is available from other lines (e.g., coloboma—see Robb et al. 2013) at same age (i.e., in terms of reducing the usage of animals).

Diplopodia-1 (dp-1) is 1 of 5 diplopodia mutations (named dp-1 through dp-5) described in chicken and was shown to be non-complementary with the other 4 mutations (Robb et al. 2011 and references therein). The morphological severity of the wing and leg defects was assessed and ordered as such: dp-3 < dp-1 < dp-4 < dp-2 (Taylor 1972). The dp-1 mutation originated within a Single Comb White Leghorn commercial stock, was acquired by researchers at the University of California in 1941, and later incorporated into a congenic inbred line (a.k.a. UCD-Dp-1.003) in the late 1980s to early 1990s. This mutation is characterized by dwarfism, craniofacial defects (e.g., cleft palate), exposed visceral organs (omphalocele), micromelia, and polydactyly ranging from the addition of 1–5 pre-axial digits (pre-appendage) affecting the legs and wings (Taylor and Gunns 1947; Landauer 1956; Robb et al. 2011; Robb 2012). Digit pattern conformation analyses in dp-1 embryos showed that a majority (~80%) of mutants had 6 digits on each foot arranged in 2 sets of 3, in a 2ʹ-1ʹ-1-2-3-4 digit conformation (Landauer 1956). (The primed numbers [e.g., 1ʹ] refer to the extra, duplicated digit; 1-2-3-4 is the foot pattern of a normal embryo.) During the characterization studies conducted by Taylor and Gunns (1947), dp-1 was found to be an autosomal recessive, embryonic lethal mutation affecting both sexes. Over half a century after this mutation was initially studied for mode of inheritance, Robb et al (2011) mapped the dp-1 mutation to a 708-kb region on GGA 1 (p arm) using an SNP array (formerly 720 kb based on the galGal3 assembly).

Herein we describe our research conducted as part of a long-term strategy to identify the causative element responsible for the chicken dp-1 phenotype, which involved breeding, genetic, genomic and developmental tools of classical and modern biology. Fine-mapping and identification of recombinant individuals narrowed the region to 586 kb. To advance the aim to discover the causative element, a targeted genomic capture enrichment (a.k.a. capture array [CA]) was paired with next-generation sequencing (NGS) to sequence, in its entirety, the region linked to the UCD-Dp-1.003 mutation. Bioinformatics of the NGS results identified sequence variants (e.g., SNPs, micro-indels, gaps) specific to the dp-1 mutation, while multiple pairwise-line comparative genomic analyses and validation tests identified those unique (UCD-Dp-1.003-specific) elements, which could be responsible for the mutant phenotype. Lastly, RNA in situ hybridization was utilized to assess the spatial and temporal expression patterns of known genes in the linked region.

Methods

Chicken Genetic Lines

The individuals utilized for this study were from 2 genetic lines, the developmental mutant-congenic inbred line UCD-Dp-1.003 (aka Dp-1.003) and the inbred (F > 0.99) parent background line UCD-003 (Abplanalp 1992; Robb et al. 2011). Animals were under the care and supervision of trained staff and as per an approved protocol by the UC Davis Institutional Animal Care and Use Committee (Protocol # 15439). Associated with the studies described below, an in-depth assessment of mutant (n = 56) phenotypic variation, in the context of the inbred congenic line background, was conducted.

Sample Collection for Fine-Mapping

Dp-1.003 embryos were incubated to E10, a stage of development such that the phenotypes (normal: +/+ or +/−; mutant: −/−) could be easily and accurately discerned (Figure 1). Adult or embryonic blood samples were collected according to Robb et al. (2011) and pin feathers were collected from day-of-hatch chicks (Robb 2012). DNA (with an RNAse step) was isolated from blood and tissue (feather pulp, brain) sources using the DNeasy® Blood & Tissue kit (Qiagen) and examined by electrophoresis to determine integrity. See Robb and Delany (2012b) for preparation of samples used in the capture array (CA) technology.

SNP-Genotyping, Fine-Mapping Analysis, and CR Identification

Eight SNPs (rs14931758–rs13990802; Table 1), in complete linkage disequilibrium with the dp-1 trait (Robb et al. 2011), were used to fine-map the region. DNAs (10 ng/µL) isolated from Dp-1.003 homozygous normal (+/+, n = 54), heterozygous normal (+/−, n = 106), and mutant (−/−, n = 78) samples were used to determine and/or confirm genotypes at the 8 loci using standard polymerase chain reaction (PCR) conditions. Reactions were amplified using Phire® Hot Start II DNA Polymerase (Thermo Scientific) and purified by QIAquick® Spin Kit (Qiagen) under manufacturer conditions. Amplicons were sequenced (Davis Sequencing, Davis, CA) using ABI 3730 DNA sequencers (Applied Biosystems) and were analyzed for genotype-specific SNP differences. The SNP analysis and causative region (CR) identification (maximum and minimum CR: CRmax and CRmin, respectively) were defined as previously described (Robb et al. 2011). A subset of the mutant samples (n = 20) were utilized in the variant validation portion of this study described below.

Table 1.

Diplopodia-1 fine-mapping primers used to identify carrier status, CR size, and recombination events

SNP	Positiona	Primers (5ʹ–3ʹ)
		Forward	Reverse	Product size (bp)b
rs13989579	185432685	GAGATGTGGAGGCCAAAAGT	CTAATGCCAACCTGGCTTCT	282
GGaluGA060936	185482293	CACAAAAATGGTGAATTCTTAACAT	GTTAGAGTTCACAAAGTAGAAGTCTGG	250
rs14931758	185587553	CAGTGACCAGCAGAAGGACA	TGGCCTGGTTAATGTGTGAA	192
rs14931949	185671027	GCCAAGTATCAACTCCCATC	CCACTCTGAAAAGGAGCATT	214
rs13989873	185695823	TTCTGCTATTCTGTGGCTCA	TGGAAATCAGAGTTGAAGCA	224
rs14932143	185805128	GTTACTTTCAGCAAGCACAG	GCACCATTTCAACACATCTA	250
rs13989974	185830327	CTGATGTGCTCAGGTAACAC	AGGAGACACTGATTTCCTTG	194
rs13990135	185909526	TCAAGAACTTGCTGGCATAA	GGTTCCTTAGTCTCCCCAGT	192
rs13990785	186281274	CCACAAATGCATTACCTGAG	TCCCCTACAGACCCTATTGA	169
rs13990802	186295291	ATCAACGCAGGGTCTCAACT	GGCTGTACAAATCTGATCATCATTAC	260

Chromosomal location (bp) of SNP on GGA 1; positions are based on the December 2015 Gallus gallus assembly (galGal5).

PCR fragment size was determined by 3 methods: 1) using the UCSC genome browser (http://genome.ucsc.edu/), 2) sizing by gel electrophoresis, and 3) DNA sequencing.

Capture Array and SOLiD™ Sequencing

Genomic Enrichment Services of SeqWright, Inc., were utilized to design a custom NimbleGen capture array to enrich for the 585 451-bp chromosomal region (GGA 1: 185,695,823-186,281,274 [galGal5]; formerly 189,915,923-190,511,744 [galGal3]) linked to the UCD-Dp-1.003 developmental mutation (Robb et al. 2011; Robb and Delany 2012b). Our strategy for the capture array development was to design probes for the 586 kb CRmax identified post-60K SNP array and subsequent fine-mapping. Two dp-1 mutant samples (dp1-166F and dp1-168F) were utilized in the CA/NGS after genotypic confirmation at the 8 SNP loci listed in Table 1. Both samples exhibited a reduced CRmax, indicated by SNP analyses. Specifically, dp1-166F was utilized because the mutant reduced the region at the 5′ end, whereas dp1-168F narrowed the region at the 3′ end of the linked region (Table 2). Full details of the array setup and methods can be found in Robb and Delany (2012b). (Note that the same array was utilized to sequence 2 other UCD-congenic developmental mutations, coloboma [co.003] and wingless-2 [wg-2.331]; Robb and Delany 2012b; Robb et al. 2013; Webb et al. 2018.) The other 2 congenic lines, which mapped to different chromosomes (GGA Z and 12), were utilized as controls for GGA 1 comparisons (Robb et al. 2011; Robb 2012; Robb and Delany 2012b).]

Table 2.

SNP fine-mapping shows a reduced CR via assessment of linked molecular markers in the UCD-Dp-1.003 genetic line

Original 60K SNP array mutant samplesa

Post-60K SNP array mutantsb

Post-60K SNP array heterozygotesc

Overall reduced genotyped

Normal (+/+)e

Sample ID No:f

n = 2

n = 3

n = 5

n = 7g

n = 3h

n = 2

n = 2

n = 2

n = 3

n = 6

SNP ID

SNP

Positioni

Genotype:j

rs13556505

A/G

177307260

AG

GG

AA

AG

GG

AG

GG

GG

GG

GG

AA

AA

rs13989579

A/G

185432685

GG

GG

AA

AG

GG

AG

GG

GG

GG

GG

AA

AA

GGaluGA060936

A/G

185482293

AA

AA

GG

AG

AA

AG

AA

AA

AA

AA

GG

GG

rs14931758

A/G

185587553

GG

GG

AA

AG

GG

AG

GG

GG

GG

GG

AA

AA

rs14931949

A/C

185671027

CC

CC

CC

AC

CC

AC

CC

CC

CC

CC

AA

AA

rs13989873

A/G

185695823

GG

GG

GG

AG

GG

AG

AG

AG

GG

GG

AA

AA

rs14932143

A/C

185805128

CC

CC

CC

CC

CC

CC

AC

AC

AC

AC

CC

AA

rs13989974

A/G

185830327

AA

AA

AA

AA

AA

AA

AG

AG

AG

AG

AA

GG

rs13990135

A/G

185909526

AA

AA

AA

AA

AA

AA

AG

AG

AG

AG

AA

GG

rs13990785

A/G

186281274

GG

GG

GG

GG

AG

AG

AG

AA

AG

GG

AA

AA

rs13990802

A/G

186295291

GG

GG

GG

GG

AG

AG

AG

GG

AG

GG

GG

GG

rs13990804

A/G

186316405

AA

AA

AA

AA

AG

AG

AG

AA

AG

AA

AA

AA

60K CRmin:k

610,247

Individual CRmin:l

848,589

848,589

610,247

476,146

8,602,266

104,398

620,582

213,703

511,277

104,398

104,398

60K CRmax:k

707,738

Individual CRmax:l

8,988,029

9,447,557M

707,736

599,466

9,474,014

585,449

1,145,378

610,245

1,120,582

585,451

585,449

Samples were assessed for their 60K SNP genotyping pattern to identify recombination events, thereby reducing the size of the linked region. This CRmax decreased by 122-kb post-60K SNP array analysis. Reducing the size of the CR is beneficial as it: 1) eliminates candidate genes in question, 2) eliminates genomic regions (as other elements could be the cause of the mutation), and 3) allowed for the identification of samples with a narrowed CR to be used in the capture array technology.

Samples originally genotyped using the Illumina 60K SNP genotyping array (Robb et al. 2011).

A total of 78 mutant embryos were collected post-60K SNP analysis (yr 2009–2011). These were genotyped using the 60K SNPs (Table 1), but only 12 displayed recombination events. Two of these samples, dp1-166F and dp1-168F, were utilized in the targeted genomic capture enrichment technology.

A total of 73 heterozygous chicks were hatched post-60K SNP analysis (yr 2009–2011). All were genotyped using the 60K SNPs (Table 1) and 13 displayed recombination events.

The overall reduction, combining the fine-mapping results of both mutant and heterozygous samples, is displayed. Three SNPs remain linked to the Dp-1.003 mutation, thereby maintaining a CRmin and CRmax of 104,398 and 585,449 bp, respectively.

SNP genotype observed in normal, control samples.

Please note that all individuals used in reproduction were Dp-1.003 carriers (+/−) except for female 18539F (used in 15 matings) and male 10775M (used in 1 mating) which were UCD-003 (congenic background).

Note: targeted genomic capture enrichment technology sample, dp1-166F, showed this 60K SNP genotype pattern.

Note: targeted genomic capture enrichment technology sample, dp1-168F, showed this 60K SNP genotype pattern.

Chromosomal location (bp) of SNP on GGA 1; positions are based on the December 2015 Gallus gallus assembly (galGal5). SNPs were identified using the 60K SNP array (Robb et al. 2011).

Genotypes shaded in gray indicate a region which is no longer linked to the Dp-1.003 mutation.

UCD-Dp-1.003 minimum and maximum CR (candidate gene region; CRmin and CRmax) identified by Robb et al. (2011) using the chicken 60K SNP genotyping array.

Minimum and maximum CR (candidate gene region) identified for each sample.

The recombination rate at this position of GGA 1 is 2 cM/Mb (Groenen et al. 2009). In order to best estimate the CRmax, 500 kb was subtracted from the 5′ sequence coordinate of the CRmin.

Unique Variant Identification and Causative Element Analysis

The CA/NGS data were analyzed to identify SNPs, micro-indels (1–3 nt), gaps, and chromosomal rearrangements. Reads were aligned to the Gallus gallus NCBI UCD-001 Red Jungle Fowl (RJF) reference genome sequence (WASHUC2, May 2006). Reference-assisted assembly and variant identification procedures are described fully in Robb and Delany (2012b); Dp-1.003-specific variant analysis is briefly reviewed below.

The discrimination of non-causative polymorphic elements (naturally occurring in the introgressed region) from the potentially causative polymorphic elements was accomplished by a comparative approach. Polymorphisms were identified by comparing the Dp-1.003 causative (sequenced) region to variants identified in other chicken genetic lines including: the other inbred congenic lines (wg-2, Webb et al. 2018, and co, Robb et al. 2013) utilized on the capture array, the chicken reference genome (UCD-001), or any other previously reported polymorphism (e.g., within NCBI [ex. dbSNP] and the UCSC genome browser). Any polymorphism (e.g., SNP, indel, etc.) found in the Dp-1.003 GGA 1 CR, not found within any of the other sources, was deemed dp-1-specific (i.e., unique) and therein considered as a possible causative element for the dp-1 mutation. Newly identified, shared variants were submitted to NCBI (accessions: ss472337944–ss472340673).

To further eliminate normal variation, not causing the dp-1 phenotype, unique sequence elements were assessed in a set of new mutant samples, the progeny from different matings. DNAs isolated from 20 mutant (−/−) and 2 control (UCD-003, +/+, F > 0.999) individuals were used to assess the dp-1-specific elements identified by pairwise comparisons described above within the newly reduced 260 kb region on GGA 1 identified by fine-mapping of the 2 mutant samples utilized in the CA.

Specifically, all polymorphisms located within an exon, UTR, or splice site within the original 586 kb were assessed for linkage. Primers were designed to flank the variant using sequences available through the UCSC genome browser (http://genome.ucsc.edu/) (primer sequences are available upon request). Elements were PCR amplified (standard conditions) and amplicons were sequenced. Each element was then evaluated using polymorphic marker pattern analyses described in Robb et al. (2011). Briefly, a variant can be eliminated as causative if it is present in control (+/+) samples or if the polymorphism is absent in additional mutant samples (beyond the 2 individuals sequenced as part of the CA/NGS study). Using this strategy, unique SNPs, micro-indels, and gaps found within an exon or splice site (thus potentially impacting mRNA and protein) of known and predicted genes (Table 3)were tested in an additional 20 mutant and 2 control (UCD-003) samples to verify if a particular variant remained linked to the dp-1 mutation (and thus could be potentially causative). Table 5 outlines the number of each element found at each stage of the bioinformatics analysis.

Table 3.

Thirteen genes found within the Dp-1.003 260 475 nt CR identified by fine-mapping techniques

Candidate gene	Synteny chromosomea		Gene functiond	Knockout mouse phenotypee	Affected anatomical systemse	Associated diseasesf
	Mouseb	Humanc
MRE11 g	Mmu 9	Hsa 11	Component of MRN (Mre11a, Rad50, Nbs1) complex, which is involved in DNA repair in animals. Specifically, it has a role in DNA recombination, double-strand break repair, and maintenance of telomere length and meiosis. Disruption results in defective sister chromatid cohesion.	Mortality/aging, embryogenesis, behavior, tumorigenesis; lethal: embryos developmentally arrest at e9.5 with a poorly defined body axis, no heart development	Reproductive, cardiovascular	Human: Ataxia telangiectasia-like disorder (ATLD); breast, lymphoid, and colorectal cancers; Alzheimer’s disease; Cornelia de Lange Syndrome I
						Mouse: Associated with several identified QTLs (e.g., tuberculosis severity 2, and others)
GPR83	Mmu 9	Hsa 11	G-protein coupled receptor; endogenous GPR83 plays a role in the induction of Foxp3 expression in peripheral T cells. Gpr83 may be involved in the central control of the energy metabolism.	Immune, hematopoietic, lean, resistant to diet-induced obesity and resistant to diet-induced glucose intolerance	Normal regulatory T cell development	Human: Associated with several identified QTLs (e.g., blood pressure, body weight, prostate tumor susceptibility, and others)
						Mouse: Associated with several identified QTLs (e.g., cholesterol 6, Salmonella enteritidis susceptibility 6, and others)
PANX1	Mmu 9	Hsa 11	Plasma membrane protein which mediates “find-me” signal release from apoptotic cells. Associated with the initiation and propagation of calcium waves; may regulate ER Ca2+ homeostasis.	Nervous, immune, and respiratory systems, behavior	Impaired macrophage recruitment, ATP release by apoptotic astrocytes, hippocampal neurons, and thymocytes	Human: Associated with several identified QTLs (e.g., insulin level, chronic obstructive pulmonary disease, and others)
						Mouse: Tumor-suppressive effects
HEPHL1	Mmu 9	Hsa 11	Hephaestin-like protein 1. Function inferred by similarity: may be involved in copper transport and homeostasis as well as may function as ferroxidase	In progress	In progress	Human: Associated with several identified QTLs (e.g., prostate tumor susceptibility, serum adiponectin level, and others)
						Mouse: Associated with several identified QTLs (e.g., haloperidol-induced catelepsy 5, cholesterol 6, and others)
CP	Mmu 3	Hsa 3	Ceruloplasmin (ferroxidase enzyme) helps move iron from organs/tissues into blood stream by preparing iron for incorporation into transferrin (which transports iron to RBCs); also binds copper in plasma.	Behavior, immune system, homeostasis, hematopoietic	Vision, cardiovascular, nervous system, liver, pigmentation	Human: Aceruloplasminemia; Wilson’s disease; Obsessive compulsive disorder; Schizophrenia
						Mouse: Associated with several identified QTLs (e.g., organ weight, HIV-associated nephropathy1, and others)
HEPH	Mmu X	Hsa X	A multi-copper oxidase critical for intestinal iron absorption as it regulates transport of dietary iron from epithelial cells of the intestinal lumen into the circulatory system.	Embryogenesis, cardiovascular, hematopoietic, digestive/alimentary, nervous system	Impaired iron transport to gut and placenta, small and pale at birth, exhibit hypochromic anemia (disappears with age)	Human: Colorectal cancers
						Mouse: Severe microcytic, hypochromic anemia; age-related macular degeneration
TM4SF1a	Mmu 3		A cell surface antigen belonging to the tetraspanin family, which is characterized by the presence of 4 hydrophobic domains. Members of this transmembrane family are cell-surface proteins, which typically mediate signal transduction events involving the regulation endothelial cell movement and proliferation and growth; commonly associated with cancer cell proliferation.	Embryonic lethal at day 9 due to vasculogenesis failure	Cardiovascular	Human: Breast, colon, lung, and ovarian carcinomas
						Mouse: NA
VSTM5	Mmu 9	Hsa 11	Cell adhesion membrane glycoprotein, which regulates neuronal morphology and synapse formation during cortical development in the brain	—	—	Human: Unknown
						Mouse: Unknown
C11orf54	Mmu 9	Hsa 11	Exhibits ester hydrolase activity on the substrate p-nitrophenyl actetate	In progress	In progress	Human: Associated with several identified QTLs (e.g., blood pressure, body weight, insulin level, and others)
						Mouse: Associated with several identified QTLs (e.g., cholesterol 6, tuberculosis severity 2, and others)
MED17	Mmu 9	Hsa 11	Serves as a scaffold for the assembly of a functional preinitiation complex with RNA polymerase II and the general transcription factors; without MED17 and TFIID, SP1 cannot be activated.	In progress	In progress	Human: Microcephaly postnatal progressive with seizures and brain atrophy (MCPHSBA)
						Mouse: Associated with several identified QTLs (e.g., haloperidol-induced catelepsy 5, cholesterol 6, and others)
SNORD6	Mmu 9	Hsa 11	Small nucleolar RNA, C/D box 6 involved in a variety of important processes such as RNA splicing, regulation of transcription factors or RNA Polymerase II, and maintaining telomeres. Exact role unknown.	—	—	Human: Associated with several identified QTLs (e.g., body weight, insulin level, prostate tumor susceptibility, and others)
						Mouse: Unknown
SNORA8	Mmu 9	Hsa 11	Small nucleolar RNA, H/ACA box 8 involved in a variety of important processes such as RNA splicing, regulation of transcription factors or RNA Polymerase II, and maintaining telomeres. Exact role unknown.	—	—	Human: Associated with several identified QTLs (e.g., chronic obstructive pulmonary disease, body weight, and others)
						Mouse: Unknown
SMCO4	Mmu 9	Hsa 11	Single-pass membrane and coiled-coil domain-containing protein 4; hypothetically found in membrane and interacts with carbon nanotube, tributyltin oxide, troglitazone	—	—	Human: Unknown
						Mouse: Unknown

UCSC genome browser (http://genome.ucsc.edu/) was utilized to identify the location of each gene and the percent to the identity to human protein.

Synteny alignment was made to the mouse reference genome (December 2011; GRCm38/mm10) in the UCSC genome browser.

Synteny alignment was made to the human reference genome (December 2013; GRCH38/hg38) in the UCSC genome browser.

References for the gene function are as follows: MRE11: Petrini et al. (1995), Trujillo et al. 1998, Gatei et al. 2000, Barber et al. 2008; GPR83: Lu et al. (2007), Müller et al. (2013); PANX1 (aka Pannexin 1): Vanden Abeele et al. (2006), Chekeni et al. (2010), Kienitz et al. (2011); HEPHL1: Chen et al. (2010); CP: Holmberg and Laurell (1948), Hellman and Gitlin (2002); HEPH: Frazer et al. (2001), Vulpe et al. (1999), Prohaska (2011), C11orf54 (aka mouse: 4931406C07Rik): Manjasetty et al. (2006); MED17: Baek et al. (2006); SNORD6: Dieci et al. (2009); SNORA8: Dieci et al. (2009). SMCO4: https://rgd.mcw.edu/rgdweb/report/gene/main.html?id=1319742; TM4sf1a: Zukauskas et al. (2011). VSTM5: Lee et al. (2016).

Information for the knockout mouse phenotype and affected anatomical systems were obtained from the Mouse Genome Informatics (MGI, http://www.informatics.jax.org/), the International Knockout Mouse Consortium (IKMC, http://www.knockoutmouse.org/), and the Rat Genome Database (RGD) PhysGen Knockouts (http://rgd.mcw.edu/). A dash indicates that the particular gene was not present in the database.

References for the associated disease are as follows: MRE11: human: Stewart et al. (1999), Fukuda et al. (2001), Jacobsen et al. (2004), Sjöblom et al. (2006), mouse: http://rgd.mcw.edu/ ID=735478; GPR83: human: http://rgd.mcw.edu/ ID=1352259, mouse: http://rgd.mcw.edu/ ID=733083; PANX1: human: http://rgd.mcw.edu/ ID=1347653, mouse: Lai et al. (2007); HEPHL1: human: http://rgd.mcw.edu/ ID=1345614, mouse: http://rgd.mcw.edu/ ID=1622024; CP: human: Scheinberg and Gitlin (1952), Gitlin (1998), Wolf et al. (2006), Virit et al. (2008), mouse: http://rgd.mcw.edu/ ID=10384; HEPH: human: Brookes et al. (2006), mouse: Hahn et al. (2004), Griffiths et al. (2005); C11orf54: human: http://rgd.mcw.edu/ ID=1603055, mouse: http://rgd.mcw.edu/ ID=1319924; MED17: human: Kaufmann et al. (2010), mouse: http://rgd.mcw.edu/ ID=1322409; SNORD6: human: http://rgd.mcw.edu/ ID=1602744, mouse: NA; and SNORA8: human: http://rgd.mcw.edu/ ID=1602779, mouse: NA; TM4SF1: human: Marken et al. (1992).

Note that the 3′ end of MRE11 (11,462 nt inclusive of 14 conserved—human, mouse, frog, opossum–exons) remains linked to the Dp-1.003 mutation.

Table 5.

Diplopodia-1.003 SNPs, micro-indels, and sequence gaps: number and genomic location

CR sizea	Originally identified in CAb				Unique variantsc
	SNPs	Insertionsd	Deletionsd	Gapse	SNPs	Insertionsd	Deletionsd	Gapsef
585,451 nt	2593	150	133	581	778	82	64	581
260,477 nt	1509	73	63	296	696	47	39	296

CR size = 585,451 nt (based on the galGal5 coordinate alignment) is the CR identified by the 60K SNP array, which was subsequently utilized for capture array probe creation. 260,477 nt is the CR identified by fine-mapping analysis of recombinant individuals and CA/NGS bioinformatic analyses.

The original number of variants identified after bioinformatics analyses, prior to multiple pairwise-line comparisons to identify unique variants.

Unique variants are those specific to Dp-1.003 only after multiple pairwise genomic comparisons. See Methods section for details as to unique variant identification.

Insertions and deletions range from 1 to 3 nt in length (within the paper referred to as micro-indels).

Sequence gaps (DNA that was not captured for sequencing in the CA/NGS) are those gaps greater than 4 nt in length, as identified by alignment to the 585,451 nt RJF reference genome obtained from NCBI (see Methods section).

Gaps listed are those found within the region.

RNA In situ Hybridization

Eleven candidate protein-encoding genes (MRE11, GPR83, PANX1, HEPHL1, CP, HEPH, TM4SF1A, VSTM5, C11ORF54, MED17, and SMCO4) within the dp-1 260 kb fine-mapped region (post-CA/NGS: 260,475 bp) were examined for their expression in normal (+/+) chick embryos covering developmental stages HH14 to 27 (Hamburger and Hamilton 1951). Two small nucleolar RNA genes, SNORD6 (63 nt) and SNORA8 (137 nt), could not be assessed using the methodology employed. A probe of at least 300 nt is needed to obtain enough signal to detect a moderate to abundantly expressed RNA. The RNA in situ hybridization procedures were followed as described in Darnell et al. (2007); in some cases RNA probes were created from cDNA (EST) clones (Table 4) acquired from BBSRC ChickEST Database (http://www.chick.manchester.ac.uk/) using methods adapted from Nieto et al. (1996) or by PCR amplification. The expression patterns for these genes at developmental stages other than those shown here can be found at http://geisha.arizona.edu/geisha/ and in Robb (2012). Embryos used were under the care and supervision of trained staff and as per University of Arizona regulations (Protocol # 08-133).

Table 4.

Probes used in whole-embryo in situ hybridization: Analysis of UCD-Diplopodia-1.003 candidate gene expression

Geneab	Accession number	Forward primer (5ʹ–3ʹ)	Reverse primer (5ʹ–3ʹ)	Probe length (bp)
MRE11	NM_204778	TTATAAAAGCTTGCGG CCGCA GAATATACTCC AGCCCAAAAC GAACAG	GCTCTAGAAATTAACCCTCAC TAAA GGGAATGTCAT CATCTGAATCACCCG	931
GPR83	NM_001256138.1	AACAGAAGCCTTG AGGAGGTG	AATTAACCCTCACTAAAGGAA GC GAGGTTTCAGAGGGTG	467
PANX1	XM_015280513.2	GAGAACATTGCTGT GTCCGAG	AATTAACCCTCACTAAAGGAG TCTAAACTCCCAGAAGGATACG	1087
HEPHL1	XM_015280492.2	GTCAGGCTCTAC TACATCGCTG	AATTAACCCTCACTAAAGGG TCAAAT ACATCTCCTCTGTGGC	885
CP	XM_015291853.2	TTCCGAAATAAT GCCAGCCGCC	AATTAACCCTCACTAAAGGC TTC CTCTGCTCTGCTTTTTGG	1000
HEPH	XM_420165.6	GCTGCCATCT ATGAAGTCCGTC	AATTAACCCTCACTAAAGGG TCC CTTTGTTGTCAGGTTGC	1013
TM4SF1A	NM_001277830.1	TCCTGTATTTCC CCAATGGACG	AATTAACCCTCACTAAAGGCC ACACACAAAGCCAATGAGAC	469
VSTM5	XM_003640569.4	CAGCCCAACATCA ACGCAAC	AATTAACCCTCACTAAAGGCC AGT CATTTCAGACAGCATCATC	659
C11ORF54	NM_001277277.1	GTGTTCCATACCTCATACCACTTG	AATTAACCCTCACTAAAGG TCAGCCA TAAGCCCTTCAGTTAC	940
MED17	NM_001006280.1	CCTGTCTCATCAA GAGCAAGAAC	AATTAACCCTCACTAAAGGG TTATTC CATCAACCCCAAGCC	921
SMCO4	NM_001164340.3	TGCTCCCCCATTTTGGTGAAC	AATTAACCCTCACTAAAGGGG CAGA GAGTTTCTTTTGGCT	290
Genea	EST IDc	Clone IDc	BLAT coordinatesd	Probe length (bp)e
MRE11	603507591F1	ChEST434f24	chr1:185695679–185700608	~750
PANX1	603767553F1	ChEST695j13	chr1:185742472–185742920	~850
C11ORF54	603114284F1	ChEST65h21	chr1:185861531–185865290	~900
MED17	603799139F1	ChEST767n9	chr1:185857532–185861163	~900
SMCO4	603851519F1	ChEST848f5	chr1:185946385–185974927	~850

A gene identified in any species within the Dp-1.003-linked 260,475 nt region was used in whole-embryo in situ hybridization*. Note: C11orf54 encodes for protein Fn5; Med17 is a.k.a. Crsp6.

Probes established using PCR.

ESTs were purchased from Source BioScience UK Limited geneservice (Cambridge, UK), through the BBSRC ChickEST database (http://www.lifesciences.sourcebioscience.com/). Clones were selected on carbenicillin plates (50 μg/mL) prior to growth in LB broth+ carbenicillin (50 μg/mL) and clone purification (using Qiagen’s Plasmid Purification Kit). ESTs were sequenced prior to use in RNA ISH to confirm clone identity.

UCSC genome browser (http://genome.ucsc.edu/) was utilized to identify the location of each gene and the EST percent identity to chicken mRNAs previously identified.

Each EST clone was inserted and amplified in the pBluescript II KS+ vector, 3.0 kb (Stratagene). The estimated size of each EST was determined through standard restriction enzyme digest (NotI and EcoRI) and subsequent gel electrophoresis. NotI was used to cleave the vector for sense-strand RNA creation using T3 polymerase. Similarly, EcoRI was used, paired with T7 polymerase to generate anti-sense-strand RNA.

*Two small nucleolar RNA genes found within the region (SNORD6 [63 nt] and SNORA8 [137 nt]) were too small to test by the methodology employed. A probe of at least 300 nt is needed to obtain enough signal to detect a moderate to abundantly expressed RNA by these methods.

Results

Mutant Phenotype

Assessment of the UCD-Dp-1.003 congenic inbred line indicated that ~70% of the mutants (n = 56) display ≤11 total digits on both feet (mode: 5 digits on each foot [n = 27]), while ~30% display ≥12 digits (total for both feet). In 1959, it was reported that 17% of dp-1 mutant embryos displayed asymmetrical numbers of toes; of these, 83% had more toes on the right leg (Abbott 1959b). We found a similar pattern, with 27% of dp-1 mutant embryos displaying an asymmetrical number of digits between the right and left legs. However, of these asymmetrical cases only 33% had more toes on the right leg. Note that the number of digits on the foot (ranging from 4 to 6 on each foot; 9–13 total) was found to be relatively constant within a mating pair. We also observed phenotypic variability in the dp-1.003 mutants with regard to the amount of exposed viscera (Figure 1G–J), degree of dwarfism, and cleft palate (Figure 1K–N), as compared to normal embryos at the same stage of development.

Fine-Mapping Strategies Reduce Mutant-Associated Region Size and Indicate Dp-1.003 GGA 1 Recombination Rate

A 60K SNP genotyping array was initially employed to map the CR for the dp-1 phenotype to GGA1 with a CRmin and CRmax of 610,247 and 707,738 bp (formerly 622 and 720 kb, respectively, based upon the galGal3 assembly) (Robb et al. 2011). In this study the dp-1 mutation was further fine-mapped (region narrowed) using carrier (+/−) individuals and mutant (−/−) progeny via assessment with 8 SNPs (Table 1) to identify their genotypic pattern and detect recombination events. Based upon recombinant analysis of new progeny the CRmin and CRmax were reduced to 104,398 and 585,451 bp, respectively (Table 2).

A breeding scheme was initiated, which included the mating of those heterozygous individuals with a reduced linked region (Table 2), to decrease the CR size in subsequent generations of carrier and mutant progeny. Of the 78 mutant embryos and 106 day-of-hatch chicks (+/−) collected post-60K SNP analysis, 12 and 13, respectively, showed a decreased CR relative to that identified in the 60K SNP array (Table 2). From these results the recombination rate within the GGA 1 CR for the dp-1.003 congenic line was calculated to be ~0.4 cM/Mb (LOD score [Z] = 21.7), thereby indicating high linkage disequilibrium in this region and contrasting with the rate of 1.0–1.5 cM/Mb calculated for this region in other chicken genetic lines (Groenen et al. 2009; Elferink et al. 2010).

Targeted Genomic Capture Enrichment Technology

A targeted-sequence genomic enrichment capture array (CA) paired with NGS was employed to sequence the entire 586 kb region (GGA 1: 185,695,823-186,281,274) linked to the dp-1 mutation by recombinant analysis post-60K array. Bioinformatic analyses were then carried out to identify all dp-1-specific variants (e.g., SNPs, indels). Those elements found within an exon or splice site of known/predicted genes (Table 5) were further assessed for causation in the validation portion of this study. The CA/NGS methods (capture array format and procedures) and results (number of reads, coverage, etc.) for dp-1 are reported in Robb and Delany (2012b). Briefly, a total of 21.0 M dp-1 reads were generated in the CA/NGS sequencing efforts. Those reads covered 96.9% of the 586 kb targeted region with an average fold coverage of 107.2x (Robb and Delany 2012b).

SNP Identification within the Dp-1.003 Congenic Line ~586 kb CR and ~260 kb CR

586 kb CR (GGA 1: 185,695,823-186,281,274) SNPs

A total of 2593 SNPs were found within the GGA 1 targeted region. Multiple pairwise genomic comparisons led to the elimination 1815 SNPs as these variants were shared across one or both of the other developmental congenic lines, which served as controls for GGA 1 on the array (Table 5). The transition:transversion ratio was calculated to be 1.8:1, which is 1.2x lower compared to that reported for other chicken breeds (2.2:1) (Sherry et al. 2001) and compared to the other control congenic lines used on the array (Coloboma.003: 2.2:1; Wingless-2.331: 2.1:1) for this region. Similarly, the average SNP density of the dp-1 region (any SNP identified relative to the RJF reference genome) was 4.4 SNP/kb, which is also 1.2x lower than that reported in other domestic chicken lines (5.1–5.8 SNP/kb) (ICPMC 2004).

A total of 172 SNPs were found to be heterozygous and could be eliminated as causative since mutants are homozygous recessive and require 2 of the same allele to be affected. Assessment of the heterozygous SNP pattern led to the identification of the recombination breakpoints for mutants dp1-166F and dp1-168F (Table 2), which were utilized in the CA because of their reduced linked regions. CA SNP sequence analyses further reduced the CRmax to 260,477 bp (herein referred to as 260 kb) (Figure 2B), specifically moving the coordinate location from chr1: 186,281,274 to chr1: 185,956,300. The new CRmax now spans chromosome 1 from position 185,695,823 to 185,956,300.

Figure 2.

Chromosome 1 CR associated with the chicken diplopodia-1 mutation. A. Causative/Linked Region Identified by the 60K SNP Array: 8 polymorphic SNP markers (red lines) throughout the 707,738 bp CR identified by utilization of the Illumina 60K chicken iSelect SNP genotyping array. SNP markers rs14931758 and rs13990802 denote the boundaries of this region (GGA 1: 185,587,533-186,295,271). A total of 24 genes (green/orange boxes) are encoded within this region. The high-priority candidate gene, MRE11 (Robb et al. 2011), is marked with an orange box. B. Causative/Linked Region Identified by Analysis of the Capture Array Results (CA): A 260,475 bp CR was established through SNP fine-mapping and capture array analyses. Assessment of the capture array data identified a unique dp-1 polymorphic SNP at GGA 1: 185,956,300 (blue line in both panels) which decreased the CR by 447 kb. This region encompasses 13 genes. Note: Coordinates and values shown in the figure represent those of the December 2015 Gallus gallus assembly (galGal5).

~260kb CR (GGA 1: 185,695,823-185,956,300) SNPs

Within the 260 kb CR, 1509 SNPs were identified with 696 of those found to be unique to dp-1 (Table 5) by the pairwise comparison analyses. Of these SNPs, 382 are non-genic and 314 are found within genes. Specifically, 23 SNPs were found in exons or UTRs, 21 at splice sites, and 270 within introns (see Table 5 and Supplementary Table 1). The unique SNPs within an exon or at a splice site were assessed for codon and amino acid changes within 6 reading frames (+1, +2, +3, –1, –2, −3); none generated nonsense mutations (synonymous SNPs: 528; non-synonymous SNPs: 168). Those SNPs found within an exon or splice site of a known gene were further assessed in the validation portion of this study described below.

Micro-indel (1–3 nts) within the UCD-Dp-1.003-Associated CRs

~586 kb CR (GGA 1: 185,695,823-186,281,274) micro-indels

The number and location of micro-indels (insertions or deletions relative to the reference genome) of 1–3 nucleotides in length were investigated. In total, 283 indels (150 insertions, 133 deletions) were placed within the original 586 kb GGA 1 CR with an average distance between micro-indels of 2.11 kb (Table 5). On average, the distance between insertions was 3.98 kb and the average distance between deletions was 4.38 kb. This value (0.48 micro-indels/kb), is slightly higher than that reported previously (~0.4 short indels/kb) (ICPMC 2004). The deletion:insertion ratio (1:1.1) within the 586 kb region for this congenic line is lower than that observed in other breeds for the macrochromosomes (1.6:1) (Brandström and Ellegren 2007). Moreover, the average micro-indel size for Dp-1.003 was 1.3 bp. Note that the average distance between micro-indels and deletion:insertion ratio was calculated here for only the 586 kb region, while values from other studies were calculated for the entire chromosome 1 or macrochromosomes (GGA 1–5) in the case of deletion:insertion ratio.

~260 kb CR (GGA 1: 185,695,823-185,956,300) micro-indels

Within the reduced 260 kb CR, 86 dp-1 specific micro-indels (47 insertions and 39 deletions) were found (Table 5). A total of 31 and 28 insertions and deletions, respectively, were found external to a gene (non-genic), while 21 insertions and 11 deletions were found within a gene (see Table 5 and Supplementary Table 1 for details). Upon frameshift analysis of the 260 kb region, 5 of the 86 micro-indels generated non-synonymous mutations. The functional impact of micro-indels within an exon, UTR, or splice site was also examined. One insertion (+CC), located at the MRE11 donor splice site (3′ end of exon 13), would alter protein structure by failing to remove a 856 nt intron, thereby adding 285 amino acids to the polypeptide and causing a frameshift in the downstream sequence. A second micro-indel (1 nt deletion), present in the 3′ UTR of C11ORF54 could possibly alter the expression of this gene via miRNA-targeted translational inhibition or mRNA degradation. The micro-indels found within an exon, UTR, or splice site of a known gene were further assessed in the validation portion of this study, see below.

Identification of Gaps within the 260 kb CR

Lastly, alignment data indicated 296 gaps (i.e., putative large deletions ≥4 nts) across the 260 kb CR. The average and maximum gap sizes are 73 and 1932 nt, respectively. Those gaps localized to a known exon or splice site (Table 6 and Supplementary Table 1) were validated for legitimacy (i.e., is the gap a large deletion or the result of a sequencing error), for linkage, and potential causation in the validation portion of this study, see below.

Table 6.

Genomic location of dp-1-specific SNPs, micro-indels, and sequence gaps within the 260 kb CR

	SNP locationa			Insertion locationb			Deletion locationc			Gap locationd
Genes	Exon	Splice site	Intron	Exon	Splice site	Intron	Exon	Splice site	Intron	Exon	Splice site	Intron
MRE11	0	3	18	0	1	0	0	0	1	1	0	11
GPR83	1	0	1	0	0	0	0	0	0	1e	0	1
PANX1	0	0	7	0	0	2	0	0	1	1f	0	15
CP	3	0	50	0	0	2	0	0	0	0	0	14
HEPH	1	1	50	0	0	2	0	0	0	0	0	14
HEPHL1	5	1	66	0	0	2	0	0	0	0	0	20
TM4SF1	3e	2	40	0	0	2	0	0	3	1e	0	10
VSTM5	3	0	16	0	0	1	0	0	0	0	0	3
MED17	7e	14	54	0	0	3	0	0	1	3ef	0	23
C11ORF54	3	2	32	0	1	2	1e	1	2	1f	0	13
SNORD6	0	0	0	0	0	0	0	0	0	0	0	0
SNORD8	0	0	0	0	0	0	0	0	0	0	0	0
SMCO4	1e	0	45	0	0	7	0	0	2	1f	0	26

All other unique SNPs (n = 454) are non-genic (not found within a gene).

All other unique insertions (n = 26) are non-genic.

All other unique deletions (n = 27) are non-genic.

All other gaps (n = 554) are non-genic.

One variant is located in the 3′ UTR of the gene.

One variant is located in the 5′ UTR of the gene.

Validation of Elements Identified by Capture Enrichment and Assessment of Causation

It was a goal of this study to use advanced technologies to identify, and without bias, all sequence polymorphisms in the CR. Following that, an equally important aspect was to eliminate those polymorphic elements not involved in the dp-1 mutation representing normal variation in the introgressed region. Finally, it was our goal to identify the top priority candidate elements or genes for future study. As such, the unique SNPs, micro-indels, and gaps found within an exon or splice site (thus potentially impacting mRNA and protein) of known and predicted genes (Table 3) were verified in an additional 20 mutant and 2 control ((UCD-003)) samples to confirm that a particular variant remained linked to the dp-1 mutation. All validated elements (SNPs, micro-indels, variants found in gaps) have been submitted to NCBI (Accessions: ss475871243–ss475871302; ss475875438; ss475871304; ss475871305; ss475875441; ss475875442; ss475871308–ss475871332). Individual variant validation results are described below.

Validation of Exonic, UTR, and Splice Site SNPs, Micro-indels (1–3 nts), and Gaps in the 260 kb CR

A total of 49 SNPs were localized to an exon, UTR or splice site of a known or predicted gene (Table 6, Supplementary Table 1). The additional mutant sample analysis revealed that 41 of the 49 SNPs remain linked to the Dp-1.003 mutation. These variants reside within 7 (HEPHL1, CP, HEPH, TM4SF1A, VSTM5, MED17, and C11ORF54) of the 13 known/predicted genes in the 260 kb region.

Four micro-indels were found within an exon, UTR, or splice site (Table 6, Supplementary Table 1). Assessment of new mutant samples eliminated 2 of these variants—the insertion in MRE11 (+CC) and the deletion in C11ORF54, which upon frameshift analysis would cause a disruption in the polypeptide sequences produced. Thus, 2 micro-indels remain linked to the mutation. One resides within the 3′ UTR of C11ORF54 and the other within a splice site of C11ORF54.

Nine gaps (>4 nt) were identified within an exon, UTR, or splice site (Table 6, Supplementary Table 1). Of these, one gap identified no longer is localized to chromosome 1 based upon the new genome assembly update. Sequencing of the remaining 8 gaps identified the presence of sequence in 3 gaps identical to control DNA; thus, these genome gaps were likely due to a sequencing error during completion of the CA. Furthermore, validation efforts reliably identified sequence in 3 gaps; however, upon DNA alignment, 6 SNPs and one insertion reside within these regions. These newly identified variants remain linked to the Dp-1.003 mutation and are found within an exon and splice site of MED17 and within an exon of C11ORF75. Moreover, 2 gaps remain linked due to sequencing validation complications attributed to the high G-C percentage and polymerase slippage during sequencing. These gaps (aka ≥4 nt deletions) are found within the 3′ UTR of PANX1 and the 5′ UTR/splice site of C11ORF54 (Supplementary Table 1).

Developmental Gene Expression of Candidates by RNA In situ Hybridization

To assist in our prioritization of the likelihood that a gene is involved in the dp-1 mutation, the expression patterns of 11 protein-coding genes (MRE11, GPR83, PANX1, HEPHL1, CP, HEPH, TM4SF1A, VSTM5, C11ORF54, MED17, and SMCO4) mapping within the 260 kb CR were studied during early embryonic stages in normal (+/+) chicken embryos. Eight of these genes did not show expression during early embryonic developmental stages (Figure 3D–K), while one gene (PANX1; Figure 3C) displayed a positive staining pattern of expression but not in the regions of interest (e.g., limbs, pharyngeal arches/clefts, facial prominences). Interestingly, 2 genes, MRE11 and GPR83, which are less than 5000 bp apart, show positive expression results in which appropriate, localized expression was observed during stages of embryogenesis wherein organogenesis, limb, and craniofacial development occur (Figure 3A,B).

Figure 3.

Priority candidate genes MRE11 and GPR83 show strong expression in structures leading to the formation of the face and limbs in normal chick embryos. RNA in situ hybridization was utilized to assess the expression of 11 candidate protein-encoding genes in the 260 kb region using standard procedures. Analysis of the results identified 3 genes with positive expression (A–C), with 2 (A,B) showing appropriate dp-1 spatio-temporal expression and the third showing expression in regions not affected in the dp-1 mutant embryo. Additionally, 8 genes showed negative expression (D–K). A. MRE11 (HH24): high levels of expression in the limb buds, somites, pharyngeal arches/clefts, and brain. B. GPR83 (HH24): punctate expression in the limb buds, somites, and pharyngeal arches/clefts. C. PANX1 (HH22) expression is observed in the ventral neural tube (not shown: pancreatic expression identified via ventral cross-section). D.HEPHL1 (HH24), only background is present. E.CP (HH25), only background staining is present. F.HEPH (HH24), only background staining is present. G.TM4SF1A (HH24), only background is present. H.VSTM5 (HH24), only background staining is present. I.C11ORF54 (HH23), only background staining is present. J.MED17 (HH24), only background staining is present. K.SMCO4 (HH24), only background staining is present. L. SPRY1 (HH24)—positive control: strong, highly specific expression in the limb buds, somites, pharyngeal arches/clefts, and brain. M. TNNI1 (HH24)—positive control: strong expression in the somites, atrium/ventricle, and pharyngeal arches. N. No probe (HH24)—negative control (no gene probe used): controls for background staining that may arise independent of the presence of an antisense probe. The expression patterns for these genes at developmental stages other than those shown in this figure can be found at http://geisha.arizona.edu/geisha/.

The dp-1 mutant phenotype displays abnormalities within several main anatomical features (i.e., somites, mandibular arch, head and limb regions). At HH24, a stage in which the limb buds and maxillary arches are very distinct, strong MRE11 expression is found within the developing limbs, maxillary process, mandibular and hyoid arches, optic cup, somites, ventral neural tube, and brain (mesencephalon and telencephalon) (Figure 3A). Additionally, GPR83 expression at the same stage can be clearly observed in the limbs, maxillary process, mandibular and hyoid arches, optic cup, somites, ventral neural tube, ventricle of the heart and brain (telencephalon) (Figure 3B). The expression pattern for these genes at additional developmental stages can be found at http://geisha.arizona.edu/geisha/.

Discussion

The chick embryo shares many features (anatomy, function, organization, cellular composition, and molecular pathways) in common with higher vertebrates, including mammals, during development (Le Douarin 2004; Tickle 2004; Stern 2005). This allows for valuable comparison between similar chick (e.g., dp-1) and human malformations (e.g., diplopodia). Despite the advanced state of the human genome sequence, not all heritable disorders have been linked to specific genes. Therefore, the chick embryo, specifically the developmental mutant dp-1 described herein, is an excellent system in which to search for and identify a causative element, with the hope that this information can then be applied and developed further for screening and therapeutic applications of human malformations. Moreover, such animal model research adds knowledge to known developmental pathways and can uncover new molecular networks and gene involvement in normal vertebrate development.

Our aim was to generate new information on the genomics of the candidate region as well as the molecular, developmental, and morphological processes that contribute to the dp-1 phenotype. Dp-1 mutants develop normal digits (digits 2 through 4 are normal on the foot; digits 2 and 3 on the wing are normal); however, the most anterior digit, digit 1, is accompanied by 2–3 digit supernumeraries (Figure 1) (MacCabe and Abbott 1974; MacCabe et al. 1975). The dp-1 mutant also exhibits truncated limbs (micromelia), dwarfism, exposed visceral organs, and cleft palate. These conditions are seen individually or combined in reported cases of human congenital disorders (reviewed in Vlahovic et al. 2015). Interestingly, the extent of phenotypic variation observed in mutants of the dp-1.003 congenic inbred line was reduced compared to initial reports on the dp-1 mutation. We hypothesize that through the establishment of the congenic line (breeding the dp-1 mutation on the highly inbred UCD-003 background, F > 0.999), variant sequence elements (e.g., promoter, epistatic influences, etc.) contributing to the phenotypic variability observed in the 1940s through 1970s were eliminated, thereby creating the relatively uniform phenotype currently observed.

To date, dp-1 has been well-studied as to processes involved in normal and mutant limb development (Taylor and Gunns 1947; Landauer 1956; Abbott 1959a, 1959b; Taylor et al. 1959; MacCabe and Abbott 1974; Taylor 1974; MacCabe et al. 1975; Rodriguez et al. 1996). Tissue interaction studies revealed the defect to be of mesodermal origin (MacCabe and Abbott 1974; MacCabe et al. 1975). Rodriguez et al. (1996) studied several genes in the dp-1 mutants known to be involved in limb formation and digit outgrowth (SHH, HOXD, BMP2, FGF4). The expression patterns for SHH and HOXD were similar to normal embryos. However, BMP2 and FGF4 expression was expanded. Noteworthy, none of these genes reside within the dp-1-linked CR of GGA 1 (this study and Robb et al. 2011); thus, the causative element of dp-1 may be operating within a novel pathway or is a critical upstream element in a known pathway.

In this study, we assessed the RNA expression patterns of 11 protein-coding genes (MRE11, GPR83, PANX1, HEPHL1, CP, HEPH, TM4SF1A, VSTM5, C11ORF54, MED17, and SMCO4) within the 260 kb fine-mapped region. Of these, 2 proximally close genes (MRE11 and GPR83) showed strong expression (Figure 3A,B) during the appropriate stages of development wherein main morphological features (limbs, facial prominences, pharyngeal arches/clefts, somites) affected in the dp-1 mutation start to and become elaborated. To determine if the function of these genes is conserved across vertebrates, we performed Mre11 RNA in situ hybridization in mouse and obtained similar results (Supplementary Figure 1). Histological analyses reported for GPR83, in mouse are also consistent with what was observed here (http://www.informatics.jax.org/marker/MGI:95712). In other model organisms, such as the zebrafish, gene expression was reported in the pectoral fin buds and pharyngeal arches for mre11 (Thisse and Thisse 2004); no studies on gpr83 in the zebrafish have been published to date. Moreover, human fetuses at 10–20 weeks of gestation shows elevated MRE11 and GPR83 expression in the same tissues (NCBI: https://www.ncbi.nlm.nih.gov/gene; geneID: 4361 and 10888). Although a gene might not be expressed at these developmental stages (as in the case of the other 9 non-expressing genes), hypothetically, a polymorphism could induce (mis)expression. It is noteworthy, however, that even at later stages of development, the non-expressing genes found in this study are also not expressed in the mutant-relevant tissue locations (NCBI: https://www.ncbi.nlm.nih.gov/gene; European Bioinformatics Institute Expression Atlas: https://www.ebi.ac.uk/gxa/home; The Protein Atlas:https://www.proteinatlas.org/).

The primary goal of our genomics research was to discover and establish those variants, which might be involved in dp-1 through sequencing, linkage analysis/fine-mapping and variant validation in additional mutants or invalidation in non-mutants. Interestingly, splice and/or exonic variants remain linked in 8 genes (PANX1, HEPHL1, CP, HEPH, TM4SF1A, VSTM5, C11ORF54, MED17) but not in 5 genes (MRE11, GPR83, SNORD6, SNORA8, SMCO4). Although exonic and splice variants were eliminated as being causal in both MRE11 and GPR83, a dp-1-linked regulatory element could impact their gene expression (i.e., cis-regulatory element) and/or modify genes in other molecular pathways (i.e., trans-regulatory element, e.g., altering BMP2 or FGF4 expression). This type of gene action is possible as mutations in regulatory elements commonly affect the expression pattern of a gene and notably has already been reported in the dp-1 mutant with the expanded expression of BMP2 and FGF4 (Rodriguez et al. 1996).

Based on comparative vertebrate biology (mouse, human, zebrafish), pathway role (Supplementary Figure 2), and the spatio-temporal expression profiles in this work we suggest MRE11 as a candidate of high interest for the dp-1 mutation. Until recently, GPR83 had only been identified for its role in peripheral T-cell development (Lu et al. 2007). However, our RNA in situ results (appropriate spatio-temporal expression, see Figure 3B) and recent functional studies (Müller et al. 2013) combined suggest an additional role of GPR83 in limb and craniofacial development, and a potential role in the dp-1 mutant phenotype. Below we outline more on the function of MRE11 and GPR83, associated human diseases, molecular pathway(s) of interest, as well as the next steps and implications of this work.

MRE11’s Potential for a Contribution to the Dp-1 Phenotype

MRE11 (meiotic recombination 11 homologue) is a protein with known involvement in DNA double-stranded break repair, homologous recombination, nonhomologous end joining, and telomere length maintenance (Lee and Paull 2005; Deng et al. 2009; Zha et al. 2009). Altered regulation and mutations in MRE11 were shown to produce human pleiotropic disorders (e.g., Nijmegen breakage syndrome (Matsumoto et al. 2011), Ataxia Telangiectasia-like disorder (Stewart et al. 1999)) with affected developmental systems similar to that observed in dp-1 mutants. Moreover, mutations in genes with similar functions to MRE11 (e.g., NIPBL, SMC1L1) are associated with Cornelia de Lange syndrome (CdL) (Krantz et al. 2004; Tonkin et al. 2004; Kaur et al. 2005; Musio et al. 2006; Barber et al. 2008), which shares features in common with the dp-1 phenotype such as limb and digit abnormalities, growth retardation, cleft palate, and organ malformations (Jackson et al. 1993; Kline et al. 2007; Liu and Krantz 2009).

Evaluation of downstream targets of MRE11 suggest 2 potential molecular mechanisms that could contribute to an aberrant phenotype such as dp-1 given potential MRE11 misexpression; these are highlighted in more detail in Supplementary Figure 2 (and corresponding Supplementary Material). Briefly, under normal circumstances MRE11 recruits ATM (a serine/threonine kinase) to a double-stranded break to initiate repair of DNA damage. In turn, ATM subsequently activates DNMT1 (Uziel et al. 2003), a DNA methyltransferase that preferentially methylates specific genomic regulatory regions, including that of the NOGGIN promoter (Shamma et al. 2013). Interestingly, NOGGIN is an antagonist of BMPs which bind directly to BMP receptors (Supplementary Figure 2A) and, as a result, plays an important role in many developmental processes for limb, craniofacial, and skeletal formation (Reshef et al. 1998; Wijgerde et al. 2005; Bayramov et al. 2011) by creating morphogenic gradients. Thus, misregulation of NOGGIN via DNMT1 methylation can lead to irregular BMP spatial and temporal expression (thereby altering craniofacial and limb development) (outlined in Supplementary Figure 2B,C), as was reported by Rodriguez et al. (1996) in their studies of dp-1 mutants (Figure 1B–F, K–N). Additionally, ATM phosphorylates histone H2A variant, H2AX (Stiff et al. 2006), thereby inducing chromatin remodeling and subsequent apoptosis/senescence of the cells in the interdigit space of vertebrate embryos (Montero et al. 2016). Failure to trigger programed cell death results in interdigital webbing (Supplementary Figure 2A), another phenotype observed in dp-1 mutants (Figure 1B–F). Thus, given these genetic networks, as well as the MRE11 mutant human developmental phenotypes, we hypothesize that abnormal expression of MRE11 (and as a result, effects on downstream pathways) is a likely contributor of the dp-1 phenotype.

GPR83’s Potential for Contribution to the Dp-1 Phenotype

GPR83 (G Protein-Coupled Receptor 83) is an orphan protein with an identified, but not well-understood role in regulatory T cell development via its interaction with Foxp3 (Lu et al. 2007). Additional expression analysis in the rat brain led researchers to propose a second role for GPR83 in central processing of glucocorticoids (Sah et al. 2005). Specifically, GPR83 regulates metabolic pathways by binding to ghrelin’s receptor (Ghsr1a) (Müller et al. 2013). When Gpr83 is knocked out, mice show lower body weight, fat mass, food intake, glucose tolerance, and insulin sensitivity, as compared to control mice on the same diet (Müller et al. 2013), phenotypes which are consistent with that of a knockdown/out ghrelin receptor. These results therefore suggest that altered/mutated Gpr83 disturbs native Ghsr1a function.

Interestingly, loss of Ghsr1a in mouse leads to short stature through reduced hypothalamic–pituitary–growth hormone activity (a phenotype observed in the dp-1 mutant, Figure 1) as well as a hypoglycemic state (Pantel et al. 2006; Sun et al. 2006, 2008). Interestingly several studies outline the presence of both cleft lip/palate in children with neonatal hypoglycemia (De Leon and Stanley 2017). Moreover, neonatal hypoglycemia has been associated with Beckwith–Wiedemann syndrome, characterized by in midline abdominal wall defects (omphalocele) (Queensland Clinical Guidelines 2013) similar to that observed in the dp-1 mutant (Figure 1G–J), as well as dysmorphic facial features (https://ghr.nlm.nih.gov/). Note that other congenital abnormalities are also found in neonatal hypoglycemic patients that are not observed in the dp-1 mutants. Thus, based upon our spatio-temporal in situ expression results and comparative mouse/human studies, we suggest that GPR83 should also be considered a priority candidate gene.

Final Thoughts and Research Implications

Although Mre11 knockout results in early lethality (Buis et al. 2008) and Gpr83 knockout mice survive to adulthood with no morphological abnormalities (Müller et al. 2013), our proposed hypothesis (misexpression via a polymorphism in the 260-kb region) of one, or perhaps even both genes deserve future analysis especially with regard to intronic elements within MRE11 (18 SNPs, 1 indel) and GPR83 (1 SNP) that remain linked to the phenotype in addition to a SNP in the intergenic region between these 2 proximally close genes.

Performing such research is of course a complicated effort. Moreover, the proximity of the 2 high priority candidate genes (5000 bp apart), makes it even more difficult (and yet interesting given the mode of inheritance) to further eliminate variants using breeding techniques and recombination/linkage analyses. And, although much has been discovered regarding gene regulation, status/knowledge of regulatory motifs and their impacts remains obscure. Functional studies will be essential to explore the remaining variants and their contributing role in MRE11 and GPR83 expression toward the dp-1 phenotype.

Through this work we successfully discovered unique genetic variants specific to the dp-1 mutation, yet despite the advances of genomics and allied fields, this work affirms what other such research has shown, that finding the exact sequence element responsible for even a simply-inherited trait (single-gene autosomal recessive pattern) is not easily accomplished. In addition to more studies focusing on MRE11 and GPR83, evaluation of elements found within the 260 kb-linked region, beyond those found within exons, UTRs and splice sites reported here, will be a contribution to uncovering the causative sequence. Identification of the causal mutation for dp-1 will be valuable to not only the field of developmental biology, shedding light on the mechanisms involved in normal/abnormal development, but also to the fields of genetics and genomics by elucidating the pathways involved in the regulation and expression of various genes and involved pathways. Ultimately, the application of this research in the medical field and the implementation of unique genetic lines as biomedical models to study cleft palate, dwarfism, polydactyly, and omphalocele, could lead to screening systems for congenital malformations as well as research opportunities for therapeutic remedies.

Funding

This work was supported by the U.S. Department of Agriculture-National Institute of Food and Agriculture Multistate Research NC-1170 (CA-D*-ASC-6414-RR), the NRSP-8 National Animal Genome Research Support Program (CA-D*-ASC-7233-RR) and the UC Davis (UCD) John and Joan Fiddyment endowment (M.E.D.). The in situ hybridization studies were supported by National Institutes of Health grant HD064559 (P.B.A.). We gratefully acknowledge the UCD Department of Animal Science for graduate student fellowship support and the UCD Jastro-Shields Graduate Research Award (E.A.O.). We appreciate the infrastructure and the poultry genetic resources supported by the UC Davis Department of Animal Science, College of Agricultural and Environmental Sciences, and the California Agricultural Experiment Station.

Data Availability

We have deposited the primary data underlying these analyses as follows:

-Validated sequence variants: NCBI accessions: ss475871243-ss475871302; ss475875438; ss475871304; ss475871305; ss475875441; ss475875442; ss475871308–ss475871332.

-RNA in situ hybridization gene images: GEISHA—Gallus Expression in Situ Hybridization Analysis (http://geisha.arizona.edu/geisha/).

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.