G C A T T A C G G C A T genes Article Deviations from Mendelian Inheritance on Bovine X-Chromosome Revealing Recombination, Sex-of-Offspring Effects and Fertility-Related Candidate Genes Samir Id-Lahoucine 1 , Joaquim Casellas 2 , Pablo A. S. Fonseca 1 , Aroa Suárez-Vega 1 , Flavio S. Schenkel 1 and Angela Cánovas 1, * 1 Centre for Genetic Improvement of Livestock, Department of Animal Biosciences, University of Guelph, Guelph, ON N1G 2W1, Canada 2 Departament de Ciència Animal i dels Aliments, Universitat Autònoma de Barcelona, Bellaterra, 08193 Barcelona, Spain * Correspondence: acanovas@uoguelph.ca Abstract: Transmission ratio distortion (TRD), or significant deviations from Mendelian inheritance, is a well-studied phenomenon on autosomal chromosomes, but has not yet received attention on sex chromosomes. TRD was analyzed on 3832 heterosomal single nucleotide polymorphisms (SNPs) and 400 pseudoautosomal SNPs spanning the length of the X-chromosome using 436,651 genotyped Holstein cattle. On the pseudoautosomal region, an opposite sire-TRD pattern between male and female offspring was identified for 149 SNPs. This finding revealed unique SNPs linked to a specific- sex (Y- or X-) chromosome and describes the accumulation of recombination events across the pseudoautosomal region. On the heterosomal region, 13 SNPs and 69 haplotype windows were identified with dam-TRD. Functional analyses for TRD regions highlighted relevant biological functions responsible to regulate spermatogenesis, development of Sertoli cells, homeostasis of Citation: Id-Lahoucine, S.; Casellas, J.; Fonseca, P.A.S.; Suárez-Vega, A.; endometrium tissue and embryonic development. This study uncovered the prevalence of different Schenkel, F.S.; Cánovas, A. Deviations TRD patterns across both heterosomal and pseudoautosomal regions of the X-chromosome and from Mendelian Inheritance on Bovine revealed functional candidate genes for bovine reproduction. X-Chromosome Revealing Recombination, Sex-of-Offspring Effects Keywords: transmission ratio distortion; chromosome X; recombination; bovine and Fertility-Related Candidate Genes. Genes 2022, 13, 2322. https:// doi.org/10.3390/genes13122322 Academic Editor: Scott Pratt 1. Introduction Significant deviations from Mendelian inheritance, well known as transmission ratio Received: 26 October 2022 distortion (TRD), have been observed in a number of diploid organisms [1–5]. It is widely Accepted: 6 December 2022 recognized that the TRD phenomenon is linked to a broad range of biological mechanisms Published: 9 December 2022 associated with basic processes of life that underlie meiosis, gametogenesis, sperm and Publisher’s Note: MDPI stays neutral ova fertility and viability at developmental stages of the reproductive cycle (e.g., embryos, with regard to jurisdictional claims in postnatal metabolism and growth, etc.) [6,7]. The importance of the above-mentioned published maps and institutional affil- events coupled with the current widespread use of DNA technologies, such as genotyping iations. of breeding stock on a routine basis in some livestock industries (e.g., beef and dairy cattle and swine), emphasizes the relevance and potential usefulness of assessing TRD and opening new opportunities for research and extension initiatives. Although TRD remains an unclear and ambiguous phenomenon [8,9], it can be viewed as the ultimate consequence Copyright: © 2022 by the authors. of several genetic factors arising at different stages of the reproductive process and early Licensee MDPI, Basel, Switzerland. neonatal life. Indeed, some cases of TRD, such as the absence of homozygosity, are already This article is an open access article distributed under the terms and used to target recessive phenotypes that affect reproduction, for example, decreased fertility conditions of the Creative Commons and embryo loss in different livestock species [5,10–13]. Attribution (CC BY) license (https:// Despite the recent interest in targeting TRD regions across the whole genome or, alter- creativecommons.org/licenses/by/ natively, evaluating for the absence of homozygous haplotypes in livestock 4.0/). species [4,5,10,11,13–15], TRD on sex chromosomes has not yet received attention. In Genes 2022, 13, 2322. https://doi.org/10.3390/genes13122322 https://www.mdpi.com/journal/genes Genes 2022, 13, 2322 2 of 23 fact, sex chromosomes have been discarded in a wide range of genetic studies of livestock species due to their large heterosomal (i.e., non-pseudoautosomal) part, which requires adaptation of standard analytical approaches. Nevertheless, two studies in humans have found a specific sex-offspring TRD on the X chromosome [16,17]. In this study, we aimed to investigate the inheritance of single nucleotide polymor- phisms (SNPs) and haplotypes from parents to offspring in order to characterize and understand the biology of TRD for the bovine sex chromosomes. The main objectives of this research were: (1) to generalize TRD models with male- and female-offspring-specific TRD components, (2) to determine the pseudoautosomal and heterosomal parts of the sex chromosomes using trios (sire-dam-offspring genotypes), (3) to characterize TRD on both heterosomal and pseudoautosomal parts of the Holstein cattle sex chromosomes, and finally, (4) to annotate genes and to evaluate the functional consequences of TRD regions via functional analyses. 2. Materials and Methods 2.1. Genotype Data Data used in this study comprised a total of 436,651 Holstein dairy cattle genotypes provided by the Canadian Dairy Network, a member of Lactanet (Guelph, ON, Canada). These animals were sampled from all available Holstein genotypes (>1 million genotypes; October 2017) and were genotyped within the first 90 days after birth to minimize offspring selection issues [13,18]. Animals were genotyped by various SNP genotyping arrays ranging from 2900 to 777,962 SNPs (Supplementary Table S10) and mapped to the UMD3.1 Bos taurus genome assembly. All SNP markers mapped in sex chromosomes and genotyped for at least 10 sire-dam-offspring trios were selected for analyses, resulting in a total of 6577 SNPs distributed across the X-chromosome (Supplementary Table S1). The number of trios available for SNPs ranged from 10 to 340,332 with an average of 18,980 trios. The maximum number of genotyped sires and dams was 5917 and 130,501, respectively. 2.2. Imputation In order to carry out haplotype analyses, 974 SNPs were selected from the heterosomal part of the X-chromosome (Supplementary Table S2). The number of raw genotypes ranged from 23,704 to 402,267 with an average of 141,821 genotypes across 974 SNPs. The length of the X-chromosome covered by these markers was between 29,245 and 137,012,896 bp. Data was phased and imputed using a Fortran program using familial information (i.e., parent- offspring trios and sibling pairs) with variable sliding windows ranging from 2 to 200 SNPs. 2.3. Male- and Female-Offspring Specific-TRD Taking an autosomal (or pseudo-autosomal) chromosome region as a starting point, the analysis of TRD for male (♂ ) and female offspring (♀ ) can be carried out by modeling the probability of transmission of alleles from heterozygous parents [4,14]. Considering male offspring as an example, the probability of transmission of the A allele from heterozygous parents (A/B) can be stated as: P(A) = 1 − P(B) = 0.5 + α♂ and P(B) = 1 − P(A) = 0.5 − α♂ , where: α♂ is the TRD parameter ranging between −0.5 and 0.5. As described by Casellas et al. (2017) [4], this model can be easily expanded to account for sire- (αs♂ ) and dam-specific (αd♂ ) TRD. The same applies to female offspring to obtain α♀ , αs♀ and αd♀ . The previous models can be generalized to the analysis of TRD phenomena in the non-autosomal (i.e., heterosomal) region of the X chromosome, but restricted to dam- specific TRD. 2.4. Statistical Analyses Two different approaches were used for analyzing TRD across the bovine X-chromosome. Transmission ratio distortion was evaluated using both SNP-by-SNP (a total of 6577 SNPs) Genes 2022, 13, 2322 3 of 23 and a sliding windows approach of 2-, 4-, 10- and 20-SNP across 974 imputed SNPs as described by Id-Lahoucine et al. (2019) [13]. For sliding windows, the biallelic-haplotype procedure was applied. All the analyses were performed within a Bayesian framework using an adapted version of TRDscan software [13] with a unique Monte Carlo Markov chain of 110,000 iterations, where the first 10,000 iterations were discarded as burn-in. The statistical significance of the TRD was tested using a BF [19]. 2.5. Pseudoautosomal Characterization Mendelian inconsistency was estimated under a pseudoautosomal and heterosomal inheritance models using trios of genotypes (from dam, sire and offspring) in order to deter- mine the nature of SNPs across the X-chromosome. Both analyses of LD and recombination events were performed. For recombination events, the sum of recombination probabilities was calculated for markers in pseudoautosomal parts of the X-chromosome using the package hsphase in R [20]. In addition, for better characterization of the pseudoautosomal regions on X- and Y- chromosomes, the regions with pseudoautosomal inheritance models were initially selected to perform a “nucleotide BLAST” analysis on Basic Local Alignment Search Tool (BLAST; https://blast.ncbi.nlm.nih.gov/Blast.cgi accessed on 10 December 2020) against the bovine reference genome. The sequence was retrieved using the bovine UMD3.1 genome assembly from NCBI database. After this step, a short sequence of 20 bp around the genomic position of the markers was extracted to obtain the nucleotide context around the TRD markers. Using the genome data viewer, regions with homology between the short sequence around the markers and the Y-chromosome were investigated. The whole assembly of the bovine Y-chromosome is completely available for visualization only for the bovine genome assembly Btau_4.6.1 and, therefore, this assembly was used for the data visualization of the Y-chromosome. It is important to reinforce that neither for UMD 3.1 or ARS-UCD 1.2 an assemble of bovine Y-chromosome is available, which results in the impossibility to use these assemblies for the abovementioned analysis. 2.6. Functional Analyses The genes mapped around the TRD regions identified by both SNP- and haplotype- based methods were annotated using R software (Version 3.2.3) and the R package: Genomic functional Annotation in Livestock for positional candidate Loci (GALLO) [21]. The genes mapped around the TRD regions identified using the SNP approach were annotated using an interval of 100 Kilobase (Kb) upstream and downstream from the SNP position. For the haplotype approach, the genes were annotated using the start and end coordinates from each haplotype based on the UMD 3.1 assembly from the bovine reference genome. Subsequently, in addition to the bovine Ensembl ID and official associated Gene Symbol, the respective orthologous genes in human and mouse were also annotated for each positional candidate (Ensembl ID and Gene Symbol) using the R package biomaRt [22]. In this step, only the human ortholog showing a similarity higher than 75% with the bovine annotated gene were retained for subsequent analyses. The highest |α| and the respective log10 (BF) of the associated marker (SNPs or haplotypes with TRD) were assigned for each positional candidate gene. Functional analysis were performed using the IPA (Ingenuity Pathways Analysis) software (QIAGEN Redwood City, www.qiagen.com/ingenuity 15 December 2020) to identify the canonical metabolic pathways, diseases and functions associated with the list of candidate genes harboring deviations from Mendelian inheritance on the bovine X-chromosome (including both SNPs and haplotypes). In addition, functional analysis was performed to explore the existence of signaling networks connecting those genes [23]. The significantly enriched metabolic pathways were identified using a significance threshold of false-discovery rate (FDR) < 0.05. Genes 2022, 13, 2322 4 of 23 3. Results and Discussion 3.1. Characterization of X-Chromosome Using Parent-Offspring Genotyped Trios The mammalian sex chromosomes are characterized by two main types of regions: a short region of sequence homology between the X- and the Y-chromosomes, known as the pseudoautosomal region and a longer heterosomal (sex-linked) region specific for each sex chromosome. The pseudoautosomal region in cattle is located in the long arm of the X-chromosome and the short arm of the Y-chromosome. Its physical length is estimated to be between 5–9 Mb [24]. Other studies specified a size of 8 Mbp [25] and 5.7 Mb on the ARS-UCD1.2 reference genome assembly [26]. Within the pseudoautosomal region, pairing and recombination occurs between X- and Y-chromosomes in males as it does for autosomal chromosomes. On the other hand, the heterosomal regions do not pair and recombine during meiosis in males, due to the reduced homology between the chromosomes. In this study, the patterns of Mendelian incompatibilities, estimated under pseudoautosomal and heterosomal inheritance models for each individual SNP together with the patterns of their adjacent SNPs, were used to confirm their nature (i.e., pseudoautosomal or heterosomal). The short arm of the X-chromosome and part of the long arm showed high incompatibility under the pseudoautosomal inheritance model whereas null Mendelian inconsistencies under the heterosomal inheritance model (29,245–137,034,696 bp; Figure 1). Alternatively, the extreme part of the long arm of the X-chromosome (143,865,210–148,816,634 bp) dis- played high and null Mendelian inconsistencies under heterosomal and pseudoautosomal inheritance models, respectively. This opposite behavior supported its pseudoautosomal nature described in the bovine genome [24,25,27]. However, unique individual SNPs showed Mendelian incompatibilities discrepant with their adjacent SNPs, suggesting possi- ble genotyping artifacts or inaccurate annotation. These discrepant individual SNPs were Genes 2022, 13, x FOR PEER REVIEW 5 of 22 removed from further TRD analyses. On the other hand, small regions of adjacent SNPs between these two clear extremes (heterosomal and pseudoautosomal parts) fitted either a heterosomal or a pseudoautosomal inheritance model. 1. Mendelian Figure 1. Figure Mendelian inconsistencies across bovine inconsistencies across bovine X-chromosome based on X-chromosome based on the the pseudoautosomal pseudoautosomal the heterosomal inheritance model (positive value and black color) and the heterosomal inheritance model (negative values and grey color). From the From thepreliminary preliminaryscan, scan,5935 5935andand780780 SNPs SNPs were were shown shown as exhibiting as exhibiting heteroso- heterosomal mal pseudoautosomal and and pseudoautosomal inheritance, inheritance, respectively. respectively. In order In order to guarantee to guarantee a reasonable a reasonable statis- statistical power, only SNPs with at least 100 trios (sire-dam-offspring) of genotypes tical power, only SNPs with at least 100 trios (sire-dam-offspring) of genotypes and and aa Mendelian inconsistency < 1% were maintained for TRD analyses (Supplementary Table S1). After this editing, the number of SNPs was reduced to 3832 and 400 with heterosomal and pseudoautosomal nature, respectively. The two main parts of the X-chromosome in- cluded 3680 heterosomal SNPs (BTX:29,245–137,034,696 bp) and 322 pseudoautosomal Genes 2022, 13, 2322 5 of 23 Mendelian inconsistency < 1% were maintained for TRD analyses (Supplementary Table S1). After this editing, the number of SNPs was reduced to 3832 and 400 with heterosomal and pseudoautosomal nature, respectively. The two main parts of the X-chromosome included 3680 heterosomal SNPs (BTX:29,245–137,034,696 bp) and 322 pseudoautosomal SNPs (BTX:143,865,210–148,816,634 bp). The undetermined regions of the X-chromosome (BTX:137,034,696–143,865,210 bp) displayed different inheritance patterns in the follow- ing order: 8 pseudoautosomal SNPs (BTX:137,181,696–137,476,239 bp), 9 heterosomal SNPs (BTX:137,531,272–137,717,498 bp), 54 pseudoautosomal SNPs (BTX:137,951,394– 139,042,034 bp), 30 heterosomal SNPs (BTX:139,126,989–140,104,202 bp), 16 pseudoautoso- mal SNPs (BTX:140,116,340–140,380,961 bp) and 113 heterosomal SNPs (BTX:140,445,475– 143,832,372 bp). The total sum of these parts was 141.6 Mb and 6.6 Mb for pseudoautosomal and heterosomal regions, respectively. The physical length of pseudoautosomal regions is in accordance with Das et al. (2009) [24], i.e., between 5–9 Mb. It must be noted that the pseudoautosomal part of the cattle X-chromosome is demarcated by the pseudoautosomal boundary, a border where the sequence similarity between the X- and the Y-chromosomes decreases (to 80–50%) and regions specific to each sex chromosome begin [28]. Couldrey et al. (2016) [25], reported 2 regions (presumably pseudoautosomal) with male heterozygosity greater than 1%, extending from 137 to 141 Mb and from 144 Mb to the end of the chromosome, providing evidence of the different inheritance patterns observed between 137 and 144 Mb. In addition, McClure et al. (2018) [29] characterized 101 pseudoautosomal SNPs for the X-chromosome from 137,330,600 to 148,482,298 bp. Among these 101 pseudoautosomal SNPs, 69 were in accordance with our results based on Mendelian compatibilities. However, 10 SNPs were assumed heterosomal in nature according to our results whereas were pseudoautosomal by McClure et al. (2018) [29]. Specifically, these 10 SNPs displayed ≥ 9% and ≤0.31% Mendelian inconsistencies under pseudoautosomal and heterosomal inheritance models, respectively, and their heterosomal nature was supported by the similar pattern of the corresponding adjacent SNPs. These SNPs are also reported with null male heterozygosity (≤0.02) and moderate-to-high female heterozygosity (≥0.29) by McClure et al. (2018) [29]. 3.2. Deviations from Mendelian Inheritance on the Pseudoautosomal Part of the Sex Chromosomes Decisive evidence with Bayes factor (BF; ≥100) according to Jeffreys’ scale [30] was identified for TRD across the Holstein sex chromosomes. The pseudoautosomal SNPs on the sex chromosomes revealed an important pattern of sire-TRD with specific sex-offspring components. A total of 180 SNPs was characterized by significantly displaying oppo- site sire-TRD between offspring sexes (Figure 2A; Supplementary Table S4). SNPs that exhibited specific sire-TRD via one single offspring sex (only male or female) were also ob- served, being 26 SNPs with sire-TRD via male and 19 SNPs with sire-TRD via female (Tables 1 and S4). In addition, 3 SNPs displayed an opposite sire-TRD between off- spring sexes and also a significant dam-TRD pattern only via one single offspring sex (Supplementary Table S4). On the other hand, additional 5, 2 and 1 SNPs were also iden- tified with BF ≥ 100 for overall TRD (specific offspring sex), dam-TRD and sire-TRD (unspecific offspring sex), respectively. Among them, only one individual SNP showed a probability of random TRD ≤ 0.001 (BTX:147,104,359; sire-TRD = −0.04 with BF = 1020.18 , 625 heterozygous sires and 1458 under-represented offspring). Genes 2022, 13, x FOR PEER REVIEW 7 Genes 2022, 13, 2322 6 of 23 Figure 2. Schematic representation of the long arm of X-chromosome (Xq) and short arm of Y- chromosome (Yp) and sire-TRD patterns. (A) Magnitude of sire-TRD with opposite sign between Figure male- (blue color) and2.female-offspring Schematic representation (red color)ofonthe long arm of X-chromosome pseudoautosomal (Xq) and short SNPs. (B) Schematic repre-arm of Y- mosome (Yp) and sire-TRD patterns. (A) Magnitude of sire-TRD sentation of X and Y peudoautosomal regions (PARX and PARY) and peudoautosomal boundaires with opposite sign between m (PABX and PABY). (C) Density of SNP available for analyses (pseudoautosomal with red color and repres (blue color) and female-offspring (red color) on pseudoautosomal SNPs. (B) Schematic tion of X and Y peudoautosomal regions (PARX and PARY) and peudoautosomal bound heterosomal with blue color). The pink and blue circles represent the first two markers with α = 0.5 (PABX and PABY). (C) Density of SNP available for analyses (pseudoautosomal with red colo (BTX:143,865,210 and BTX:143,870,595), respectively. The red circle represents the last marker with heterosomal with blue color). The pink and blue circles represent the first two markers with α an opposite sire-TRD effect between (BTX:143,865,210 male and female offspring and BTX:143,870,595), (BTX:148,789,832). respectively. The orange The red circle represents and the last marker green rectanglesancorrespond to the markers opposite sire-TRD with opposite effect between male andsire-TRD effects female between offspring male and female (BTX:148,789,832). The orange offspring and dam-TRD for malecorrespond green rectangles (orange) and female to the (green) markers offspring. with oppositeThe blue diamonds sire-TRD represent effects between male and fe the region withoffspring highest homology and dam-TRD between for the maleX- (orange) and Y-chromosomes (BTX:144,251,086–144,399,458). and female (green) offspring. The blue diamonds r It is important sent the region to highlight that with this ishighest homology a schematic image between and somethe X- where genes and Y-chromosomes suppressed. (BTX:144,251 144,399,458). It is important to highlight that this is a schematic image and some genes where pressed. Genes 2022, 13, 2322 7 of 23 Table 1. Specific sex-offspring for sire-TRD in Holstein X-chromosome. Male-Offspring Female-Offspring Coordinates Frequency n◦ Hetero 1 n◦ Hetero n◦ Info 2 Sire-TRD n◦ Hetero n◦ Hetero n◦ Info Sire-TRD in Base Pairs BF BF Sires Dams Offspring (αs♂ ) Sires Dams Offspring (αs♀ ) 137,464,119 0.63 376 2307 10,297 0.02 0.79 275 1625 7100 −0.04 104.2 138,080,118 0.46 699 4368 19,985 −0.01 0.04 1281 20,752 48,445 0.03 1022.8 144,824,704 0.31 267 1691 8037 0.01 0.03 365 2258 8646 0.04 103.4 144,927,188 0.83 7 88 228 0.36 103.8 8 46 107 −0.14 0.63 145,093,201 0.24 403 2284 10,759 0.00 0.02 361 1879 9007 −0.04 104.2 145,322,565 0.67 280 1532 9045 0.02 1 364 1946 9432 −0.08 1027.4 146,268,179 0.14 216 1351 6317 −0.15 1058.9 334 5575 13,299 −0.02 3.98 146,403,456 0.08 14 51 155 −0.11 0.79 13 19 115 0.29 105.9 147,269,713 0.85 56 162 809 0.18 109.2 33 132 442 0.01 0.08 147,433,504 0.52 302 2111 10,749 −0.07 1018.8 406 2793 11,668 0.02 0.63 147,587,823 0.58 98 290 1456 0.11 108.3 66 209 764 −0.03 0.13 148,398,942 0.22 477 2871 14,803 0.02 10 851 14,732 37,427 0.05 1031.7 148,498,697 0.34 267 1685 8671 −0.04 102.9 338 2173 8922 0.00 0.02 148,665,523 0.61 269 1853 8777 0.05 105.4 351 2335 8873 −0.01 0.10 148,682,952 0.36 98 311 1556 −0.11 108.8 67 225 795 0.03 0.13 αs♂: male-sire-TRD, αs♀: female-sire-TRD, BF: Bayes factor. 1 Number of heterozygous. 2 Number of informative offspring (from heterozygous parents). 3.3. Opposite Sire-TRD between Male- and Female-Offspring Genes 2022, 13, 2322 Opposite sire-TRD between male- and female-offspring were observed across 8 of 23 the whole pseudoautosomal part of sex chromosomes with a specific pattern (Figure 2). This TRD pattern was observed with a preferential transmission of one specific allele to male- 3.3. Opposite offspring and theSire-TRD between opposite Male- allele and Female-Offspring to the female-offspring. Figure 3 shows how strong the sire-TRD Opposite sire-TRD was between between male- and male- and female-offspring female-offspring were observed among different matings across for thethe SNP whole pseudoautosomal part of sex chromosomes with a specific (BTX:144,617,673) with the highest BF. Taking into account the approximate empirical null pattern (Figure 2). This TRD pattern was observed with a preferential transmission of one specific allele to male- distribution of TRD [13] at ≤0.001% margin error for random TRD, the number of regions offspring and the opposite allele to the female-offspring. Figure 3 shows how strong the with this pattern reduced from 180 to 149 SNPs spanning the length of the pseudoautoso- sire-TRD was between male- and female-offspring among different matings for the SNP mal(BTX:144,617,673) part of the sex chromosome. with the highestHere, the important BF. Taking into accountdiscovery resides,empirical the approximate as shown in Fig- null uredistribution 2A, in the strong of TRDsire-TRD observed [13] at ≤0.001% for SNPs margin at the error for beginning random TRD, the of number the pseudoautosomal of regions region and its marginal reduction with the distance until being with this pattern reduced from 180 to 149 SNPs spanning the length of the pseudoauto- null on the extreme of the X-chromosome. somal part of the It issex important chromosome. to highlight Here, the that this pattern important is primarily discovery resides, attributed as shown into the mainFigure 2A, in the strongpart pseudoautosomal sire-TRD of theobserved chromosomesfor SNPs(i.e., at the beginning of the pseudoautoso- 143,865,210–148,816,634 bp), where mal region and its marginal reduction with the distance until the first 8 SNPs (143,865,210–144,085,178) displayed TRD ranging from 0.47 to 0.50 and being null on the extreme withof full the X-chromosome. linkage to one sex It is (i.e., important to highlight that this null recombination). Thispattern is primarily result indicates attributed a strong linkage to the main pseudoautosomal part of the chromosomes (i.e., 143,865,210–148,816,634 bp), of some alleles to X- or Y-chromosomes at the beginning of pseudoautosomal region. Ad- where the first 8 SNPs (143,865,210–144,085,178) displayed TRD ranging from 0.47 to 0.50 dition, and the withdecay of opposite full linkage to one sire-TRD sex (i.e., nullbetween male- andThis recombination). female-offspring result indicates witha strongthe dis- tance describe linkage the accumulated of some alleles to X- or recombination Y-chromosomes events at theacross the pseudoautosomal beginning of pseudoautosomal part of theregion. cattle genome and suggest a consistent level of these events. The Addition, the decay of opposite sire-TRD between male- and female-offspring with results of recombina- tionthe analyses distancesupported describe the these patterns recombination accumulated showing an absence of recombination events across events at the the pseudoautosomal part of the beginning cattle of this genome and suggest pseudoautosomal a consistent genomic regionlevel and of these consistent a more events. Thelevel results of of recom- recombination analyses supported these patterns showing an bination events started as TRD magnitude decreased (Figure 4). Therefore, specific sex- absence of recombination events at offspring TRDthe revealed beginninganother of this pseudoautosomal genomic region way to study recombination and and a moredisequilibrium linkage consistent level of recombination events started as TRD magnitude decreased (Figure 4). Therefore, (LD) in sex chromosomes. Thus, whereas marker alleles in strong LD to specific sex chro- specific sex-offspring TRD revealed another way to study recombination and linkage dise- mosome (X or(LD) quilibrium Y) were in sexobserved chromosomes. at the beginning Thus, whereasofmarker the pseudoautosomal alleles in strong LDregion, to specificthis LD decreases sex chromosome (X or Y) were observed at the beginning of the pseudoautosomal region, be- at longer distances from the beginning of the pseudoautosomal region until ingthis nullLD ondecreases the extreme of the at longer X-chromosome. distances This expected from the beginning partial sex linkage of the pseudoautosomal on pseu- region doautosomal until being nullregions on thewas previously extreme reported in humans of the X-chromosome. [31]. Moreover, This expected partial sexthis result linkage on sug- pseudoautosomal regions was previously reported in humans gested that novel mutations generated at the beginning of the pseudoautosomal part will [31]. Moreover, this result suggested be linked to athat novel sex specific mutations for long generated time than at the beginning potential of the pseudoautosomal mutation occurred on thepart extreme will be linked of the X-chromosome. to a specific sex for long time than potential mutation occurred on the extreme of the X-chromosome. Figure 3. Distribution of male and female offspring for informative mattings for the SNP with highest Figure 3. Distribution BF with of male opposite sire-TRD and female between offspring for informative mattings for the SNP with high- sex offspring. est BF with opposite sire-TRD between sex offspring. Genes 2022, 13, x FOR PEER REVIEW 9 of 22 Genes 2022, 13, 2322 9 of 23 Figure 4. Figure Overlappingbetween 4. Overlapping betweenopposite oppositesire-TRD sire-TRD(right (righty-axis) y-axis)and andsum sumofof thethe probability probability of re- of recombination combination (left (left y-axis) y-axis) in inthethemain main pseudoautosomal pseudoautosomal part partofofthethe bovine X-chromosome bovine X-chromosome (143,865,210–148,816,634 bp). (143,865,210–148,816,634 bp).The The x-axis shows x-axis showsthe the genomic coordinated genomic in Mega coordinated in bases Mega(Mb). basesThe red The (Mb). red andand blueblue dotsdots correspond correspond to theto the sire-TRD sire-TRD magnitude magnitude in the in the male male and andoffspring, female female offspring, respec- respectively. tively. This opposite sire-TRD pattern could reinforce the high recombination rate in the cattle pseudoautosomal This oppositeregion, as itpattern sire-TRD was also described could reinforce on other mammalian the high species,rate recombination which in the cattle pseudoautosomal region, as it was also described on other mammalianonspecies, exceed the genome average by 10–20 times [32]. Notice that the recombination events the sex chromosomes in males are expected to be extraordinarily restricted to a short genomic which exceed the genome average by 10–20 times [32]. Notice that the recombination segment (i.e., pseudoautosomal regions), appearing to be a recombinational hotspot in events on the sex chromosomes in males are expected to be extraordinarily restricted to a male meiosis [31]. In addition, the crossover in the pseudoautosomal region was described short genomic segment (i.e., pseudoautosomal regions), appearing to be a recombina- to be essential for the proper disjunction of X- and Y-chromosomes in male meiosis and a tional deletionhotspot in male meiosis of pseudoautosomal [31]. region was Inreported addition, to the crossover result in the pseudoautosomal in male sterility [33]. However, region was described biologically, to be essential the mechanisms by which for this the proper mandatory disjunction crossoverof X- and Y-chromosomes is achieved remains in male meiosis unknown and a deletion [33]. Moreover, of pseudoautosomal the results of estimated LD region was reported using squared to result correlation in male between sterility [33]. However, biologically, the mechanisms by which this mandatory alleles also supported the high recombination rate in the cattle pseudoautosomal region. A crossover is achieved decay of LDremains unknown in the main [33]. Moreover, pseudoautosomal region theinresults of estimated comparison LD using squared to the heterosomal and correlation betweenwas autosomal regions alleles also supported observed the high (Supplementary recombination Figure rate in the S1). Furthermore, cattle pseu- a relatively small decay ofregion. doautosomal LD on heterosomal A decay of region LD in intherelation to autosomal regionregion main pseudoautosomal was also in observed comparison (Supplementary to the heterosomal Figure andS1). autosomal regions was observed (Supplementary Figure S1). Fur- On the thermore, other hand, a relatively 15 pseudoautosomal small SNPs (from 180) decay of LD on heterosomal with region in opposite relation tosire-TRD autosomal patterns were in the undetermined regions region was also observed (Supplementary Figure S1). (i.e., BTX:137,181,696–137,476,239 bp, BTX: 137,951,394–139,042,034 bp and BTX:140,116,340–140,380,961 bp). In On the other hand, 15 pseudoautosomal SNPs (from 180) with opposite sire-TRD contrast to the main pseudoautosomal part, we did not observe a similar pattern (i.e., the decreased mag- patterns were in the undetermined regions (i.e., BTX:137,181,696–137,476,239 bp, nitude of TRD with distance). It is important to emphasize that opposite sire-TRD be- BTX:137,951,394–139,042,034 bp and BTX:140,116,340–140,380,961 bp). In contrast to the tween male- and female-offspring not necessarily imply recombination between X- and main pseudoautosomal part, we did not observe a similar pattern (i.e., the decreased mag- Y-chromosomes. Therefore, this pattern could also be explained by different linkage of nitude of X- alleles to TRDandwith distance). It is important to emphasize that opposite sire-TRD between Y-chromosome. male- and female-offspring not necessarily imply recombination between X- and Y-chro- mosomes. Therefore, this pattern could also be explained by different linkage of alleles to X- and Y-chromosome. Genes 2022, 13, 2322 10 of 23 3.4. Specific Sex-Offspring Sire-TRD Across the X-chromosome, we identified 45 initial SNPs that are exhibiting sire-TRD via one single sex whereas null TRD in the opposite sex. After discarding possible random TRD at ≤0.001%, according to the approximate empirical null distribution of TRD [13], 8 and 7 SNPs showed male- and female-offspring sire-TRD, respectively (Table 1). As introduced previously, similar patterns of TRD on only one sex (specific sex-offspring TRD) were already identified in the human X-chromosome using microsatellite loci [16]. In fact, related patterns were also described for some disease (e.g., retinoblastoma) and were hypothesized to be associated with defective imprinting [34]. In addition, the sire-TRD pattern could be due to the fact that DNA sequences remain methylated for a much longer period in the male germline than in the female germline [35]. Within this context, an adequate imprinting may be required for viability of male or female embryos and the identified SNPs could be mapped in potential genes with differential imprinting patterns. It is notably important to mention that some of these SNPs may be in the same genes as distance between most of them was small (9 SNPs were at < 500 kb from the adjacent significant SNP). The different magnitudes of TRD could also support imprinting for some of them, as TRD ranged from 0.04 to 0.36 in male-offspring and between 0.03 to 0.29 in female-offspring. For these SNPs, the number of under-represented offspring ranges from 47 to 1014 (1 SNP with ≥1000) for males and from 54 to 2183 (3 SNP with ≥1000) for females. The number of heterozygous sires ranges from 7 to 302 and from 13 to 1281 for male- and female-offspring, respectively. The minor allele frequencies of the detected SNPs ranges from 0.08 to 0.48. It is important to mention that moderate-to-high TRD magnitude signals were observed in regions with low or extremely low frequencies, suggesting the effects of TRD selection as observed in autosomal chromosomes as well [13]. Notably, the moderate frequencies of some SNPs in this case could be explained by the null TRD observed in one parent and/or offspring (male or female). Nevertheless, mechanisms such as imprinting could have the same consequences since the not preferentially transmitted allele is still observable in the offspring generation. 3.5. Deviations from Mendelian Inheritance on the Heterosomal Part of the X-Chromosome 3.5.1. SNP-by-SNP Analysis Among 3832 heterosomal SNPs, 12 SNPs were significantly detected with decisive evidence (BF ≥ 100) for dam-TRD (Table 2). Among them, 7 SNPs were significantly observed via male, 2 via female and 3 via both (i.e., non-sex-offspring TRD). However, these SNPs were reduced to 9, 8 and 2 considering a chance of random TRD of ≤0.1, ≤0.01 and ≤0.001%, respectively. The SNP with the highest BF displayed 3675 under-represented offspring with αd = 0.02. The maximum magnitude observed for dam-TRD was |0.12| and displaying 287 under-represented offspring. In general, the prevalence of TRD across the heterosomal regions was small in comparison to pseudoautosomal and autosomal regions. This observation could be related to the biological differences between sexes, such as the processes of germ cell formation in males versus females (e.g., cell divisions rates in spermatogenesis versus oogenesis) [36]. Genes 2022, 13, 2322 11 of 23 Table 2. Significant dam-TRD, male- and female-dam-TRD among heterosomal SNPs in Holstein X-chromosome. A × AB 2 B × AB Coordinates in Offspring-Sex n◦ Hetero 1 Frequency TRD Parameters Base Pairs Dams AA/A 3 AB B A AB BB/B M/(F) 4,134,094 - 391 111 46 88 209 150 623 0.15 (0.18) αd = −0.12 (BF = 1013.3 ; ≤0.001 4 ) 40,239,653 - 42,334 37,552 23,430 11,719 7816 15,685 22,229 0.64 (0.62) αd = 0.02 (BF = 1020.6 ; ≤0.01) 85,463,982 - 13,963 294 241 98 3401 9983 14,552 0.03 (0.04) αd = −0.02 (BF = 109.2 ; ≤0.01) M 640 854 0 668 58 0 39 0.94 αd♂ = 0.06 (BF = 104.14 ; ≤0.01) 2,230,488 F 477 632 628 0 0 20 18 0.93 αd♀ = 0 (BF = 0.06; >5) M 81 88 0 56 42 0 20 0.52 αd♂ = 0.13 (BF = 107.2; ≤1) 3,873,466 F 48 35 25 0 0 19 32 0.49 αd♀ = −0.01 (BF = 0.12; >5) M 8389 10,984 0 10,371 1785 0 1673 0.86 αd♂ = 0.01 (BF = 234.4; ≤0.01) 5,451,082 F 31,944 25,802 25,820 0 0 3570 3826 0.85 αd♀ = 0 (BF = 0.01; >5) M 1364 846 0 710 931 0 795 0.46 αd♂ = 0.04 (BF = 103.2 ; ≤1) 22,030,046 F 1067 603 591 0 0 713 672 0.43 αd♀ = 0.01 (BF = 0.04; >5) M 12,977 5812 0 6482 13,607 0 14,685 0.34 αd♂ = −0.02 (BF = 1014.3 ; ≤0.001) 65,136,071 F 49,960 13,780 14,055 0 0 32,898 33,032 0.37 αd♀ = 0 (BF = 0.01; >5) M 1627 425 0 513 1658 0 1943 0.21 αd♂ = −0.04 (BF = 104.93 ; ≤1) 106,467,781 F 1128 305 309 0 0 1234 1265 0.21 αd♀ = −0.01 (BF = 0.03; >5) M 13,163 7425 0 7888 12,371 0 13,080 0.38 αd♂ = −0.01 (BF = 104.78 ; ≤0.01) 143,832,372 F 52,914 19,253 18,869 0 0 30,256 30,222 0.44 αd♀ = 0.0 (BF = 0.01; >5) M 785 1192 0 1146 53 0 41 0.88 αd♂ = 0.01 (BF = 0.05; >5) 67,099,097 F 504 936 586 0 0 28 37 0.91 αd♀ = 0.11 (BF = 1014.5 ; ≤0.01) M 767 1116 0 1082 56 0 41 0.95 αd♂ = 0.01 (BF = 0.04; >5) 125,033,467 F 613 880 708 0 0 45 40 0.95 αd♀ = 0.05 (BF = 363.1; ≤0.1) M: male, F: female, BF: Bayes factor, αd : dam-TRD, αd♂ : male-dam-TRD, αd♀ : female-dam-TRD. 1 Number of heterozygous. 2 Sire × dam mating genotypes. 3 Offspring genotype. 4 Probability of random TRD. Genes 2022, 13, 2322 12 of 23 3.5.2. Haplotype Analysis Across 974 SNPs of the heterosomal part of the X-chromosome, haplotypes displaying TRD were detected with decisive evidence of BF (≥100). After discarding regions with <10 heterozygous dams and <50 informative offspring, the number of regions was 944, being 44, 71, 282 and 597 detected for 2-, 4-, 10- and 20-SNP haplotypes, respectively. The number of regions at <0.001% margin of error reduced to 14, 26, 138 and 319, respectively. When correcting for multiple tests with kernel smoothing, the number of regions reduced to 69 haplotype windows explaining the observed TRD among the corresponding linked regions. The TRD was observed in both male- and female-offspring and with strong magnitudes ranging up to |0.48|. The top 20 regions according to BF are summarized in Table 3 and the complete list in the Supplementary Table S5. It must be emphasized that, even though 69 haplotypes were characterized, the maximum number of under-represented offspring observed was 193, showing a reduced relevance in comparison to autosomal chromosomes. Table 3. Significant dam-TRD among heterosomal part in Holstein X-chromosome for haplotype analyses. n◦ Hetero 2 A × AB 3 B × AB Frequency Coordinates in Base Pairs 1 Dam-TRD BF Dams AA/A 4 AB B A AB BB/B % 134,602,363:135,936,247 4 10 −0.49 1019.80 0 0 0 0 0 72 0.004 134,841,753:136,063,445 11 −0.49 1019.80 0 0 0 0 0 72 0.004 134,624,675:135,970,186 10 −0.49 1019.51 0 0 0 0 0 71 0.004 119,781,376:121,600,055 11 −0.44 1019.13 0 0 0 1 4 91 0.007 123,894,981:125,905,926 101 −0.22 1018.12 0 0 0 46 72 310 0.109 053,138,058:058,932,871 18 −0.38 1016.11 0 0 0 3 10 103 0.021 056,821,723:059,606,508 12 −0.47 1013.65 0 0 0 0 1 56 0.005 056,846,129:059,878,165 12 −0.47 1013.65 0 0 0 0 1 56 0.005 129,911,571:132,268,728 15 −0.44 1013.11 0 0 0 0 3 62 0.007 130,070,306:132,334,180 15 −0.44 1013.11 0 0 0 0 3 62 0.007 130,112,536:132,620,448 15 −0.44 1013.11 0 0 0 0 3 62 0.007 055,893,493:059,564,066 44 −0.24 1012.55 0 0 0 31 40 199 0.058 129,070,088:131,766,182 12 −0.45 1011.64 0 0 0 1 1 53 0.006 129,145,800:131,806,987 12 −0.45 1011.64 0 0 0 1 1 53 0.006 129,193,110:131,960,783 12 −0.45 1011.64 0 0 0 1 1 53 0.006 129,718,409:132,109,716 12 −0.45 1011.64 0 0 0 1 1 53 0.006 129,887,133:132,211,189 12 −0.45 1011.64 0 0 0 1 1 53 0.006 023,761,259:026,607,222 30 −0.29 1011.07 0 0 0 5 25 118 0.025 056,962,509:060,125,574 13 −0.40 1010.59 0 0 0 2 4 62 0.008 100,533,235:105,003,475 14 −0.38 1010.47 0 0 0 5 3 67 0.011 BF: Bayes factor, αd : dam-TRD. 1 Number of SNPs on the haplotype window was 20. 2 Number of heterozygous dams. 3 Sire × dam mating genotypes. 4 Offspring genotype. 3.6. Recessive TRD Pattern Recessive TRD pattern (i.e., lethality only in homozygous state) was targeted on pseu- doautosomal SNPs (in both male- and female-offspring) and heterosomal SNPs/haplotypes (only female-offspring) on sex chromosomes. Whereas no recessive TRD patterns were found in SNP-by-SNP analyses, haplotype analyses revealed eight haplotypes (Tables 4 and S6). These regions were detected with the genotypic TRD model [3,37] exhibiting additive-TRD ranging from −0.34 to −0.56. The unfavorable effects of the reces- sive alleles were compensated by positive dominance-TRD effects (ranging from 0.21 to 0.28) in heterozygous offspring, allowing for their survival. The corresponding number of non-observed homozygous female-offspring was between 5 to 8 (Table 4). In comparation to autosomal chromosomes, Holstein haplotype 3 was initially identified by 7 non-observed homozygous offspring from heterozygous carrier sire by heterozygous carrier maternal grandsire matings [10]. Genes 2022, 13, 2322 13 of 23 Table 4. Haplotype windows with recessive TRD pattern on heterosomal part of Holstein X-chromosome. Coordinates in Base Pairs AA × AB 3 BB × AB (n. SNPs 1 ) n◦ Hetero 2 Dams 4 Frequency % TRD Parameters AA AB AB BB 307,557: 854,184 (10) 247.00 2 8 300 281 0.40 αg = −0.34 (BF = 103.35 ); δg = 0.21; (BF = 103.17 ) 14,897,141:16,123,681 (10) 452.00 1 8 559 579 0.60 αg = −0.47 (BF = 1012.96 ); δg = 0.21; (BF = 106.64 ) 23,895,384:27,387,980 (20) 541.00 1 7 674 634 0.90 αg = −0.40 (BF = 1012.08 ); δg = 0.24; (BF = 1010.46 ) 30,116,443:30,459,667 (4) 574.00 0 6 683 689 0.90 αg = −0.51 (BF = 1019.6 ); δg = 0.25; (BF = 1011.37 ) 48,750,092:52,808,614 (20) 592.00 0 5 647 669 1.10 αg = −0.50 (BF = 1017.26 ); δg = 0.23; (BF = 108.99 ) 59,564,066:61,635,072 (20) 525.00 0 8 550 549 1.80 αg = −0.56 (BF = 1019.61 ); δg = 0.28; (BF = 1011.44 ) 109,041,641:110,425,524 (20) 438.00 1 8 513 539 1.20 αg = −0.48 (BF = 1012.4 ); δg = 0.21; (BF = 105.82 ) 127,431,846:129,887,133 (20) 406.00 2 9 573 551 0.80 αg = −0.37 (BF = 108.37 ); δg = 0.21; (BF = 106.78 ) BF: Bayes factor, αg: additive-TRD, δg : dominance-TRD. 1 Number of SNPs on the haplotype window. 2 Number of heterozygous dams. 3 Sire × dam mating genotypes. 4 Female- offspring genotype. Genes 2022, 13, 2322 14 of 23 3.7. Functional Impact around the Pseudoautosomal Boundaries Identified Using Opposite TRD Patterns between Male and Female Offspring Among the SNPs with highest opposite sire-TRD between male- and female-offspring, the 20 bp sequence around BTX:143,928,317 SNP was identified within the LOC790278 pseu- dogene on the Y-chromosome, which has a corresponding version on the X-chromosome (SHROOM2 gene). The pseudogene LOC790278 is mapped on the heterosomal region of Y-chromosome (based on the Btau 4.6.1 annotation), before the SINE element responsible to reduce the homology between the X- and Y-chromosomes in the boundary between the heterosomal and pseudoautosomal [38]. Additionally, the region with the highest homology between the X-chromosome region and the Y-chromosome used in the BLAST analysis is also located in the predicted heterosomal part of the Y-chromosome. This region comprised 148.3 Kb (BTX:144,251,086–144,399,468) and has more than 99% similarity with the Y-chromosome (Supplementary Table S3). It is important to highlight that, as shown in Figure 2, even the sequential structure (regarding the gene order) of the homologous region and the BTX:143,928,317 SNP remains in the predicted heterosomal region of Y-chromosome. Within this region of homology between the pseudoautosomal region of the X-chromosome and predicted heterosomal region of the Y-chromosome, no TRD was detected due to the absence of markers exactly in this genome range, as shown in Figure 2C. As described by Raudsepp et al. (2012) [27], the SHROOM2 is a X-specific gene in cattle. However, following the gene annotation from NCBI genome data viewer, this gene is mapped to the boundary between the pseudoautosomal and the heterosomal regions, due to the location of GPR143, the gene responsible to define the beginning of the pseudoautosomal region [38]. The results obtained here indicate a possible miss-assembly of the pseudoautosomal region of Y chromosome in Btau 4.6.1 or a complex rearrangement of the pseudoautosomal Y region, which must be better studied due to the functional importance of these regions. Additionally, the functional interpretation of the candidate markers in these regions can be affected since the causal variant could be mapped to both chromosomes. 3.8. Enriched Metabolic Pathways, Diseases and Biological Functions Associated with the Positional Candidate Genes Displaying Deviations from Mendelian Inheritance on the X-Chromosome Functional annotation of positional candidate genes displaying deviations from Mendelian inheritance on the bovine X-chromosome was performed for the 55 detected regions: 15 SNPs with specific sex-offspring sire TRD (pseudoautosomal region), 12 SNPs with dam-TRD, 20 haplotypes with dam-TRD and 8 haplotypes with recessive TRD (het- erosomal region). A total of 408 genes were annotated in an interval of 200 Kb (100 Kb up- and downstream) for the individual SNPs or within the haplotype coordinates. A total of 113 genes were annotated in the interval harboring the individual SNPs (among them, 60 genes had TRD signals within the gene coordinates) and 296 within the haplotype coordinates. It is important to highlight that some genes were annotated in both SNP and haplotype approaches. The human and mouse orthologs were used for those genes without an assigned gene symbol in the bovine reference genome. A total of 295 genes on TRD regions (out of the initial 408) were annotated. This list of 295 genes was used to perform the functional analysis using the IPA (Ingenuity Pathways Analysis) software to identify the enriched canonical metabolic pathways (Supplementary Table S7), diseases, and bio- logical functions (Supplementary Table S8). The top 10 metabolic pathways, diseases, and functions identified are presented in Table 5. Despite no metabolic pathways were signifi- cant after adjustment for multiple testing (FDR < 0.05), the following canonical metabolic pathways were in the list of the top 10 enriched pathways at p-values < 0.05: citrulline degradation, spermine biosynthesis, phosphoribosyl diphosphate biosynthesis I, arginine biosynthesis IV, proline biosynthesis II (from arginine), urea cycle, pentose phosphate path- way (non-oxidative branch), glycine cleavage complex, phosphatidylcholine biosynthesis I, and acetyl-CoA biosynthesis I (pyruvate dehydrogenase complex). Among these metabolic pathways, the spermine biosynthesis plays a crucial role in male fertility. The gene asso- ciated with this pathway, in the current study, was the SMS (Spermine Synthetase) gene. Genes 2022, 13, 2322 15 of 23 Deletions in the X-chromosomal loci harboring the mice orthologous (SpmS) are associated with several impaired development phenotypes, such as sterility and profound postnatal growth retardation [39]. Another enriched pathway was the phosphatidylcholine biosyn- thesis, which was associated with the PCYT1B gene in this study, which had a haplotype exhibiting TRD. The disruption of function of the first enzyme in the phosphatidylcholine biosynthesis, the choline kinase α, causes embryonic lethality in mice [40]. Female knockout mice for the PCYT1B gene exhibited ovarian tissue disorganization with progressive loss of follicular formation and oocyte maturation [41]. Table 5. Top 10 canonical metabolic pathways and diseases and functions identified using the list of positional candidate genes mapped around (100 Kb up- and downstream) from candidate SNPs or within haplotype coordinates showing significant TRD. p-Value False Discovery Rate Genes Canonical Metabolic Pathway Citrulline Degradation 0.0072 0.603 OTC Spermine Biosynthesis 0.0145 0.603 SMS PRPP Biosynthesis I 0.0282 0.603 PRPS1 Arginine Biosynthesis IV 0.0427 0.603 OTC Proline Biosynthesis II 0.0427 0.603 OTC (from Arginine) Urea Cycle 0.0427 0.603 OTC Pentose Phosphate Pathway 0.0427 0.603 TKTL1 (Non-oxidative Branch) Glycine Cleavage Complex 0.0427 0.603 GCSH Phosphatidylcholine 0.0490 0.603 PCYT1B Biosynthesis I Acetyl-CoA Biosynthesis I (Pyruvate 0.0490 0.603 PDHA1 Dehydrogenase Complex) Diseases and Functions ADGRG2, ALG13, AMELX, AMER1, ANOS1, AP1S2, ARHGEF9, ARMCX5GPRASP2/GPRASP2, ARR3, BTK, CA5B, CDKL5, CNKSR2, COL4A6, FANCB, FGF13, FLNA, FRMD7, GLA, GPR143, HNRNPH2, IGSF1, KDM6A, X-linked hereditary disease 3.21 × 10−47 1.6 × 10−44 KIF4A, MBTPS2, MECP2, MID2, mir221, NDUFA1, OPN1LW, OTC, PCDH19, PDHA1, PDK3, PHEX, PHKA2, PIGA, PIH1D3, PLP1, POLA1, PRPS1, RNF113A, RPS6KA3, RS1, SERPINA7, SMPX, SMS, SRPX2, TIMM8A, TSPAN7, UBE2A, ZC4H2, ZNF674 AP1S2, CNKSR2, HNRNPH2, KDM6A, KIF4A, MECP2, MID2, PDHA1, RNF113A, X-linked mental retardation 1.47 × 10−15 3.7 × 10−13 RPS6KA3, SMS, SRPX2, TSPAN7, UBE2A, ZC4H2, ZNF674 Genes 2022, 13, 2322 16 of 23 Table 5. Cont. p-Value False Discovery Rate Genes ARMCX5GPRASP2/GPRASP2, COL4A6, X-linked hearing loss 1.29 × 10−8 2.2 × 10−6 PRPS1, SMPX Syndromic X-linked AP1S2, CNKSR2, HNRNPH2, KDM6A, 1.97 × 10−8 2.5 × 10−6 mental retardation MECP2, SMS, SRPX2, UBE2A, ZC4H2 AP1S2, CNKSR2, FLNA, HNRNPH2, KDM6A, KIF4A, MECP2, MID2, PDHA1, Familial mental retardation 4.28 × 10−8 4.3 × 10−6 RNF113A, RPS6KA3, SMS, SRPX2, TSPAN7, UBE2A, ZC4H2, ZNF674 X-linked nonsyndromic 1.46 × 10−6 1.2 × 10−4 COL4A6, PRPS1, SMPX sensorineural hearing loss AP1S2, CDKL5, CNKSR2, FLNA, HNRNPH2, KDM6A, KIF4A, MECP2, Autism or 1.95 × 10−6 1.4 × 10−4 MID2, PDHA1, RNF113A, RPS6KA3, intellectual disability SMS, SRPX2, TSPAN7, UBE2A, ZC4H2, ZNF674 AP1S2, CNKSR2, FLNA, HNRNPH2, Familial syndromic 2.57 × 10−6 1.6 × 10−4 KDM6A, MECP2, SMS, SRPX2, intellectual disability UBE2A, ZC4H2 Childhood onset 3.63 × 10−6 2.0 × 10−4 ARHGEF9, CDKL5, SRPX2 epilepsy syndrome AP1S2, ARHGEF9, CA5B, CNKSR2, FLNA, HNRNPH2, KDM6A, KIF4A, Cognitive impairment 3.49 × 10−5 1.7 × 10−3 MECP2, MID2, PDHA1, RNF113A, RPS6KA3, SMS, SRPX2, TSPAN7, UBE2A, VEGFD, ZC4H2, ZNF674 Using the list of 295 positional candidate genes, 500 diseases and biological functions were enriched (FDR < 0.05). Among the top 10 enriched diseases and biological functions identified, a strong association between the genes and X-linked syndromic phenotypes was observed. It is important to highlight that several X-linked syndromes caused by structural alterations (copy number variations, translocations, inversions, breakpoints, etc.; Table 5) display a combination of neurological and reproductive alterations, such as X-fragile, Turner, Klinefelter, Kallmann, and Rett syndromes [42–44]. In addition, among the other significant enriched diseases and biological functions, the enrichment analysis indicated that the candidate genes are associated with relevant biological processes associated with embryo and cellular development, such as cell viability, angiogenesis, development of cardiovascular tissue, migration of endothelial cells, growth of connective tissue, abnormal morphology of body cavity and neuritogenesis (Figure 5). The genes associated with these processes are mapped in regions where different TRD patterns were identified, all in heterosomal regions. The relationship between the TRD markers and the positional candidate genes is shown in Supplementary Table S9. Among these positional and functional candidate genes, only PLP1, DRP2 and BTK did not show a direct association with embryonic development or reproductive processes base on the literature review performed for each gene. different TRD patterns were identified, all in heterosomal regions. The relationship be- tween the TRD markers and the positional candidate genes is shown in Supplementary Genes 2022, 13, 2322 Table S9. Among these positional and functional candidate genes, only PLP1, DRP217andof 23 BTK did not show a direct association with embryonic development or reproductive pro- cesses base on the literature review performed for each gene. Figure Figure5.5.Interaction Interactionnetwork network between between the positional candidate the positional candidategenesgenesaround aroundTRD TRDregions regions onon bo- bovine vine X-chromosome and biological functions and diseases identified. For the most relevant biologi- X-chromosome and biological functions and diseases identified. For the most relevant biological cal functions, the edges between the genes and the processes were colored in order to highlight the functions, the edges between the genes and the processes were colored in order to highlight the connection. These biological functions were: cell viability (green), angiogenesis (pink), development ofconnection. Thesetissue cardiovascular biological (red),functions migrationwere: cell viabilitycells of endothelial (green), angiogenesis (blue), growth of(pink), development connective tissue of cardiovascular (cyan), tissue (red), migration abnormal morphology of endothelial of body cavity (orange) cells and (blue), growth of(purple). neuritogenesis connective Thetissue (cyan), positional abnormalgenes candidate morphology of bodyincavity were colored (orange) function of theand neuritogenesis predominant (purple). (highest |α|)TheTRDpositional observedcandidate for the associated genes were marker, colored where: pink represents in function dam-TRD among of the predominant (highest in|α|) heterosomal regionsfor TRD observed (RNF128, PIGA, the associated PLP1, BMX, marker, KDM6A, where: pinkCSNK1A1, represents VEGFD, dam-TRD DRP2 amongandinbta-mir-221); heterosomalyellow represents regions (RNF128,male- PIGA,and fe- PLP1, male-dam-TRD in heterosomal regions (FLNA, FGF13 and MECP2); blue represents BMX, KDM6A, CSNK1A1, VEGFD, DRP2 and bta-mir-221); yellow represents male- and female-dam- recessive TRD pattern TRD ininheterosomal heterosomalregions regions(FLNA, (TNMD, SLITRK4 FGF13 and BTK).blue and MECP2); Therepresents gene TSC22D3 was colored recessive in redin TRD pattern because two different TRD patterns were observed for the markers mapped in the interval used to heterosomal regions (TNMD, SLITRK4 and BTK). The gene TSC22D3 was colored in red because two annotate the gene. These TRD patterns were dam-TRD and recessive-TRD, both in heterosomal re- different TRD patterns were observed for the markers mapped in the interval used to annotate the gions. gene. These TRD patterns were dam-TRD and recessive-TRD, both in heterosomal regions. 3.9. 3.9.Functional FunctionalCandidate CandidateGenes GenesIdentified Identifiedininthe theCandidate CandidateRegions Regionswith withDeviations Deviationsfrom from Mendelian MendelianExpectations Expectations The Thegenes genesassociated associatedwith withregions regionswithwithdam-TRD dam-TRDwere: were:RNF128, RNF128,PIGA, PIGA,PLP1, PLP1,BMX, BMX, KDM6A, KDM6A,CSNK1A1, CSNK1A1,VEGFD, VEGFD,DRP2 DRP2andandbta-mir-221. bta-mir-221.TheTheRNF128 RNF128gene genewas wasidentified identifiedasas differentially differentially expressed [45]between expressed [45] betweenhigh high and and lowlow fertility fertility heifers heifers during during mid-luteal mid-luteal phase phase of theof the estrous estrous cycle. cycle. Mutations Mutations in thegene in the PIGA PIGAare gene are associated associated with nocturnal with nocturnal he- hemoglobin- uria [46]. Additionally, moglobinuria high rates [46]. Additionally, highofrates earlyofembryonic lethality early embryonic and low lethality andchimerism low chimer- (in the(in ism surviving animals) the surviving are observed animals) in mice are observed in knockedout mice knockedout for PIGA The BMX [47]. [47]. for PIGA The BMX gene has ahas gene crucial rolerole a crucial during the the during endothelial endothelial cellcell migration migration and angiogenesis and angiogenesis [48] [48]acting actingas asa downstream a downstream RAP1 RAP1effector in VEGF-induced effector in VEGF-induced endothelial endothelialcell activation [49]. [49]. cell activation An inser- An tion/deletion polymorphism in KDM6A is significantly associated with litter size and growth traits in goats, indicating a possible association between this gene and the survival ratio during embryonic development and post-natal development [50]. The CSNK1A1 acts as a regulator of DNA damage response in embryonic stem cells [51]. Additionally, the expression of CSNK1A1 is increased by the injection of miR-320 in mice oocytes [52]. Genes 2022, 13, 2322 18 of 23 Interestingly, knockout mice for mmu-miR-320 show a reduced proportion of oocytes that develop into two-cell and blastocyst stage embryos [52]. The VEGFD gene has a critical role during embryonic vasculogenesis and angiogenesis in zebrafish, reinforced by the microinjection of zVEGFD mRNA into one-cell-stage embryos, which results in severe misguidance of intersegmental vessels and abnormal connection between dorsal aorta and caudal vein [53]. It is important to highlight that the previously described gene, BMX, acts in the same VEGF-induced endothelial cell activation pathway. In addition, regarding the VEGF-induced endothelial cell activation pathway, the zebrafish ortholog of bta-mir-221 is required for the tip cell proliferation and migration in mosaic blood vessels during angiogenesis, a process which is regulated by VEGF [54]. Among the genes corresponding to male- and female-dam-TRD regions, 3 genes were highlighted: FLNA, FGF13 and MECP2. Mutations in FLNA gene cause total gene disruption, resulting in embryonic lethality, with severe cardiac structural defects, revealing the important role of FLNA during cardiac development [55]. In mice, the targeted ablation of FGF13 results in embryonic lethality [56]. The MECP2 gene is essential for embryonic development [57], acting in the differentiation of neuronal progenitors, which, in the absence of MECP2, do not show regular increase rates of the neuronal nuclei during development [58]. The positional and functional candidate genes associated with markers with significant recessive TRD patterns were: TNMD, SLITRK4 and BTK. Despite the association between TNMD and connective and cardiovascular tissue growth and differentiation, knockout mice show normal overall morphology at birth [59]. The knockout of Sirh7/Ldoc1 in mice exhibit an overproduction of placental progesterone and placental lactogen, resulting from abnormal placental cell differentiation and maturation, subsequently, causing delayed parturition [60]. The TSC22D3 gene was associated with markers displaying two different TRD patterns (dam-TRD and recessive-TRD). Deficient mice for TSC22D3 show severe testis dysplasia followed by infertility [61], reinforcing the association of this gene with male fertility-related biological processes. 3.10. Validation of Functional Candidate Genes Previously Identified in Independent Populations Using Genome-Wide Association Studies’ Meta-Analysis for Fertility Traits in Cattle In addition to the previous reported results, 15 out of the 408 positional candidate genes were previously identified as functional candidate genes for cattle male fertility traits in a systematic review performed by our group [62]. These genes were: MECP2, L1CAM, IRAK1, KIF4A, ATP6AP1, G6PD, POLA1, FLNA, IKBKG, HCFC1, EMD, IL13RA2, PAK3, DCX, and TEX11. Some of these genes are responsible for regulating important metabolic pathways or biological processes associated with fertility. Among these processes, it is important to highlight the control of the progression of spermatogenesis, control of ciliary activity, development of Sertoli cells, DNA integrity in spermatozoa, and home- ostasis of testicular cells. Additionally, the impact on fertility traits of these 15 functional and positional candidate genes was also evaluated in other species. The MECP2 gene is associated with the development of Rett syndrome in humans [63]. In human and mice, variants located in the MECP2 gene were previously associated with precocious puberty, undescended testes, and micropenis [64,65]. During the endometriosis devel- opment, L1CAM is considered a pathogenic factor in humans due to its overexpression in ovarian and endometrial carcinomas [66]. Interestingly, IRAK1 gene is differentially expressed in impaired spermatogenesis [67]. The function of KIF4A transcript is crucial for the chromosome positioning and spindle formation, playing an important role during chromosome segregation in mouse oocytes [68,69]. Additionally, KIF4A variants were already associated with uterine capacity in beef cattle [70]. The ATP6AP1 gene encodes the main accessory protein of V-ATPase, which is a crucial regulator of intra-organellar acidification [71]. The deficiency of ATP6AP1 is associated with immunodeficiency, hep- atopathy, and abnormal protein glycosylation in prenatal and postnatal stages [72]. In Genes 2022, 13, 2322 19 of 23 addition, ATP6AP1 gene was described, in zebrafish, as responsible to mediate the dorsal forerunner cell proliferation and the left-right asymmetry [73]. The GDP6 deficiency was suggested to be associated with reduced female fertility in humans [74]. The FLNA gene is involved with the regulation of the embryonic development of several tissues, such as the cardiac, neural and muscular tissues [55,75]. Consequently, FLNA gene plays a crucial role for embryo survival. The HCFC1 gene codifies a transcriptional co-regulator responsible to control cell proliferation, which loss-of-function variants were associated with disrupted neuronal and neural progenitor cells with a deleterious effect in humans and mice [76]. In humans, copy number variants at the PAK3 gene were found in men with Sertoli-cell-only syndrome, XY gonadal dysgenesis and premature ovarian failure [74]. The TEX11 gene was already broadly investigated in genetic association studies regarding male fertility in cattle and other species. This gene was associated with male infertility in humans [75] and in cattle [76]. Indeed, variants in TEX11 gene are responsible to explain up to 13% of the total genetic variance for scrotal circumference at 12 months in beef cattle [76]. 4. Conclusions Transmission ratio distortion analyses showing deviations from Mendelian inheritance revealed 149 SNPs displaying opposite sire-TRD between male and female offspring with strong magnitude at the beginning of the pseudoautosomal region and gradually reducing along this region until being null at the extreme of the chromosomes. This finding is evidence of the strong linkage of unique SNPs to specific-sex (Y- or X-) chromosomes and depicts the accumulation of recombination events across the pseudoautosomal region of the cattle genome. Moreover, haplotypes were found displaying recessive TRD patterns in female-offspring. Functional analyses and validation analysis using independent popula- tions and methodologies showed that the TRD regions identified in this study were related to key biological processes and molecular functions responsible for regulating processes such as spermatogenesis, development of sertoli cells, homeostasis of endometrium tissue and embryonic development. Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/genes13122322/s1, Figure S1: Linkage disequilibrium for pseudoautosomal SNPs (143,865,210-148,816,634 bp, BTX), heterosomal SNPs (29,245-137,034,696 bp, BTX) and autosomal SNPs (16,947-136,227,373, BTA1); Table S1. Map file of markers of X-chromosome with number of parent- offspring genotyped trios and Mendelian inconsistency rate; Table S2. Map file of the selected markers of heterosomal part of the X-chromosome for phasing and haplotype analyses; Table S3. Data from BLAST analyses of the X- and Y-chromosomes; Table S4. Full list with detailed information of specific sex-offspring sire-TRD in Holstein X-chromosome; Table S5. Full list with detailed information of significant dam-TRD among heterosomal part in Holstein X-chromosome from haplotype analyses; Table S6. Full list with detailed information of Haplotype windows with recessive TRD pattern on heterosomal part of Holstein X-chromosome; Table S7. The enriched canonical metabolic pathways obtained from functional analysis of the identified genes; Table S8. The diseases and biological functions pathways obtained from functional analysis of the identified genes; Table S9. Detailed TRD regions with the full information of the positional candidate genes; Table S10: Number of animals genotyped for each density of the SNP array. Author Contributions: Conceptualization: S.I.-L., J.C. and A.C.; Methodology: S.I.-L., J.C., P.A.S.F., A.S.-V. and A.C.; Analysis: S.I.-L. and P.A.S.F.; Investigation: S.I.-L., J.C., F.S.S., P.A.S.F., A.S.-V. and A.C.; Resources: A.C.; Data Curation: S.I.-L. and P.A.S.F.; Writing—Original Draft Preparation: S.I.-L., P.A.S.F. and A.S.-V.; Writing—Review and Editing: S.I.-L., P.A.S.F., A.S.-V., J.C., F.S.S. and A.C.; Supervision: A.C. and J.C.; Project Administration: A.C.; Funding Acquisition: A.C. All authors have read and agreed to the published version of the manuscript. Genes 2022, 13, 2322 20 of 23 Funding: This research (FDE.13.17) was supported by the Sustainable Beef and Forage Science Cluster funded by the Canadian Beef Cattle Check-Off, Beef Cattle Research Council (BCRC), Alberta Beef Producers, Alberta Cattle Feeders’ Association, Beef Farmers of Ontario, La Fédération des Productuers de bovins du Québec, and Agriculture and Agri-Food Canada’s Canadian Agricultural Partnership. The authors acknowledge the data provided by the Canadian Dairy Network (CND), a member of Lactanet (Guelph, Ontario, Canada). This study was also supported by NSERC (Natural Sciences and Engineering Research Council). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data that support the findings of this study were provided by the Canadian Dairy Network, a member of Lactanet (Guelph, ON, Canada), which were used under a material transfer agreement and thus are not publicly available. Conflicts of Interest: The authors declare no conflict of interest. References 1. Dunn, L.C. Evidence of evolutionary forces leading to the spread of lethal genes in wild populations of house mice. Proc. Natl. Acad. Sci. USA 1957, 43, 158–163. [CrossRef] [PubMed] 2. Khatib, H.; Monson, R.L.; Schutzkus, V.; Kohl, D.M.; Rosa, G.J.M.; Rutledge, J.J. Mutations in the STAT5A Gene Are Associated with Embryonic Survival and Milk Composition in Cattle. J. Dairy Sci. 2008, 91, 784–793. [CrossRef] [PubMed] 3. Casellas, J.; Gularte, R.J.; Farber, C.R.; Varona, L.; Mehrabian, M.; Schadt, E.E.; Lusis, A.J.; Attie, A.D.; Yandell, B.S.; Medrano, J.F. Genome Scans for Transmission Ratio Distortion Regions in Mice. Genetics 2012, 191, 247–259. [CrossRef] [PubMed] 4. Casellas, J.; Cañas-Álvarez, J.J.; González-Rodríguez, A.; Puig-Oliveras, A.; Fina, M.; Piedrafita, J.; Molina, A.; Díaz, C.; Baró, J.A.; Varona, L. Bayesian Analysis of Parent-Specific Transmission Ratio Distortion in Seven Spanish Beef Cattle Breeds. Anim. Genet. 2017, 48, 93–96. [CrossRef] [PubMed] 5. Abdalla, E.A.; Id-Lahoucine, S.; Cánovas, A.; Casellas, J.; Schenkel, F.S.; Wood, B.J.; Baes, C.F. Discovering Lethal Alleles across the Turkey Genome Using a Transmission Ratio Distortion Approach. Anim. Genet. 2020, 51, 876–889. [CrossRef] 6. Huang, L.O.; Labbe, A.; Infante-Rivard, C. Transmission Ratio Distortion: Review of Concept and Implications for Genetic Association Studies. Hum. Genet. 2013, 132, 245–263. [CrossRef] 7. Fishman, L.; Willis, J.H. A Novel Meiotic Drive Locus Almost Completely Distorts Segregation in Mimulus (Monkeyflower) Hybrids. Genetics 2005. [CrossRef] 8. VanRaden, P.M.; Olson, K.M.; Null, D.J.; Hutchison, J.L. Harmful Recessive Effects on Fertility Detected by Absence of Homozy- gous Haplotypes. J. Dairy Sci. 2011, 94, 6153–6161. [CrossRef] 9. Casellas, J.; Manunza, A.; Mercader, A.; Quintanilla, R.; Amills, M. A Flexible Bayesian Model for Testing for Transmission Ratio Distortion. Genetics 2014, 198, 1357–1367. [CrossRef] 10. Hoff, J.L.; Decker, J.E.; Schnabel, R.D.; Taylor, J.F. Candidate Lethal Haplotypes and Causal Mutations in Angus Cattle. BMC Genom. 2017, 18, 799. [CrossRef] 11. Howard, D.M.; Pong-Wong, R.; Knap, P.W.; Woolliams, J.A. Use of Haplotypes to Identify Regions Harbouring Lethal Recessive Variants in Pigs. Genet. Sel. Evol. 2017, 49, 47. [CrossRef] [PubMed] 12. Id-Lahoucine, S.; Cánovas, A.; Jaton, C.; Miglior, F.; Fonseca, P.A.S.; Sargolzaei, M.; Miller, S.; Schenkel, F.S.; Medrano, J.F.; Casellas, J. Implementation of Bayesian Methods to Identify SNP and Haplotype Regions with Transmission Ratio Distortion across the Whole Genome: TRDscan v.1.0. J. Dairy Sci. 2019, 102, 3175–3188. [CrossRef] [PubMed] 13. Naumova, A.K.; Leppert, M.; Barker, D.F.; Morgan, K.; Sapienza, C. Parental Origin-Dependent, Male Offspring-Specific Transmission-Ratio Distortion at Loci on the Human X Chromosome. Am. J. Hum. Genet. 1998, 102, 3175–3188. [CrossRef] [PubMed] 14. Naumova, A.K.; Olien, L.; Bird, L.M.; Slamka, C.; Fonseca, M.; Verner, A.E.; Wang, M.; Leppert, M.; Morgan, K.; Sapienza, C. Transmission-ratio Distortion of X Chromosomes among Male Offspring of Females with Skewed X-inactivation. Dev. Genet. 1995, 62, 1493–1499. [CrossRef] 15. Id-Lahoucine, S.; Casellas, J. Impact of Incomplete Pedigree Data and Independent Culling Level Pre-Selection on the Genetic Evaluation of Livestock: A Simulation Study on Lamb Growth. Livest. Sci. 2017, 198, 76–81. [CrossRef] 16. Kass, R.E.; Raftery, A.E. Bayes Factors. J. Am. Stat. Assoc. 1995, 198, 76–81. [CrossRef] 17. Ferdosi, M.H.; Kinghorn, B.P.; van der Werf, J.H.J.; Lee, S.H.; Gondro, C. Hsphase: An R Package for Pedigree Reconstruction, Detection of Recombination Events, Phasing and Imputation of Half-Sib Family Groups. BMC Bioinform. 2014, 90, 773–795. [CrossRef] 18. Fonseca, P.A.S.; Suárez-Vega, A.; Marras, G.; Cánovas, Á. GALLO: An R Package for Genomic Annotation and Integration of Multiple Data Sources in Livestock for Positional Candidate Loci. Gigascience 2020, 9, giaa149. [CrossRef] 19. Durinck, S.; Spellman, P.T.; Birney, E.; Huber, W. Mapping Identifiers for the Integration of Genomic Datasets with the R/ Bioconductor Package BiomaRt. Nat. Protoc. 2009, 9, giaa149. [CrossRef] Genes 2022, 13, 2322 21 of 23 20. Cardoso, T.F.; Cánovas, A.; Canela-Xandri, O.; González-Prendes, R.; Amills, M.; Quintanilla, R. RNA-Seq Based Detection of Differentially Expressed Genes in the Skeletal Muscle of Duroc Pigs with Distinct Lipid Profiles. Sci. Rep. 2017, 4, 1184–1191. [CrossRef] 21. Das, P.J.; Chowdhary, B.P.; Raudsepp, T. Characterization of the Bovine Pseudoautosomal Region and Comparison with Sheep, Goat, and Other Mammalian Pseudoautosomal Regions. Cytogenet. Genome Res. 2009, 126, 139–147. [CrossRef] [PubMed] 22. Couldrey, C.; Johnson, T.; Lopdell, T.; Zhang, I.L.; Littlejohn, M.D.; Keehan, M.; Sherlock, R.G.; Tiplady, K.; Scott, A.; Davis, S.R.; et al. Bovine Mammary Gland X Chromosome Inactivation. J. Dairy Sci. 2017, 100, 5491–5500. [CrossRef] [PubMed] 23. Johnson, T.; Keehan, M.; Harland, C.; Lopdell, T.; Spelman, R.J.; Davis, S.R.; Rosen, B.D.; Smith, T.P.L.; Couldrey, C. Short Communication: Identification of the Pseudoautosomal Region in the Hereford Bovine Reference Genome Assembly ARS- UCD1.2. J. Dairy Sci. 2019, 102, 3254–3258. [CrossRef] [PubMed] 24. Raudsepp, T.; Das, P.J.; Avila, F.; Chowdhary, B.P. The Pseudoautosomal Region and Sex Chromosome Aneuploidies in Domestic Species. Sex. Dev. 2012, 6, 72–83. [CrossRef] [PubMed] 25. Galtier, N. Recombination, GC-Content and the Human Pseudoautosomal Boundary Paradox. Trends Genet. 2004, 20, 347–349. [CrossRef] [PubMed] 26. McClure, M.C.; McCarthy, J.; Flynn, P.; McClure, J.C.; Dair, E.; O’Connell, D.K.; Kearney, J.F. SNP Data Quality Control in a National Beef and Dairy Cattle System and Highly Accurate SNP Based Parentage Verification and Identification. Front. Genet. 2018, 9, 84. [CrossRef] 27. Jeffreys, H. Theory of Probability. Philos. Sci. 1940, 7, 263–264. [CrossRef] 28. Rouyer, F.; Simmler, M.C.; Johnsson, C.; Vergnaud, G.; Cooke, H.J.; Weissenbach, J. A Gradient of Sex Linkage in the Pseudoauto- somal Region of the Human Sex Chromosomes. Nature 1986, 319, 291–295. [CrossRef] 29. Lien, S.; Szyda, J.; Schechinger, B.; Rappold, G.; Arnheim, N. Evidence for Heterogeneity in Recombination in the Human Pseudoautosomal Region: High Resolution Analysis by Sperm Typing and Radiation-Hybrid Mapping. Am. J. Hum. Genet. 2000, 66, 557–566. [CrossRef] 30. Hinch, A.G.; Altemose, N.; Noor, N.; Donnelly, P.; Myers, S.R. Recombination in the Human Pseudoautosomal Region PAR1. PLoS Genet. 2014, 10, e1004503. [CrossRef] 31. Naumova, A.; Sapienza, C. The Genetics of Retinoblastoma, Revisited. Am. J. Hum. Genet. 1994, 54, 264–273. [PubMed] 32. Bourc’His, D.; Bestor, T.H. Origins of Extreme Sexual Dimorphism in Genomic Imprinting. Cytogenet. Genome Res. 2006, 113, 36–40. [CrossRef] [PubMed] 33. Graves, J.A.M. Sex Chromosome Specialization and Degeneration in Mammals. Cell 2006, 124, 901–914. [CrossRef] [PubMed] 34. Casellas, J.; Id-Lahoucine, S.; Cánovas, A. Discriminating between Allele- and Genotype-Specific Transmission Ratio Distortion. Anim. Genet. 2020, 51, 847–854. [CrossRef] 35. Van Laere, A.S.; Coppieters, W.; Georges, M. Characterization of the Bovine Pseudoautosomal Boundary: Documenting the Evolutionary History of Mammalian Sex Chromosomes. Genome Res. 2008, 18, 1884–1895. [CrossRef] 36. Wang, X.; Ikeguchi, Y.; McCloskey, D.E.; Nelson, P.; Pegg, A.E. Spermine Synthesis Is Required for Normal Viability, Growth, and Fertility in the Mouse. J. Biol. Chem. 2004, 49, 51370–51375. [CrossRef] 37. Wu, G.; Aoyama, C.; Young, S.G.; Vance, D.E. Early Embryonic Lethality Caused by Disruption of the Gene for Choline Kinase α, the First Enzyme in Phosphatidylcholine Biosynthesis. J. Biol. Chem. 2008, 279, 51370–51375. [CrossRef] 38. Gunter, C.; Frank, M.; Tian, Y.; Murti, K.G.; Rehg, J.E.; Jackowski, S. Probucol Therapy Overcomes the Reproductive Defect in CTP: Phosphocholine Cytidylyltransferase B2 Knockout Mice. Biochim. Biophys. Acta Mol. Cell Biol. Lipids 2007, 283, 1456–1462. [CrossRef] 39. Anagnostopoulos, A.V. It’s a Knockout! Male Infertility and Neuropathology. Trends Genet. 2002, 18, 8–10. [CrossRef] 40. Cariboni, A.; Maggi, R.; Parnavelas, J.G. From Nose to Fertility: The Long Migratory Journey of Gonadotropin-Releasing Hormone Neurons. Trends Neurosci. 2007, 30, 638–644. [CrossRef] 41. Fullerton, G.; Hamilton, M.; Maheshwari, A. Should Non-Mosaic Klinefelter Syndrome Men Be Labelled as Infertile in 2009? Hum. Reprod. 2010, 25, 588–597. [CrossRef] [PubMed] 42. Killeen, A.P.; Morris, D.G.; Kenny, D.A.; Mullen, M.P.; Diskin, M.G.; Waters, S.M. Global Gene Expression in Endometrium of High and Low Fertility Heifers during the Mid-Luteal Phase of the Estrous Cycle. BMC Genom. 2014, 15, 234. [CrossRef] [PubMed] 43. Johnston, J.J.; Gropman, A.L.; Sapp, J.C.; Teer, J.K.; Martin, J.M.; Liu, C.F.; Yuan, X.; Ye, Z.; Cheng, L.; Brodsky, R.A.; et al. The Phenotype of a Germline Mutation in PIGA: The Gene Somatically Mutated in Paroxysmal Nocturnal Hemoglobinuria. Am. J. Hum. Genet. 2012, 90, 295–300. [CrossRef] [PubMed] 44. Kawagoe, K.; Kitamura, D.; Okabe, M.; Taniuchi, I.; Ikawa, M.; Watanabe, T.; Kinoshita, T.; Takeda, J. Glycosylphosphatidylinositol- Anchor-Deficient Mice: Implications for Clonal Dominance of Mutant Cells in Paroxysmal Nocturnal Hemoglobinuria. Blood 1996, 87, 3600–3606. [CrossRef] 45. Pan, S.; An, P.; Zhang, R.; He, X.; Yin, G.; Min, W. Etk/Bmx as a Tumor Necrosis Factor Receptor Type 2-Specific Kinase: Role in Endothelial Cell Migration and Angiogenesis. Mol. Cell. Biol. 2002, 22, 7512–7523. [CrossRef] 46. Stoletov, K.V.; Terman, B.I. Bmx Is a Downstream Rap1 Effector in VEGF-Induced Endothelial Cell Activation. Biochem. Biophys. Res. Commun. 2004, 320, 70–75. [CrossRef] Genes 2022, 13, 2322 22 of 23 47. Cui, Y.; Yan, H.; Wang, K.; Xu, H.; Zhang, X.; Zhu, H.; Liu, J.; Qu, L.; Lan, X.; Pan, C. Insertion/Deletion within the KDM6A Gene Is Significantly Associated with Litter Size in Goat. Front. Genet. 2018, 9, 91. [CrossRef] 48. Puigvert, J.C.; Von Stechow, L.; Siddappa, R.; Pines, A.; Bahjat, M.; Haazen, L.C.J.M.; Olsen, J.V.; Vrieling, H.; Meerman, J.H.N.; Mullenders, L.H.F.; et al. Systems Biology Approach Identifies the Kinase Csnk1a1 as a Regulator of the DNA Damage Response in Embryonic Stem Cells. Sci. Signal. 2013, 6, ra5. [CrossRef] 49. Feng, R.; Sang, Q.; Zhu, Y.; Fu, W.; Liu, M.; Xu, Y.; Shi, H.; Xu, Y.; Qu, R.; Chai, R.; et al. MiRNA-320 in the Human Follicular Fluid Is Associated with Embryo Quality in Vivo and Affects Mouse Embryonic Development in Vitro. Sci. Rep. 2015, 5, 8689. [CrossRef] 50. Song, M.; Yang, H.; Yao, S.; Ma, F.; Li, Z.; Deng, Y.; Deng, H.; Zhou, Q.; Lin, S.; Wei, Y. A Critical Role of Vascular Endothelial Growth Factor D in Zebrafish Embryonic Vasculogenesis and Angiogenesis. Biochem. Biophys. Res. Commun. 2007, 357, 924–930. [CrossRef] 51. Nicoli, S.; Knyphausen, C.P.; Zhu, L.J.; Lakshmanan, A.; Lawson, N.D. MiR-221 Is Required for Endothelial Tip Cell Behaviors during Vascular Development. Dev. Cell 2012, 22, 418–429. [CrossRef] [PubMed] 52. Feng, Y.; Chen, M.H.; Moskowitz, I.P.; Mendonza, A.M.; Vidali, L.; Nakamura, F.; Kwiatkowski, D.J.; Walsh, C.A. Filamin a (FLNA) Is Required for Cell-Cell Contact in Vascular Development and Cardiac Morphogenesis. Proc. Natl. Acad. Sci. USA 2006, 103, 19836–19841. [CrossRef] [PubMed] 53. Puranam, R.S.; He, X.P.; Yao, L.; Le, T.; Jang, W.; Rehder, C.W.; Lewis, D.V.; McNamara, J.O. Disruption of Fgf13 Causes Synaptic Excitatory–Inhibitory Imbalance and Genetic Epilepsy and Febrile Seizures Plus. J. Neurosci. 2015, 35, 8866–8881. [CrossRef] [PubMed] 54. Tate, P.; Skarnes, W.; Bird, A. The Methyl-CpG Binding Protein MeCP2 Is Essential for Embryonic Development in the Mouse. Nat. Genet. 1996, 12, 205–208. [CrossRef] [PubMed] 55. Yazdani, M.; Deogracias, R.; Guy, J.; Poot, R.A.; Bird, A.; Barde, Y.A. Disease Modeling Using Embryonic Stem Cells: MeCP2 Regulates Nuclear Size and RNA Synthesis in Neurons. Stem Cells 2012, 12, 205–208. [CrossRef] 56. Docheva, D.; Hunziker, E.B.; Fässler, R.; Brandau, O. Tenomodulin Is Necessary for Tenocyte Proliferation and Tendon Maturation. Mol. Cell Biol. 2005, 30, 2128–2139. [CrossRef] 57. Naruse, M.; Ono, R.; Irie, M.; Nakamura, K.; Furuse, T.; Hino, T.; Oda, K.; Kashimura, M.; Yamada, I.; Wakana, S.; et al. Sirh7/Ldoc1 Knockout Mice Exhibit Placental P4 Overproduction and Delsayed Parturition. Development 2014, 141, 4763–4771. [CrossRef] 58. Suarez, P.E.; Rodriguez, E.G.; Soundararajan, R.; Mérillat, A.M.; Stehle, J.C.; Rotman, S.; Roger, T.; Voirol, M.J.; Wang, J.; Gross, O.; et al. The Glucocorticoid-Induced Leucine Zipper (Gilz/Tsc22d3-2) Gene Locus Plays a Crucial Role in Male Fertility. Mol. Endocrinol. 2012, 26, 1000–1013. [CrossRef] 59. De Souza Fonseca, P.A.; dos Santos, F.C.; Lam, S.; Suárez-Vega, A.; Miglior, F.; Schenkel, F.S.; de Diniz, L.A.F.; Id-Lahoucine, S.; Carvalho, M.R.S.; Cánovas, A. Genetic Mechanisms Underlying Spermatic and Testicular Traits within and among Cattle Breeds: Systematic Review and Prioritization of GWAS Results. J. Anim. Sci. 2018, 96, 4978–4999. [CrossRef] 60. Dunn, H.G.; MacLeod, P.M. Rett Syndrome: Review of Biological Abnormalities. Can. J. Neurol. Sci. 2001, 28, 16–29. [CrossRef] 61. Sanmann, J.N.; Bishay, D.L.; Starr, L.J.; Bell, C.A.; Pickering, D.L.; Stevens, J.M.; Kahler, S.G.; Olney, A.H.; Schaefer, G.B.; Sanger, W.G. Characterization of Six Novel Patients with MECP2 Duplications Due to Unbalanced Rearrangements of the X Chromosome. Am. J. Med. Genet. A 2012, 158, 1285–1291. [CrossRef] [PubMed] 62. Tsuji-Hosokawa, A.; Matsuda, N.; Kurosawa, K.; Kashimada, K.; Morio, T. A Case of MECP2 Duplication Syndrome with Gonadotropin-Dependent Precocious Puberty. Horm. Res. Paediatr. 2017, 87, 271–276. [CrossRef] [PubMed] 63. Finas, D.; Huszar, M.; Agic, A.; Dogan, S.; Kiefel, H.; Riedle, S.; Gast, D.; Marcovich, R.; Noack, F.; Altevogt, P.; et al. L1 Cell Adhesion Molecule (L1CAM) as a Pathogenetic Factor in Endometriosis. Hum. Reprod. 2008, 23, 1053–1062. [CrossRef] [PubMed] 64. Spiess, A.N.; Feig, C.; Schulze, W.; Chalmel, F.; Cappallo-Obermann, H.; Primig, M.; Kirchhoff, C. Cross-Platform Gene Expression Signature of Human Spermatogenic Failure Reveals Inflammatory-like Response. Hum. Reprod. 2007, 22, 2936–2946. [CrossRef] [PubMed] 65. Bringmann, H.; Skiniotis, G.; Spilker, A.; Kandels-Lewis, S.; Vernos, I.; Surrey, T. A Kinesin-like Motor Inhibits Microtubule Dynamic Instability. Science 2004, 303, 1519–1522. [CrossRef] 66. Tang, F.; Pan, M.H.; Lu, Y.; Wan, X.; Zhang, Y.; Sun, S.C. Involvement of Kif4a in Spindle Formation and Chromosome Segregation in Mouse Oocytes. Aging Dis. 2018, 9, 623–633. [CrossRef] 67. Neupane, M.; Geary, T.W.; Kiser, J.N.; Burns, G.W.; Hansen, P.J.; Spencer, T.E.; Neibergs, H.L. Loci and Pathways Associated with Uterine Capacity for Pregnancy and Fertility in Beef Cattle. PLoS ONE 2017, 12, e0188997. [CrossRef] 68. Jansen, E.J.R.; Martens, G.J.M. Novel Insights into V-ATPase Functioning: Distinct Roles for Its Accessory Subunits ATP6AP1/Ac45 and ATP6AP2/(pro) Renin Receptor. Curr. Protein Pept. Sci. 2012, 13, 124–133. [CrossRef] 69. Jansen, E.J.R.; Timal, S.; Ryan, M.; Ashikov, A.; Van Scherpenzeel, M.; Graham, L.A.; Mandel, H.; Hoischen, A.; Iancu, T.C.; Raymond, K.; et al. ATP6AP1 Deficiency Causes an Immunodeficiency with Hepatopathy, Cognitive Impairment and Abnormal Protein Glycosylation. Nat. Commun. 2016, 13, 124–133. [CrossRef] 70. Gokey, J.J.; Dasgupta, A.; Amack, J.D. The V-ATPase Accessory Protein Atp6ap1b Mediates Dorsal Forerunner Cell Proliferation and Left-Right Asymmetry in Zebrafish. Dev. Biol. 2015, 7, 11600. [CrossRef] Genes 2022, 13, 2322 23 of 23 71. Gloria-Bottini, F.; Meloni, G.F.; MacMurray, J.P.; Cosmi, E.; Meloni, T.; Bottini, E. Evidence of Decreased Fertility in Women Carrying the Gene for G6PD Deficiency: A Study in the Sardinian Population. In The Biodemography of Human Reproduction and Fertility; Springer: Cham, Switzerland, 2003. 72. Lian, G.; Lu, J.; Hu, J.; Zhang, J.; Cross, S.H.; Ferland, R.J.; Sheen, V.L. Filamin a Regulates Neural Progenitor Proliferation and Cortical Size through Wee1-Dependent Cdk1 Phosphorylation. J. Neurosci. 2012, 32, 7672–7684. [CrossRef] [PubMed] 73. Jolly, L.A.; Nguyen, L.S.; Domingo, D.; Sun, Y.; Barry, S.; Hancarova, M.; Plevova, P.; Vlckova, M.; Havlovicova, M.; Kalscheuer, V.M.; et al. HCFC1 Loss-of-Function Mutations Disrupt Neuronal and Neural Progenitor Cells of the Developing Brain. Hum. Mol. Genet. 2015, 24, 3335–3347. [CrossRef] [PubMed] 74. Tüttelmann, F.; Simoni, M.; Kliesch, S.; Ledig, S.; Dworniczak, B.; Wieacker, P.; Röpke, A. Copy Number Variants in Patients with Severe Oligozoospermia and Sertoli-Cell-Only Syndrome. PLoS ONE 2011, 6, e19426. [CrossRef] [PubMed] 75. Yatsenko, A.N.; Georgiadis, A.P.; Röpke, A.; Berman, A.J.; Jaffe, T.; Olszewska, M.; Westernströer, B.; Sanfilippo, J.; Kurpisz, M.; Rajkovic, A.; et al. X-Linked TEX11 Mutations, Meiotic Arrest, and Azoospermia in Infertile Men. N. Engl. J. Med. 2015, 372, 2097–2107. [CrossRef] [PubMed] 76. De Camargo, G.M.F.; Porto-Neto, L.R.; Kelly, M.J.; Bunch, R.J.; McWilliam, S.M.; Tonhati, H.; Lehnert, S.A.; Fortes, M.R.S.; Moore, S.S. Non-Synonymous Mutations Mapped to Chromosome X Associated with Andrological and Growth Traits in Beef Cattle. BMC Genom. 2015, 16, 384. [CrossRef]