Simultaneous Occurrence of Hypospadias and Bilateral Cleft Lip and Jaw in a Crossbred Calf: Clinical, Computer Tomographic, and Genomic Characterization

Simple Summary Congenital abnormalities in animals are a major concern for breeders due to the increased economic loss they entail. The aim of this article was to describe a congenital bilateral cleft lip and jaw and an abnormal opening of the penile urethra in a crossbred calf. Clinical examination, computed tomography, and whole genome sequencing were performed to describe and identify a possible cause of the abnormalities. The whole genome investigation indicates the involvement of multiple genes in the two birth defects in this case. Abstract Congenital abnormalities in animals, including abnormalities of the cleft lip and jaw and hypospadias have been reported in all domesticated species. They are a major concern for breeders due to the increased economic loss they entail. In this article, we described a congenital bilateral cheilognathoschisis (cleft lip and jaw) with campylognathia in association with penile hypospadias and preputial hypoplasia with failure of preputial fusion in a Bos taurus crossbred Piedmontese × Wagyu calf. Clinical examination, computed tomography, and whole genome sequencing were performed to describe and identify a possible cause of the abnormalities. Clinical examination revealed a bilateral cheilognathoschisis of approximately 4 cm in length and 3 cm in width in the widest part, with computer tomography analyses confirming the bilateral absence of the processus nasalis of the incisive bone and the lateral deviation of the processus palatinus towards the left side. Genomic data analyses identified 13 mutations with a high impact on the products of the following overlapped genes: ACVR1, ADGRA2, BHMT2, BMPR1B, CCDC8, CDH1, EGF, F13A1, GSTP1, IRF6, MMP14, MYBPHL, and PHC2 with ADGRA2, EGF, F13A1, GSTP1, and IRF6 having mutations in a homozygous state. The whole genome investigation indicates the involvement of multiple genes in the birth defects observed in this case.


Introduction
Congenital oral cavity anomalies are disturbances or failures of the complex processes of craniofacial development. The stomodeum, a primitive oral cavity is formed from five processes: frontonasal, medial nasal, lateral nasal, maxillary, and mandibular, derived from the neural crest mesenchyme and ectomesenchyme [1]. The oral cavity is formed by initial fusions between the medial and lateral nasal processes, followed by fusions between the lateral nasal and maxillary processes, and finally, the maxillary processes merge with the mandibular processes [2][3][4][5].
Failures in the fusion of the lateral and medial nasal processes produce cheiloschisis (lip cleft) and a failure of the medial fusion of the palatal shelves determines palatoschisis (cleft palate). Many molecular signals such as BMP (bone morphogenetic proteins), FGFs (fibroblast growth factors), SHH (Sonic Hedgehog), and WNT (Wingless-int) genes, are involved in the growth and fusion of facial processes [6], and mutations of a number of genes, e.g., BMP4, BMPR1A, SOX11, WNT9B [7], or FGF8 [8] affect the normal development of facial processes.
Facial clefts, such as congenital lip and jaw (CLJ), can be classified depending on their location into unilateral right or left side, medial or bilateral CLJ [2,9,10], with a variable expression of the phenotype, both in humans [11][12][13][14] and animals [2,10]. An additional criterion is whether they are associated or not with malformations in other organs: syndromic or non-syndromic forms.
Possible causes for the non-syndromic forms are multifactorial diseases as a result of the interaction between genetic and environmental factors. They do not follow a monogenic pattern of inheritance, but they reappear in future generations [2]. Syndromic forms, in many cases, are due to aberrations of chromosomes or are associated with monogenic diseases [19].
In cattle, the frequency of congenital defects is estimated to be 0.25% to 3.6% and can be caused by genetic factors, environmental factors, or genetic-environmental interaction during embryonic development [9,24,25].
Numerous environmental factors have been reported to be associated with facial abnormalities during embryonic development, such as chemicals, toxic substances, infections, hormonal factors etc. [2,4,9].
The other congenital abnormality seen in the calf was hypospadias (OMIA 001187-9913), a developmental abnormality in which the male urethra opens on the ventral side of the penis or on the perineum and can be accompanied by a failed development of the median perineal raphe, and of scrotum, penile aplasia or hypoplasia, hypoplasia of the corpus cavernosum, and incomplete formation of the prepuce [26][27][28]. Abnormal development of the penile urethra is due to a malfusion of the urogenital folds in the midline of normal penile development [26,28]. Depending on where the fusion of the urogenital folds failed, hypospadias can be classified as balanitic/glandular-the urethra opens in the region of the glans penis, penile-the urethra opens ventral and caudal to the glans penis and cranial to the scrotum, and could be proximal, distal, or in the mid shaft of the penis, scrotal-the urethra opens between the halves of the non-fused scrotum or at the scrotum, perineal-the urethra opens in the perineum, ventral to the anus, or anal-the urethra opens in the anal region [24,29,30]. In animals, this abnormality occurs relatively rarely in horses [31], sheep, goats, cattle [24,32,33], dogs, and cats [34,35].
In the etiology of hypospadias, there are incriminated genetic, endocrinological, environmental factors and even epigenetic alterations of the genetic material [25,28,38]. Data from the literature indicate that between 10-30% of human hypospadias cases can be at-tributed to specific gene mutations (e.g., AR gene mutation, SRD5A type II deficiency), while 70% of cases have unknown etiology, possibly a combination of genes and environmental factors [38][39][40]. Abnormal expression patterns of some transcription factors, with a key role in the spatiotemporal regulation of gene expression, may disrupt the biological processes during embryogenesis [41].
The morphogenesis of the genital tubercle involves different genes such as BMP2, BMP4, BMP7, FGF, FGF10, FGF8, HOXA13, HOXA4, HOXB6, HOXD13, TGF, and WNT5A genes [41][42][43]. Moreover, ESR1 and ESR2 genes through estrogens play a role in normal penile and clitoral development, and normal urethral development depends on the balance between androgens and estrogens [44][45][46]. Studies carried out on mice showed that the exposure of pregnant females to external estrogens induced a thin periurethral spongiosa in the offspring [47] and the loss of endogenous estrogen signaling in mice can cause a mild hypospadias in the offspring [44].
Additionally, epigenetic alterations can be involved in hypospadias development such as the methylation of the AR gene [48], the SRD5A2 gene [49], or CpG sites [28], determining the use of DNA methylation patterns to identify and evaluate new candidate genes that may be involved in the etiology of hypospadias [28].
The purpose of this study was to investigate a male Piedmontese × Wagyu calf with congenital bilateral cheilognathoschisis and hypospadias using computer tomography analysis and whole genome sequencing in order to identify the possible cause of those congenital defects. To our knowledge, based on literature data, we present a unique case of congenital defects in a Bos taurus hybrid Piedmontese × Wagyu crossbred cattle using different paraclinical examinations.

Statement of Ethics
The case study was performed in compliance with national (Romanian Law 43/2014) and European guidelines (EU Directive 63/2010) for using animals for scientific purposes, and in accordance with the internal SOP (Operational Procedures System) of USAMVBT-PG-001-R021 (Banat's University of Agricultural Sciences and Veterinary Medicine) and Ethical Statement no. 73/2020.

Animal
The male Bos taurus hybrid Piedmontese×Wagyu crossbred calf with a body weight of 32 kg was born alive after a normal pregnancy period of a heifer in September 2020, in the Didactic Farm of the University of Life Sciences ''King Mihai I" from Timişoara, Romania. Its mother was a primiparous Piedmontese crossbred (a product between a Holstein female and a Piedmontese male) and the father was from the Wagyu breed. It was the first case with such congenital abnormalities in our Research Farm and, to our knowledge, the first such documented case in the world. The calf was euthanized a few days later due to severe neonatal diarrhea and anemia, as clinically it was impossible to feed properly with colostrum and it could not be kept in a quadrupedal position. Euthanasia was performed with Ketamidor 100 mg/mL from Richter Pharma, Austria, and Xylazin Bio 20 mg/mL from Bioveta, Cech Republic and followed by T61 (embrutramide 200 mg/mL, mebezonium iodide 50 mg/mL, tetracaine hidrochloride 5 mg/mL) 1 mL/10 kg, intravenous, from Intervet International BV Wim de Korverstraat, Holland.

Computer Tomography (CT) Analyzes
The acquisition of CT scan images was performed postmortem with the body positioned such that the ventral aspect of the head was in contact with the CT table (ventral recumbency). CT scans were performed with a Siemens Somatom Definition AS 64 scanner, using conventional settings (90 kV, 100 mAs), and a slice thickness of 0.6 mm.

Whole Genome Sequencing
Genomic DNA was isolated from whole blood collected from the affected calf using a DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) for DNA extraction. Whole-genome sequencing was performed by Novogene, Cambridge, United Kindom. The sample was 41.13 µL in volume, 10.209 ng/µL concentration with a total amount of DNA of 0.41989 µg. The DNA concentration was measured using a Qubit ® DNA Assay Kit on a Qubit ® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA).
Next, the library preparation was performed. Briefly, the genomic DNA was randomly fragmented into 350 bp fragments, and then DNA fragments were end polished, A-tailed, and ligated with the adapters for Illumina sequencing, and further PCR enriched with P5 and P7 oligo primers. The PCR products as final construct of the libraries were purified (AMPure XP system), followed by a quality control test that included the size distribution using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and the measurement of molarity using real-time PCR.
Raw sequenced reads in FASTQ format were analyzed using the Bcbio-nextgen pipeline v. 1.2.9. The quality control of the raw data was performed with FastQC v. 0.11.p [50]. Read decontamination was performed with kraken2 v. 2.1.2 [51] using the options "-minimum-hit-groups 4 -confidence 0.5" and a newly-created database consisting of all the genomic sequences found in the NCBI RefSeq database [52] as of May 2022 for human, bacteria, archaea, viruses, fungi, and protozoa. The databases consisted of 101,535 sequences with a summed length of 138.7 Gbp and occupied approx. 117 GB of storage space. The decontaminated reads were mapped to the ARS-UCD1.2 reference genome [53] using the BWA software package v. 0.7.17 [54]. The FASTA file with genomic sequences and the corresponding GTF file with genomic locations for the ARS-UCD1.2 reference were acquired from Ensembl, release 105 [55]. The variant calling was performed with GATK HaplotypeCaller v. 4.2.6.1 [56]. Low-complexity regions consisting of microsatellites and simple repeats were downloaded from UCSC for the ARS-UCD1.2 reference and excluded from the variant calling process. Variant effects on gene transcripts were predicted using the Ensembl Variant Effect Predictor (VEP) tool v. 105.0 [57] and the Ensembl GTF file mentioned above.
The variants classified as having a high or moderate impact and the corresponding genes were retained for the downstream analysis, which was performed using various packages for the R programming language v. 4.1.3 [58]. A list of 126 genes implicated in hypospadias was put together following a survey of the relevant literature, as previously described [27,30,40,59]. Similarly, a list of 513 genes known to be involved in the occurrence of cheilognathoschisis (cleft lip and jaw) was obtained from the CleftGeneDB database [60], as of July 2022. After eliminating duplicate genes that were associated with hypospadias and cleft lip and jaw, a set of 559 unique target genes was retained. Finally, gene annotations were obtained via the AnnotationHub package v. 3.2.2 [61] and multiple over-representation analyses based on GeneOntology (GO) [62], Kyoto Encyclopedia of Genes and Genomes (KEGG) [63], Medical subject headings (MeSH) [64], and Reactome [65] were performed with the R packages clusterProfiler v. 4.2.2 [66], meshes v. 1.20.0 [67], and ReactomePA v. 1.38.0 [68]. In order to perform the previously mentioned over-representation analyses, common gene identifiers reported by Ensembl VEP for high impact and moderate impact variants were searched within the corresponding OrgDb and MeSHDb annotation databases from AnnotationHub for Homo Sapiens, using the select function from the AnnotationDbi package v. 1.56.2 [69]. Gene annotation was performed on sets of all genes associated with high-impact and/or moderate-impact variants (i.e., not only for the 559 target genes), in order to find significantly associated annotation terms with the specific variants found for the Bos taurus individual from this study.

Clinical Findings
At birth, the affected calf presented a cleft lip and jaw (OMIA 001714-9913) and hypospadias (OMIA 001187-9913). The bilateral cleft was characterized by a length of approximately 4 cm on each side and a width of 3 cm on the widest part. On the andrological examination, the calf showed penile hypospadias, the penile urethra ventrally opened at approximately 4 cm caudal of the normal opening of the urethra and the prepuce was incompletely developed and had failed to fuse normally (Figure 1a-d).
Animals 2023, 13, x FOR PEER REVIEW 5 of 18 database [60], as of July 2022. After eliminating duplicate genes that were associated with hypospadias and cleft lip and jaw, a set of 559 unique target genes was retained. Finally, gene annotations were obtained via the AnnotationHub package v. 3.2.2 [61] and multiple over-representation analyses based on GeneOntology (GO) [62], Kyoto Encyclopedia of Genes and Genomes (KEGG) [63], Medical subject headings (MeSH) [64], and Reactome [65] were performed with the R packages clusterProfiler v. 4.2.2 [66], meshes v. 1.20.0 [67], and ReactomePA v. 1.38.0 [68]. In order to perform the previously mentioned over-representation analyses, common gene identifiers reported by Ensembl VEP for high impact and moderate impact variants were searched within the corresponding OrgDb and MeSHDb annotation databases from AnnotationHub for Homo Sapiens, using the select function from the AnnotationDbi package v. 1.56.2 [69]. Gene annotation was performed on sets of all genes associated with high-impact and/or moderate-impact variants (i.e., not only for the 559 target genes), in order to find significantly associated annotation terms with the specific variants found for the Bos taurus individual from this study.

Clinical Findings
At birth, the affected calf presented a cleft lip and jaw (OMIA 001714-9913) and hypospadias (OMIA 001187-9913). The bilateral cleft was characterized by a length of approximately 4 cm on each side and a width of 3 cm on the widest part. On the andrological examination, the calf showed penile hypospadias, the penile urethra ventrally opened at approximately 4 cm caudal of the normal opening of the urethra and the prepuce was incompletely developed and had failed to fuse normally (Figure 1a   The anus was open physiologically. The examination of the scrotum revealed that the testicles were not fully descended and at the necropsy, the testicles with a length of Animals 2023, 13, 1709 6 of 18 2.2 cm were identified in the inguinal canal. The necropsy revealed that there was no other internal congenital abnormality in addition to the previously mentioned external defects.

Computed Tomography Findings
CT of the head ( Figure 2) revealed a bilateral cleft lip and jaw, with the bilateral absence of the processus nasalis of the incisive bone (Figure 2a-c). Only the processus palatinus was identified as being completely developed bilaterally, but both were deviated laterally toward the left side. A complete absence of the lateral aspects of the dental plate was observed bilaterally (Figure 2d).
The anus was open physiologically. The examination of the scrotum revealed that the testicles were not fully descended and at the necropsy, the testicles with a length of 2.2 cm were identified in the inguinal canal. The necropsy revealed that there was no other internal congenital abnormality in addition to the previously mentioned external defects.

Computed Tomography Findings
CT of the head ( Figure 2) revealed a bilateral cleft lip and jaw, with the bilateral absence of the processus nasalis of the incisive bone (Figure 2a-c). Only the processus palatinus was identified as being completely developed bilaterally, but both were deviated laterally toward the left side. A complete absence of the lateral aspects of the dental plate was observed bilaterally (Figure 2d). Based on cleft lip-jaw scoring system in ca le existing in literature where different orofacial structures, such as lips, the processus nasalis of the os incisivum, the dental plate with adjacent segments of the hard palate, the facial bones-os incisivum, os maxillare, os nasale, os palatinum, and the mandibles, we scored the cleft lips and accompanying malformations in the face in this case as fallow. The score was 4 (on a scale from 1-4) for the site (bilateral left and right); a score of 5 (on a scale from 0-5) for the depth of the cleft lip (complete cleft of the upper lip independent of whether the maxillary alveolus is affected), a score of 6 (on a scale from 1-7) for changes at the dental plate and hard palate (about 2/3 of the dental plate and adjacent segments of the hard palate were missing), a score of 1 (on a scale from 0-3), for curvature to the right (lateral curvature of the mandibles), and a score of 5 (on a scale from 0-9) for changes in the processus nasalis of the os incisivum (absence of the apical segment of the processus nasalis).

Whole Genome Sequencing Findings
Statistics regarding the proportion of reads mapped to the reference genome were computed. There was a total of 164.9 million mapped reads, a proportion of 98% of the total of reads. The percentage of PCR duplicate reads was 29%. The percentage of 1X coverage was 93.9%, for 5X coverage was 45%, for 10X coverage was 14.1%, with the mean genome coverage of 6.3X, as shown in Supplementary Figure S1 Based on cleft lip-jaw scoring system in cattle existing in literature where different orofacial structures, such as lips, the processus nasalis of the os incisivum, the dental plate with adjacent segments of the hard palate, the facial bones-os incisivum, os maxillare, os nasale, os palatinum, and the mandibles, we scored the cleft lips and accompanying malformations in the face in this case as fallow. The score was 4 (on a scale from 1-4) for the site (bilateral left and right); a score of 5 (on a scale from 0-5) for the depth of the cleft lip (complete cleft of the upper lip independent of whether the maxillary alveolus is affected), a score of 6 (on a scale from 1-7) for changes at the dental plate and hard palate (about 2/3 of the dental plate and adjacent segments of the hard palate were missing), a score of 1 (on a scale from 0-3), for curvature to the right (lateral curvature of the mandibles), and a score of 5 (on a scale from 0-9) for changes in the processus nasalis of the os incisivum (absence of the apical segment of the processus nasalis).

Whole Genome Sequencing Findings
Statistics regarding the proportion of reads mapped to the reference genome were computed. There was a total of 164.9 million mapped reads, a proportion of 98% of the total of reads. The percentage of PCR duplicate reads was 29%. The percentage of 1X coverage was 93.9%, for 5X coverage was 45%, for 10X coverage was 14.1%, with the mean genome coverage of 6.3X, as shown in Supplementary Figure S1.
Variant calling identified a total number of 5,598,215 variants, of which 4,951,834 were SNPs and 648,459 were indels. Predicted consequences of all variants (most severe per variant), computed by Ensembl Variant Effect Predictor (VEP), are presented in Figure 3, with the exception for upstream-gene-variants, which were 3,525,429, and intron-variants with 1,840,076 incidence.
3, with the exception for upstream-gene-variants, which were 3,525,429, and intron-variants with 1,840,076 incidence. Most of the variants, such as non-coding transcript variants, intergenic variants, synonymous variants, are considered as unlikely to change protein behavior or with low or no impact on it. Conversely, as can be seen in Figure 3, Ensembl VEP identified a low number of variants that are assumed to have a high impact in the protein (transcript ablation, splice acceptor variant, stop gained, or stop lost).
An initial quality control checks on raw sequenced reads using FastQC revealed one "failure" result for the "Per sequence GC content" test, indicating a possible contamination with DNA from a different species. Subsequently, this scenario was checked with the Kraken2 tool, which classified 21.97% of the original raw reads as one of the taxons present in the reference database, consisting of all humans, bacterial, archaeal, viral, fungal, and protozoan sequences from the NCBI RefSeq database as of May 2022. Of the remaining, unclassified reads, which amounted to 164.9 million reads, 161.6 million reads (98%) mapped to the bovine reference genome ARS-UCD1.2.
IGV (integrative genome viewer) images for ADGRA2, BMPRIB, CDH1, CCDC8, F13A1, MYBPHL, and PHC2 genes are presented in Figures 4 and 5. Most of the variants, such as non-coding transcript variants, intergenic variants, synonymous variants, are considered as unlikely to change protein behavior or with low or no impact on it. Conversely, as can be seen in Figure 3, Ensembl VEP identified a low number of variants that are assumed to have a high impact in the protein (transcript ablation, splice acceptor variant, stop gained, or stop lost).
An initial quality control checks on raw sequenced reads using FastQC revealed one "failure" result for the "Per sequence GC content" test, indicating a possible contamination with DNA from a different species. Subsequently, this scenario was checked with the Kraken2 tool, which classified 21.97% of the original raw reads as one of the taxons present in the reference database, consisting of all humans, bacterial, archaeal, viral, fungal, and protozoan sequences from the NCBI RefSeq database as of May 2022. Of the remaining, unclassified reads, which amounted to 164.9 million reads, 161.6 million reads (98%) mapped to the bovine reference genome ARS-UCD1.2.
IGV (integrative genome viewer) images for ADGRA2, BMPR1B, CCDC8, F13A1, GSTP1 and IRF6 genes are presented in Figures 4 and 5.   [30]. The positon in reads is indicated with red arrows, which contains the mismatches with the reference. In CCDC8 gene, the figure depicts an insertion of 24 nucleotides.   [30]. The positon in reads is indicated with red arrows, which contains the mismatches with the reference. In CCDC8 gene, the figure depicts an insertion of 24 nucleotides.
Animals 2023, 13, x FOR PEER REVIEW 8 of 18 Figure 4. The IGV images of new CCDC8 (a) and BMPR1B (b) high impact target gene variations and with their gene structure consequence that were not reported in literature with implication in humans or animals on CLJ or hypospadias, except for one study [30]. The positon in reads is indicated with red arrows, which contains the mismatches with the reference. In CCDC8 gene, the figure depicts an insertion of 24 nucleotides.   The list of 559 unique bovine genes that were found associated with hypospadias and CLJ in cattle and other mammals, according to literature and the CleftGeneDB database, are presented in an Excel sheet (Supplementary S1 file). The chromosomal positions are those corresponding to the ARS-UCD1.2 Bos taurus reference genome. Supplementary S1 file also contains a list of 379 genes associated with high-impact variants and the list of 4049 genes associated with moderate-impact variants. Table 1 presents the variant annotation genes known to be involved in CLJ and hypospadias in mammals with a high impact on the following overlapped genes: ACVR1, ADGRA2, BHMT2, BMPR1B, CCDC8, CDH1, EGF, F13A1, GSTP1, IRF6, MMP14, MYBPHL, and PHC2. The list of SNPs reported by Ensembl VEP as having a high impact on any of the 559 target genes are presented in an Excel sheet (Supplementary S2 file) and the list with those having a moderate impact on any of the 559 target genes are presented in the Excel sheet named Supplementary S3 file.
Gene annotations based on GO, KEGG, MeSH, and Reactome, which were found significant after over-representation analyses based on genes with moderate and high impact SNPs (379 and 4049 genes, respectively), are included in Supplementary S4 file, with the corresponding sheet names following: GO_high_moderate = Significant Gene Ontology terms for genes overlapped by variations with a high or moderate effect; KEGG_high_moderate = Significant KEGG terms for genes overlapped by variations with a high or moderate effect; MeSH_high = Significant MeSH terms for genes overlapped by variations with a high effect; MeSH_high_moderate = Significant MeSH terms for genes overlapped by variations with a high or moderate effect; Reactome_high_moderate = Significant Reactome terms for genes overlapped by variations with a high or moderate effect.
The biological terms and pathways in the gene annotation databases consist of 5 GO terms, 2 KEGG pathways, 12 MeSH terms, and 5 Reactome pathways. A total of 490 genes were assigned to the Biological Process (BP) GO term category (3 terms
During CT analysis of this case, high scores were observed for orofacial defects according to the cleft lip and jaw scoring system developed by Reinartz et al. [10] in cattle. The right campylognathia seen in this case may be a consequence of the absence of the pressing surface between the mandible and the missing segments of the os incisivum contralateral to the cheilognathoschisis.
Following the analysis of the genomic sequences, the genes with high, high/moderate, and moderate impact on the two congenital diseases were identified, as well as variations of these genes with high or moderate impact. We identified a large number of genes, 379 genes with 1103 genetic variations that can have a major impact on the phenotype. The genes overlapped by variations with high impact or moderate impact on all genes were in the number of 4177. and 221 were with impact on genes of interest, with 16,026 variants with high or moderate impact on all genes, as well as 4049 genes overlapped by variations with moderate impact, having 14,924 variations with moderate impact on the genes on which they were overlapped.
Among them, whole genome sequencing and variant annotation focused on genes known to be involved in CLJ and hypospadias in mammals identified 13 mutations with a high impact on the following overlapped genes: EGF, IRF6, BMPR1B, MMP14, ADGRA2, F13A1, CCDC8, ACVR1, BHMT2, CDH1, GSTP1, PHC2, and MYBPHL where the EGF, IRF6, F13A1, ADGRA2, and GSTP1 genes were found with homozygous mutation, and the CCDC8 and BMPR1B genes were not reported in the literature (except [30]) with implication in humans or animals CLJ or hypospadias, nor did the two genes in the Cleft Gene DataBase.
For the EGF gene (ENSBTAG00000048812), the identified mutation located in the position 6:g.15450370G>A (p.Arg105Ter) revealed a homozygous single nucleotide variant, and the consequence was a stop gained (premature stop codon). The protein encoded by EGF gene has 232 amino acids and the affected amino acid is located at position 105. EGF (epidermal growth factor), a prototypical ligand, and epidermal growth factor receptor (EGFR) are part of the EGF family, being important for cell motility, adhesion, proliferation by enhancing cell cycle progression through the G1 phase, differentiation, and survival being ubiquitously expressed in somatic tissues and in the gonads [75]. EGFR has an important role in normal craniofacial development, EGFR-/with decreased matrix metalloproteinases (MMPs) in mice embryos induced facial anomalies (hypognathia, lack of eyelids, cleft palate) [76]. MMPs, part of a protease family that comprehends adamalysins (ADAMs), disintegrin, and metalloproteinase with thrombospondin motifs (ADAMTSs), with an important role in extracellular matrix remodeling, morphogenesis, have a possible implication in the etiology of CL/P [77]. A correlation between the low expression level of EGF, EGFR, and HOXA13 in the foreskin tissue and hypospadias cases indicates a possible cause of hypospadias [78].
Another gene identified as having a potential high risk was the BMPR1B gene (ENSB-TAG00000002081) (Supplementary S2 file), a member of BMP cytokins. The heterozygous mutation identified in the BMPR1B gene was 6:g.29483476G>A (pGln10X) and the consequence was the introduction of a premature stop codon. Studies indicate that BMP signaling molecules are involved through active ACVR1 (activin A receptor type 1) in the palate epithelium and an increase of this signaling pathway causes a submucous palate [79]. BMPs, cytokines that control the function of many types of cells, exert their role through BMPR1A, BMPR1B, and BMPR2 class receptors with many other receptors in each class, and loss of function can be associated with muscular-skeletal and, cardiovascular diseases, cancer [80], or CP [81]. BMP and growth factors, including EGF, play the role of paracrine mediators for epithelial-mesenchymal interactions during the embryonic development of various organs, including the male reproductive system, EGF stimulating androgen binding activity in the male genital tract, very important hormones for the fetal development of the reproductive system during the hormone-dependent period [40,42,82]. Mutation of BMP4 and BMP7-two factors that play a role in normal urethra development, along with HOXA4 and HOXB6, were identified in Chinese patients with hypospadias, but without a clear association [83].
For the IRF6 gene (ENSBTAG00000002849), the homozygous state mutation identified was 16:g73512380-73512381insT, a frameshift variant consisting of an insertion of a nucleotide at position 162. The IRF6 gene (Interferon Regulatory Factor 6) plays an important role in normal early development. Its function is active in the cells from which the tissues of the face, head, skin, and genitals are formed. In the normal development of the epidermis, IRF6 regulates the balance between keratinocyte proliferation and differentiation, as well as in craniofacial morphogenesis, oral epithelium and breast epithelium [84,85].
Different mutations in the IRF6 locus may increase the risk of syndromic and nonsyndromic forms of orofacial clefts; in humans, for example, mutations in the IRF6 gene are associated with Van der Wounde Syndrome (VWS) (OMIM # 119300) and Polipteal Pterygium syndrome (OMIM # 119500), rare orofacial cleft syndromes associated with lip pits, and cutaneous and limb defects [86,87]. In approximately 70% of VWS cases, mutations in the IRF6 gene are responsible for the disease and Zang et al. (2020) identified a frameshift mutation in IRF6 as a possible cause (NM.006147.3, c.1088-1091delTCTA; p.Ile363ArgfsTer33) [88]. Moreover, the IRF6 gene variations are believed to affect the function of the IRF6 protein in its role as a transcription factor, which may interfere with the normal development of the face. Keratinocyte deficiency in IRF6 affects intracellular adhesion and cell colony morphology, an aspect observed clinically postoperatively in patients with cleft-lip/cleft palate by delaying wound healing [85]. An IRF6 SNP may play a major role in the pathogenesis and risk of developing non-syndromic cleft lip +/− palate in a South Indian population [89]. Loss of IRF6 in mice leads to skin, craniofacial, and limb defects. A minimum level of IRF6 in the periderm and basal epithelial cells is necessary for orofacial development and a partial genetic rescue of the IRF6 knockout embryo with the KRT14 promoter can be performed [90].
The ADGRA2 gene (adhesion G protein-coupled receptor A2)(ENSBTAG00000008814) was found as having high-risk mutation at position 27:33105944C>T (p.Gln1341Thr). The protein consists of 1343 aa, and in our case, the aa at position 1341 was affected by the gene mutation. We found no data in the literature regarding its implication in CLJ in cattle and only little information in other species [91]. Regarding the F13A1 (coagulation factor XIII A chain) gene, the data from the literature does not evidence a linkage with CL/P [92], even if it is mentioned in Cleft GeneDB [60]. The GSTP1 (glutathione S-transferase P1) was involved in the detoxification of electrophilic compounds by glutathione conjugation, I105V polymorphism in mothers and/or children with/without maternal smoking may increase the risk of CL/P [93].
Two heterozygous insertions were identified for the CCDC8 gene (ENSBTAG00000018181), the first being 18:g.53693827_53693828insCAGACAA, an insertion of seven nucleotides classified as a frameshift variant, identified at protein level p.Leu374ProfsTer5 (the leucine is the first amino acid changed in position 374 with proline, the length of the shift frame is 5) and the second 18:g53693828-53693829 ins CAGAGGGCAGAGGCCCC, also classified as frameshift variant. In general, frameshift variants have a high impact on the protein, causing in most cases the production of an altered protein 3D structure. The biochemical and cellular functions of this gene are not clearly known, but the protein encoded by the gene acts as a cofactor required for p-53-mediated apoptosis after DNA damage and may also play a role in growth through interactions with the cytoskeletal adaptor protein obscurin-like 1. Deletion of CCDC8 caused perinatal lethality, intrauterine growth restriction, and placental defects [94]. In humans, the mutation in the 19q13.32 region of CCDC8 is associated with 3M syndrome dwarfism, a rare inherited disorder characterized by low body weight and short length at birth, short stature as adults, widespread skeletal abnormalities, and unusual facial features [94,95]. The N-terminal and C-terminal of the CCDC8 proteins are rather conservative across species, suggesting their important functions.
In the CDH1 gene (ENSBTAG00000015991), the mutation located in position 18:36101086-36101087insT had as a consequence a frameshift, splice region, and intron variant. The CDH1 gene, which encodes epithelial cadherin, is also very important for the cell-cell adhesions, mobility, and proliferation of epithelial cells, being highly expressed in the frontonasal prominence and in the lateral and medial nasal prominences of human embryos, during critical stages of lip and palate development. Pathogenic variants of CDH1 gene have been implicated in syndromic and non-syndromic cleft lip-palate, as well as in hereditary diffuse gastric cancer [96,97].
Analyzing cleft lip and jaw independently by comparing the association of 13 major risk genes in this case with cleft lip and jaw reported in the literature and in the Cleft-GeneDB database, we can say that ACVR1, ADGRA2, BHMT2, BMPR1B, CDH1, EGF, F13A1, GSTP1, IRF6, MMP14, MYBPHL, PHC2, except CCDC8 can be involved in CLJ. As for the hypospadias from the 13 genes, only CCDC8, EGF, IRF6, CDH1, and GSTP1 were reported in literature and in the Gene-Disease Associations dataset [78] as being associated with this anomaly. The protein function of EGF, IRF6, or CDH1 in various tissues is to stimulate the growth of epidermal and epithelial tissues or to sustain cell-cell adhesions, so mutation in them can induce congenital abnormalities in the skin in different anatomical regions, as in this case.
Genomic data suggest the involvement of several mutated genes in the critical period of embryonic development of the affected regions. While these non-synonymous homozygous ADGRA2, EGF, F13A1, IRF6, and GSTP1 genes mutation might encode mutant proteins, they do not provide a prediction of the phenotype, but may contribute to more accurate assessments of the orofacial clefts or hypospadias. It is difficult to predict the impact on the protein, because there are many characteristics of proteins that can be modified, such as folding, protein interactions, dynamics, post-translational modifications, or solubility [98]. Knowing the origin of the parents, having no similar cases of products obtained with the same bull, we consider that most likely these mutations are not hereditary, but occurred spontaneously in the embryonic cells of the affected calf. Moreover, other genes, such as SRD5A2, AR, DGKK, NR5A1, and WNT5A, involved in the steroidogenesis of androgens and in the functionality of hormones, were not identified as having a major risk in this case, but we do not exclude the involvement of endocrine disruptors in the hormonedependent period of the development of the genital segments external as a possible cause in hypospadias. New research data indicate that DNA methylation sites in the penile foreskin are involved in hypospadias etiology and a potential role of environmental factors that can act through epigenetic alterations and even epigenetic inheritance is suggested [99,100].
Our study expands the gene spectrum toward new genes involved in cleft lip and jaw or in hypospadias, but it is possible that there are other molecular processes coordinated by other genes involved in the presented case that have not been identified. Despite these limitations, the present results provide a basis for further investigation of key genes in cleft lip and jaw or in hypospadias in animals.

Conclusions
In summary, in our study, we analyzed a unique case of congenital defects in a Bos taurus crossbred Piedmontese × Wagyu individual using different paraclinical examinations. The genomic analyses based on whole genome sequencing identified 13 mutations with a predicted high impact on the following overlapped genes: ACVR1, ADGRA2, BHMT2, BMPR1B, CCDC8, CDH1, EGF, F13A1, GSTP1, IRF6, MMP14, MYBPHL, and PHC2, of which ADGRA2, EGF, F13A1, GSTP1, and IRF6 showed homozygous mutations, and BMPR1B and CCDC8 had a low or no report in other clinical cases in bovine species regarding CLJ and hypospadias. In veterinary medicine, the genomic data for the two malformations are very rare, although there are numerous clinical reports in the literature about the existence of these congenital pathologies. Therefore, although we cannot specify exactly which of the high-risk genes identified in the study are responsible for the changes identified or if there were environmental factors involved in the appearance of these malformations, we believe that the present study expands the gene spectrum towards new genes involved in cheilognathoschisis or in hypospadias and that the obtained results increase the genetic contribution towards the knowledge of congenital defects in cattle. However, more studies are needed to better understand the specific effects of these genes variants on cheilognathoschisis or on hypospadias in bovine.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ani13101709/s1, Figure S1: Average genome coverage: (a) Coverage histogram, (b) Genome fraction, S1 file: The list of 559 unique bovine genes that were found associated with hypospadias and CLJ in cattle and other mammals, according to literature and the CleftGeneDB database (XLSX); S2 file: The list of SNPs reported by Ensembl VEP as having a high impact on any of the 559 target genes (XLSX); S3 file: The list of SNPs reported by Ensembl VEP as having a moderate impact on any of the 559 target genes (XLSX); S4 file: Gene annotations based on GeneOntology (GO), KEGG, Medical Subject Headings (MeSH), and Reactome, which were found significant after over-representation analyses based on genes with moderate and high impact SNPs (XLSX). Funding: This research paper is supported by the project "Increasing the impact of excellence re-search on the capacity for innovation and technology transfer within USAMVB Timis , oara" code 6PFE, submitted in the competition Program 1-Development of the national system of