ROH-Genomic characterization of autozygosity and recent inbreeding trends in all major breeds of US dairy cattle

Maintaining a genetically diverse dairy cattle population is critical to preserving adaptability to future breeding goals and avoiding declines in fitness. This study characterized the genomic landscape of autozygosity and assessed trends in genetic diversity in 5 breeds of US dairy cattle. We analyzed a sizable genomic data set containing 4,173,679 pedigreed and genotyped animals of the Ayrshire, Brown Swiss, Guernsey, Holstein, and Jersey breeds. Runs of homozygosity (ROH) of 2 Mb or longer in length were identified in each animal. The within-breed means for number and the combined length of ROH were highest in Jerseys (62.66 ± 8.29 ROH and 426.24 ± 83.40 Mb, respectively; mean ± SD) and lowest in Ayrshires (37.24 ± 8.27 ROH and 265.05 ± 85.00 Mb, respectively). Short ROH were the most abundant, but moderate to large ROH made up the largest proportion of genome autozygosity in all breeds. In addition, we identified ROH islands in each breed. This revealed selection patterns for milk production, productive life, health, and reproduction in most breeds and evidence for parallel selective pressure for loci on chromosome 6 between Ayrshire and Brown Swiss and for loci on chromosome 20 between Holstein and Jersey. We calculated inbreeding coefficients using 3 different approaches, pedigree-based (F PED ), marker-based using a genomic relationship matrix (F GRM ), and segment-based using ROH (F ROH ). The average inbreeding coefficient ranged from 0.06 in Ayrshires and Brown Swiss to 0.08 in Jerseys and Holsteins using F PED , from 0.22 in Holsteins to 0.29 in Guernsey and Jerseys using F GRM , and from 0.11 in Ayrshires to 0.17 in Jerseys using F ROH . In addition, the effective population size at past generations (5–100 generations ago), the yearly of inbreeding, the current and historical trends of genetic diversity. We found a historical trend of decreasing effective population size in the last 100 generations in all breeds and breed differences in the effect of the recent implementation of genomic selection on inbreeding accumulation.


INTRODUCTION
The development of dense SNP arrays has enabled the identification of patterns of homozygosity in the genome. The detection of continuous stretches of the genome with SNP in a homozygous state, known as runs of homozygosity (ROH), was first performed in humans (Gibson et al., 2006) and subsequently in domestic animal populations (Boyko et al., 2010;Ferenčaković et al., 2011;Bosse et al., 2012) to determine the current and past demographic histories of these populations. These ROH are generated when chromosomal segments inherited from each parent are derived from a common ancestor. Therefore, the characteristics of ROH can give insight into population-level dynamics of selection and inbreeding. The length of an ROH is expected to be related to the number of generations that have passed since the common ancestor, with short ROH being derived from a distant past and long ROH being of recent origin (Howrigan et al., 2011). The distribution of ROH has been extensively used to map the selective histories of populations with the knowledge that ROH are not randomly distributed across the genome (Zhang et al., 2015), and individuals from the same population share similar patterns of ROH distribution (Bosse et al., 2012). Therefore, continuous regions of the genome with a high proportion of individuals exhibiting ROH, known as ROH islands, can be used as signatures of selection. Detecting selection signatures using ROH-based metrics has uncovered genomic regions under selection within and across populations in dairy cattle. Highly autozygous genomic regions in dairy have been found to contain genes associated with economically relevant traits such as milk production (Kim et al., 2013(Kim et al., , 2015Howard et al., 2015) and udder health (Kim et al., 2015;Peripolli et al., 2018).
Genomic information has allowed the determination of an individual's realized level of inbreeding rather than the expectation obtained using pedigree records. Additionally, although inbreeding calculated from pedigree information may suffer from the nonavailability or inaccuracy of genealogical records, genomic-based inbreeding can be determined for individuals with genotype data regardless of available records. Measures of genomic inbreeding broadly fall into either markerbased or segment-based approaches. Marker-based genomic inbreeding is traditionally obtained by subtracting one from the diagonal of genomic relationship matrices. However, several methods exist to construct the genomic relationship matrix, and inbreeding coefficients obtained from them can vary in range and interpretation (Villanueva et al., 2021). In the last decade, segment-based genomic inbreeding coefficients such as the proportion of the genome composed of ROH have become consistent measures of inbreeding in human, plant, and animal populations. Furthermore, inbreeding coefficients based on ROH have been found to measure autozygosity reliably and have high correlations with traditional measures of inbreeding (Howrigan et al., 2011;Forutan et al., 2018).
Over the last 25 yr, significant advances have been made in using genomic information to predict the genetic merit of animals for use in selection (Nejati-Javaremi et al., 1997;Meuwissen et al., 2001). Genomic selection was implemented in North American dairy populations a little over a decade ago, with the first official genomic evaluations of Holstein (HO) and Jersey (JE) breeds released in 2009 . This almost immediately resulted in reductions in generation intervals and considerable improvement to many dairy traits, especially for low heritability traits such as reproductive efficiency, longevity, and health (García-Ruiz et al., 2016). Before its widespread implementation, there was much speculation about how genomic selection would affect genetic diversity in cattle. For example, genomic information allows breeders to account for the Mendelian sampling term in predicting genetic merit. This was expected to lead to a reduction in the co-selection of sibs and consequently, inbreeding (Daetwyler et al., 2007;de Roos et al., 2011). However, at the same time, the ability to predict the breeding value of young bulls would decrease generation intervals and increase the yearly rate of inbreeding, especially if the accuracy of prediction in unproven animals is low (Schaeffer, 2006;de Roos et al., 2011). Recent studies have reported increases in the rate of inbreeding after the implementation of genomic selection in several dairy cattle populations (Doekes et al., 2018;Doublet et al., 2019;Makanjuola et al., 2020;Scott et al., 2021;Ablondi et al., 2022).
Our study leveraged a sizable data set consisting of numerically large (HO and JE) and small [Ayrshire (AY), Brown Swiss (BS), and Guernsey (GU)] US dairy cattle populations with the aim to (1) characterize genome-wide and region-specific autozygosity, (2) assess the current state of pedigree and genomic inbreeding, and (3) pinpoint historical and recent changes in genetic diversity.

MATERIALS AND METHODS
Because no live animals were used in this study, Institutional Animal Care and Use Committee approval was not required.

Animals and Data
Imputed genotypes were obtained for 76,389 autosomal SNP of 4,173,679 animals of the AY, BS, GU, HO, and JE US dairy cattle breeds. Pedigree information on these animals was also made available. Methodology on variant selection and imputation are described by VanRaden et al. (2017).
Quality control of genomic data was performed in PLINK 1.9 (Chang et al., 2015) in 2 steps. First, we removed animals within each breed with a call rate of less than 95% and SNP with a minor allele frequency of less than 0.01% or a Mendelian error rate of more than 10%. Afterward, we removed SNP with less than 95% across-breed call rate. This resulted in a total of 62,546 SNP remaining for further analysis. In addition, we performed quality control of the pedigree data of genotyped animals by removing animals with several complete generations less than 3 generations and a pedigree completeness index that takes 5 generations into account less than 0.90. The pedigree completeness index was calculated according to MacCluer et al. (1983). We opted for creating 2 data sets from the available data, a complete data set that includes all animals that passed genomic quality control regardless of pedigree completeness (data set A) and a smaller data set containing those animals from data set A that passed pedigree quality control (i.e., had a sufficiently complete pedigree; data set B). We conducted all following analyses using data set A unless expressly noted (i.e., in the estimation of pedigree inbreeding). For each breed, the number of animals in the pedigree file, the 8958 number of genotyped animals, and the number of animals in each in the data set after quality control are shown in Table 1. The number of animals born per year in each data set and for each breed can be found in Supplemental Figure S1 (https: / / doi .org/ 10 . 6084/ m9 .figshare .19232739 .v2;Lozada-Soto et al., 2022).

Runs of Homozygosity
The ROH were detected using a sliding window approach in PLINK 1.9 (Chang et al., 2015). The parameters chosen to define an ROH included: (1) a window size of 20 SNP, (2) at most one heterozygote call in a window, (3) at most 2 missing calls in a window, (4) minimum physical length of an ROH of 2 Mb, (5) a maximum gap of 500 kb between consecutive SNP, (6) a minimum SNP density of at least one SNP per 100 kb, and (7) a minimum number of 60 SNP to declare an ROH. The number of SNP in a window was set at 20 according to previous results (Forutan et al., 2018). Forutan et al. (2018) found that setting the number of SNP in a window to 20 produced closer estimates to actual inbreeding under several simulated scenarios over choices of 35 or 50 SNP. Furthermore, the maximum gap between SNP and the SNP density used was chosen according to results published by Meyermans et al. (2020), which showed that these are sensible choices to achieve sufficient genome coverage in medium-density SNP data. Finally, the minimum number of SNP (l) in an ROH was chosen according to the following formula given by Purfield et al. (2012) based on a method proposed by Lencz et al. (2007) to minimize the number of ROH called by chance: where a is the false positive rate (set at 0.05 in this study), n i is the number of genotyped individuals, n s is the number of SNP, and het is the average genome-wide SNP heterozygosity. We calculated the necessary l to be 59.80 SNP in our study. The detected ROH were then grouped into 4 length categories, 2 to 4 Mb, 4 to 8 Mb, 8 to 16 Mb, and 16 Mb or larger, to calculate the average number of ROH and the average combined length of the genome in ROH from each length category.

Selection Signatures Based on Reduced Local Variability
Within breed, a coefficient of autozygosity for every marker (F L ) was calculated as the proportion of animals in which the marker was located inside an ROH. We identified noteworthy SNP as those with F L at or above the 99.5th percentile for all genome-wide markers. Adjacent noteworthy SNP were joined into larger segments known as ROH islands. The ROH islands represent regions of reduced local variability that serve as genomic signatures of past selective history. Genes located inside ROH islands were annotated using the "btaurus_gene_ensembl" database (ARS-UCD1.2 genome assembly; https: / / www .ensembl .org/ Bos _taurus/ Info/ Index) with the "biomaRt" R package (v.2.48.3;Durinck et al., 2009). We used the GALLO R package (Fonseca et al., 2020) to query the Animal QTLdb (https: / / www .animalgenome .org/ cgi -bin/ QTLdb/ index) for QTL that have been previously identified in the regions of interest. Trait enrichment analysis was performed on the annotated QTL by chromosome. Significantly enriched traits had an FDR-adjusted P-value lower than 0.05.

Measures of Inbreeding
Inbreeding was calculated using 3 distinct methods, a pedigree-based approach, a genomic marker-based approach, and a genomic segment-based approach. A pedigree inbreeding coefficient (F PED ) was calculated for the animals in data set B with PEDIG (Boichard, 2002) using the tabular method as proposed by Van-Raden (1992). The marker-based genomic inbreeding coefficient (F GRM ) was calculated as G ii -1 where G is the diagonal of the genomic relationship matrix built according to the first method outlined by VanRaden (2008) and using fixed allele frequencies of 0.5. Finally, the segment-based genomic inbreeding coefficient (F ROH ) was calculated as the proportion of the autosomal genome covered by ROH.

Estimates of Past Effective Population Size
Effective population size t generations ago N e t ( ) was estimated using a genetic algorithm (Mitchell, 1998) implemented in the GONE software (Santiago et al., 2020). This method is based on the relationship between the linkage disequilibrium observed between pairs of SNP and N e , and allows for nonlinear changes in N e , such as population bottlenecks and expansions. The GONE program was run with default parameters, including unknown phase, an average rate of recombination of 0.01 cm/Mb, and a genetic distance correction based on Haldane's function. For this analysis we used data from the last year of data collection (2020); this included 198 AY, 405 GU, 1,294 BS, 369,000 HO, and 43,200 JE animals. Because the GONE software has a maximum allowable sample size of 1,800 individuals, in the case of HO and JE, we performed the analysis using subsets of the data, 205 and 24 subsets of 1,800 individuals for HO and JE, respectively. For these 2 breeds, N e t was then obtained as the average of the estimates across subsets.

Genetic Diversity in the Era of Genomic Selection
The yearly rate of inbreeding (ΔF yr ) and effective population size (N e ) was assessed for sires and dams born in 3 periods of interest: before the advent of genomic selection (GS; period 1, P1) from 2000 to 2009, during the implementation of GS (P2) from 2010 to 2014, and after the widespread adoption of GS from 2015 to 2018 (P3). The number of sires and dams in each period for all breeds are found in Supplemental

Characterization of ROH
Descriptive statistics for the number of ROH, the average length of ROH, and the combined length of the genome in ROH for each breed are presented in Table 2. The mean number of ROH ranged from 37.24 (AY) to 62.66 (JE), the mean average ROH length ranged from 6.74 (GU) to 7.54 (BS) Mb, and the mean combined length of ROH ranged from 265.05 (AY) to 426.24 (JE) Mb. The mean number and mean combined length of ROH by length category and breed are shown in Figure  1. We observed a decrease in the mean number of ROH in all 5 breeds with increasing segment length. Jersey had the highest, and AY had the lowest mean number of ROH in all length categories. The combined length of ROH of 2 to 4 Mb was smaller than for longer ROH in all breeds. The combined length of ROH in the 3 remaining length categories (4-8 Mb, 8-16 Mb, and 16 Mb or longer) were similar. We found the largest value for mean combined length in the 16 Mb or longer category for AY (83.99 Mb), the 8 to 16 Mb category for BS (113.83 Mb) and HO (94.97 Mb), and the 4 to 8 Mb category for GU (103.38 Mb) and JE (127.85 Mb).

Selection Signatures Uncovered in Regions of Reduced Local Variability
Manhattan plots of F L for all breeds are shown in Figure 2, and the genomic coordinates and gene content of ROH islands are found in Supplemental  76-45.56) and were identified in GU and JE. The latter was also the most gene-dense region with 167 genes annotated to this genomic region. We found a 3.44 Mb genomic region shared between AY and BS on chromosome 6 (BTA6: 86.76-90.21), which contained the most homozygous SNP in both breeds. We annotated a total of 26 genes to this shared region, most notably ADAMTS3, ANKRD17, GC, and NPFFR2. We also found a highly homozygous region shared between HO and JE, this time a 2.58-Mb region in chromosome 20 (BTA20: 25.59-28.17) containing 5 genes. The most notable genes in this shared region were ISL LIM homeobox 1 (ISL1) and pelota mRNA surveillance and ribosome rescue factor (PELO). Significant regions in GU did not overlap with any other breed.
The top 10 (based on P-value) enriched traits for each breed are shown in Figure 3. We found traits related to milk production to be enriched in all breeds. Milk yield, milk component yield, or milk component percentage were enriched in AY (QTL on BTA6), BS (BTA 5 and BTA6), GU (BTA6, BTA13, and BTA19), HO (BTA10 and BTA20), and JE (BTA2, BTA3, and BTA7). Other enriched milk production traits included lactation persistency in BS (BTA6), GU (BTA13), and HO (BTA20); milking speed in BS (BTA6), HO (BTA10), and JE (BTA7); milk color in GU (BTA15); and curd firmness in JE (BTA7). Enriched body measurement traits were found in all breeds except Holstein, including BW at birth in AY (BTA5), GU (BTA11), and JE (BTA3); BW at maturity in AY (BTA4); metabolic BW in GU (BTA13); height in BS (BTA16); and rump length in AY (BTA4). Feed conversion ratio was enriched for annotated QTL on BTA 16 in BS. The length of productive life was a trait that was enriched in all breeds except GU; this was for QTL on BTA6 in AY and BS and for QTL on BTA20 in HO and JE. Traits related to reproduction were enriched in all breeds except GU, including daughter pregnancy rate in AY (BTA6) and BS (BTA6); conception rate in AY (BTA6) and BS (BTA6); gestation length in AY (BTA5); inhibin level in AY (BTA5); scrotal circumference in AY (BTA5); dystocia in HO (BTA20); and calving index in JE (BTA20). Enriched udder and body conformation traits were found for AY, BS, and JE, including udder depth in AY (BTA6) and BS (BTA6); udder attachment and structure in JE (BTA20), angularity in AY (BTA6) and BS (BTA6); and rump angle in AY (BTA8). The sole enriched pigmentation trait was facial pigmentation in JE (BTA7). Traits related to health were enriched in all breeds. Many of these were related to the susceptibility to disease, including mastitis in AY (BTA6) and BS (BTA6); ketosis in AY (BTA6) and BS (BTA6); bovine tuberculosis in AY (BTA8); and bovine spongiform encephalopathy in GU (BTA6). Additional enriched health traits included immunoglobulin G level in AY (BTA4), HO (BTA20), and JE (BTA20); tick resistance in BS (BTA5) and JE (BTA3); and SCS in AY (BTA6) and BS (BTA6). One of the most interesting results was the high number of enriched traits related to meat and carcass found in GU, with around 45% of enriched traits belonging to this type. In this breed, we found QTL related to meat production on BTA6, BTA13, BTA15, and BTA19, with the majority of enriched traits associated with the content of SFA and MUFA.

Trends in Pedigree and Genomic Inbreeding
Within-breed means for animals born from 2000 to 2020 for pedigree (F PED ) and genomic inbreeding (F GRM and F ROH ) can be found in Table 3. The average inbreeding ranged from 0.06 (AY and BS) to 0.08 (HO and JE) using F PED , from 0.22 (HO) to 0.29 (GU and JE) using F GRM , and from 0.11 (AY) to 0.17 (JE) using F ROH . The Pearson correlation coefficient between inbreeding mea- sures (results not shown) for animals born from 2000 to 2020 ranged from high to very high in all breeds, with the correlation between F PED and F GRM ranging from 0.58 (GU) to 0.77 (HO), between F PED and F ROH ranging from 0.59 (GU) to 0.76 (HO), and between F GRM and F ROH ranging from 0.95 (GU) to 0.98 (HO). The within-breed average for F PED , F GRM , and F ROH by year of birth (2000-2020) is shown in Figure 4. We found a general trend of increased average inbreeding in each population for all 3 measures across time. However, we can observe substantial heterogeneity between breeds for the magnitude of change in inbreeding across time. For example, although JE animals have similar levels of pedigree inbreeding and higher average levels of genomic inbreeding than HO (  ing rates for GU dams were low and nonsignificant across periods and inbreeding measures. For AY dams, we found an increase in the ΔF yr after the implementation of GS, from a negative rate in P1 (∆F yr PED = −0.17%; ∆F yr GRM = −0.19%) to a positive rate in P2 (∆F yr PED = 0.25%; ∆F yr GRM = 0.28%

Effective Population Size
Effective population size estimates for each breed for the last 5 to 100 generations are shown in Figure 5, expressed in terms of calendar years, assuming one generation equals 5 calendar years. For both HO and JE, estimates of N e t were highly similar across data subsets, with the coefficient of variation calculated within generation ranging from 0.04 to 0.32 (average of 0.15) and from 0.02 to 0.21 (average of 0.09), respectively. Across all breeds, genetic diversity, as measured by N e, has di-  breeds in both magnitude and in the timing of population bottlenecks. The estimates for N e for GU and JE for 100 generations ago (circa 1520) were the lowest across all breeds, 2,640 and 1,792 individuals for GU and JE, respectively. For HO, the estimates of N e hovered around 4,000 individuals until 30 generations ago (circa 1870), after which we observed a fast decline, with estimates remaining under 500 individuals from 25 generations ago (circa 1895) onwards. The estimates of N e for AY and BS were the highest 100 generations ago, 14,024 and 14,595 individuals, respectively. However, although the N e of AY started slowly declining around 85 generations ago (circa 1595), the N e of BS remained above 10,000 individuals until about 30 generations ago (circa 1870), when a fast decline in genetic diversity was observed. After this, N e estimates for BS have remained under 500 individuals from 25 to 5 generations ago (circa 1895 to 1995), similar to what was observed for Holsteins.
We used the estimates of the yearly rate of inbreeding coupled with the generation intervals calculated at each period to estimate each population's effective population size in the 3 recent time periods (P1, P2, and P3). The generation intervals and estimates of N e are included in Supplemental Tables S3 and S4 (https: / / doi .org/ 10 .6084/ m9 .figshare .19232739 .v2; Lozada-Soto et al., 2022), respectively. This measure was only estimated when the yearly rate of inbreeding was positive and different from zero. The effective population size after the adoption of GS (P3) was found to be low for JE sires (46-51 individuals) and BS dams (31-64 individuals), and very low for HO sires (13-18 individuals), HO dams (20-25 individuals), and BS sires (12-18 individuals). The population with the largest estimate of N e at P3 was JE dams, with an effective population size ranging from 99 to 276 individuals.

Patterns of Homozygosity
Because selection is a significant force in shaping genome-wide homozygosity, we expect differences in the ROH landscapes of populations due to differences in their selective history. Multiple studies have identified differences in ROH characteristics between cattle populations (Purfield et al., 2012;Ferenčaković et al., 2013;Marras et al., 2015). We found substantial heterogeneity between breeds in terms of the number, the average length, and the combined length of ROH. We observed an inverse relationship between the number of segments and segment length, which is in line with previous findings (Forutan et al., 2018;Huson et al., Lozada-Soto et   2020). Forutan and colleagues (2018) reported relative frequencies of ROH that ranged from 43.5% for segments shorter than 2 Mb to 4.7% for segments longer than 16 Mb. Compared with the other breeds, the higher average number of short ROH (2-4 Mb) and the lower mean average ROH length found in JE and GU suggest more ancient origins to the current levels of homozygosity. This is consistent with the insular origin of these 2 breeds in the Channel Islands. Similarly, Huson et al. (2020) found a higher average number of ROH segments of short length (0.5-4 Mb) in JE and GU than HO cattle.

Genomic Signatures of Selection
The frequency of ROH in particular genomic regions allows us to infer and compare the selective history of populations. We found overlap between the significantly homozygous regions identified in our study and those previously identified in US populations of Holstein and Jersey cattle. Specifically, the ROH islands identified on BTA10 and BTA20 in HO and BTA2, BTA3, and BTA7 in JE overlapped with previously identified re-gions in these breeds (Kim et al., 2013(Kim et al., , 2015. Trait enrichment for QTL identified in these genomic regions suggest that historical selection for health and milk production might have contributed to the consistently high levels of autozygosity observed across studies.
We found the genomic coordinates of ROH islands to be mostly heterogeneous across populations. There were only 2 overlaps between ROH islands in different breeds, one on chromosome 6 between AY and BS and the other on chromosome 20 between HO and JE.
The 3.44 Mb overlapped region on BTA6 between AY and BS contained several genes with known functions associated with traits selected for in dairy cattle; these were ADAMTS3, ANKRD17, GC, and NPFFR2. These genes have been associated with productive life (ADAMTS3 and NPFFR2; Mészáros et al., 2014;Zhang et al., 2016), milk and milk component yield (ADAMTS3, ANKRD17, GC, and NPFFR2; Costa et al., 2019;Jiang et al., 2019), clinical mastitis (GC and NPFFR2; Sahana et al., 2014;Freebern et al., 2020), and subclinical ketosis (ANKRD17; Soares et al., 2021). Trait enrichment based on QTL annotated in this chromosome corroborated the high selective pressure for   loci on BTA6 affecting health, milk and milk component yield, and productive life in AY and BS. Current weights given to these trait types in selection indices created by the US Ayrshires Breeders' Association and the Brown Swiss Association range from 25 to 35% and 22 to 30% for milk component yield traits, 4% and 4 to 6% for health traits, and 0 and 8% for productive life, in AY (U.S. Ayrshire Breeders' Association, 2022a) and BS (Brown Swiss Association, 2022a), respectively. The ISL1 and PELO genes identified in the 2.57 Mb overlapped region on BTA20 between HO and JE have been associated with reproductive and developmental traits. The genomic region harboring these genes has been previously associated with age at puberty and calving difficulty in beef cattle (Fernández et al., 2014;Purfield et al., 2020). Furthermore, Hayes et al. (2021) found the ISL1 gene to be significantly associated with maturity in cattle. Previous studies have reported low and positive genetic correlations (r g = 0.22-0.28) between age at reproductive milestones (puberty, first calving) and calving difficulty (Bennett and Gregory, 2001;Berry and Evans, 2014). We speculate that strong selective pressure to decrease age at puberty, age at first calving, calving difficulty, or related traits has led to high homozygosity in this genomic region in HO and JE. Breed-specific selection indices offered to Holstein and Jersey breeders have weights on fertility and calving traits that range from 2 to 12% and 1 to Traits enriched in GU were mostly split between the milk or meat production trait categories. A possible explanation for this finding could be a pleiotropic effect of genes located in the genomic regions targeted by selection in GU (Gutiérrez-Gil et al., 2015). For example, the single ROH island on  Mb) identified in GU partially overlaps with a selection signature observed in multiple dairy breeds (Stella et al., 2010), and which also contains SNP with significant effects on fatty acid composition in Japanese Black cattle (Uemoto et al., 2011).

Comparison of Inbreeding Levels Among Dairy Breeds.
We calculated 3 different measures of inbreeding for animals born between 2000 to 2020, a pedigree-based measure (F PED ), a genomic markerbased measure (F GRM ), and a genomic segment-based measure (F ROH ). In all breeds, the magnitude of average inbreeding was largest for F GRM . The larger magnitude of this coefficient compared with F PED and F ROH can be explained because inbreeding calculated using a genomic relationship matrix tracks allelic homozygosity without distinguishing if the alleles are identical by descent (actual inbreeding) or not. The positive and high correlations between F PED and genomic measures and the positive and very high correlation between F GRM and F ROH were in the range of previous studies in North American dairy cattle Forutan et al., 2018). VanRaden et al. (2011) reported correlations between F PED and F GRM (using 0.5 allele frequencies) in 3 breeds (HO, JE, and BS) that ranged from 0.59 to 0.68. Forutan et al. (2018) reported a correlation between F PED and F ROH (calculated with a window size of 20 SNP) of 0.70, a correlation between F GRM and F PED of 0.64, and a correlation between F ROH and F GRM of 0.94, in Holsteins.
We observed little consistency in the ranking of breeds from least to most inbred (in terms of average inbreeding) across measures, with the sole exception of JE. The JE breed had the highest average F ROH and was tied with HO (AY) for highest average F PED (F GRM ). A previous comparison  of North American HO, JE, and BS found the highest withinbreed means for F PED and F GRM in JE. In Canadian JE, (Makanjuola et al., 2020) found consistently higher genomic inbreeding in these breeds compared with HO and similar levels of pedigree inbreeding in both breeds. Yearly averages of genomic inbreeding indicate that the inbreeding observed in JE originated before the year 2000 and has remained relatively consistent in recent years, the latter being a possible result of management strategies in this breed to curb further inbreeding accumulation.
Historical and Recent Changes in Genetic Diversity. The rate of inbreeding of a population is a valuable metric to understand the effect of management decisions and monitor genetic diversity changes across time. This study calculated the yearly rates of pedigree and genomic inbreeding in sires and dams during 3 periods known for being distinct in both the strategy and intensity of selection.
The small population of genotyped AY and GU individuals in this study prevented us from obtaining precise estimates of the yearly rate of inbreeding during the 3 periods of interest. Most estimates of the inbreeding rate in these populations had large confidence intervals that overlapped with zero. However, we did find noticeably large point estimates in the rate of genomic inbreeding in AY and GU sires in the most recent period (P3) that suggest a decrease in genetic diversity after the adoption of GS. In addition, the determination of the effects of GS in these 2 breeds is further complicated by the more recent implementation of this technology, in 2013 for AY (Wiggans and Cooper, 2013) and 2016 for GU (Cooper et al., 2016). We understand that as genotyping rates increase, future estimates of genetic diversity for these numerically small breeds will benefit from increased precision and therefore inform selection and mating strategies.
Inbreeding trends in the 3 most numerous North American breeds (BS, HO, and JE) showed a marked increase in the yearly rate of inbreeding after implementing and adopting GS. Jersey sires and dams were the only populations studied that went from a decrease to an annual increase in the rate of inbreeding after the implementation of GS. Previous studies have uncovered a similar effect of GS on the genetic diversity of JE populations (Makanjuola et al., 2020;Scott et al., 2021). Makanjuola et al. (2020) reported a shift from negative to positive rates of inbreeding after the introduction of GS in Canadian JE. Scott and colleagues (2021) reported an increase in the rate of inbreeding of Australian JE bulls after GS that ranged from a 0.75fold to 51.59-fold increase. Results for the HO breed showed the most apparent distinction in the rate of inbreeding between periods and the 2 subpopulations studied (sires and dams). First, we observed larger rates of inbreeding in dams than sires in P1, followed by a reversal of this trend during P2 and P3. The higher rates of inbreeding dams in P1 are most likely due to the genotyping strategy used. Second, the significant differences observed between the rates of inbreeding during P2 and P3 indicate an escalation in selective pressure exerted in HO after 2015. Yearly rates of genomic inbreeding during P3 exceed values of 0.5% in both sexes, reaching a maximal value of 1.50% using ∆F yr GRM in sires. A previous assessment of genetic diversity in Canadian Holstein cattle found yearly inbreeding rates that surpassed 0.30 and 0.50% for pedigree and genomic measures, respectively (Makanjuola et al., 2020). Additionally, studies of Dutch (Doekes et al., 2018), French (Doublet et al., 2019), Australian (Scott et al., 2021), and Italian (Ablondi et al., 2022) populations of Holstein-Friesian cattle have also revealed a loss in genetic diversity following their respective adoption of GS strategies.
The effective population size reflects the level of genetic diversity of populations. Therefore, populations with small effective sizes pose a severe challenge to conservation efforts and are vulnerable to decreased adaptability to new environments or selection criteria. The consensus for the minimum effective population size to maintain fitness in a population ranges from 50 to 100 individuals (Meuwissen, 2009). The estimates of effective population size have been found to vary between livestock species and breeds. Leroy et al. (2013) reported estimates of effective population size using the inbreeding rate of successive generations that ranged from 35 to 204 in cattle, 21 to 375 in sheep, 33 to 257 in horses, and 15 to 510 animals in dogs.
The estimates of effective population size for past generations showed a decreasing trend in the genetic diversity of the 5 US dairy cattle breeds. This result was expected, as the period studied, from 100 generations ago (circa 1520) to 5 generations ago (circa 1995), covers the events of breed formation, importation into the United States, creation of US breed associations and specific breed standards, and improvements in assisted reproductive technologies. Our results generally agreed with the known history of these US dairy cattle breeds. For BS and HO, we observed a sizable population bottleneck 30 generations ago (circa 1870), which corresponds to the approximate time in which the herd books of these 2 breeds were established in the United States, 1880 for BS (Brown Swiss Association, 2022b) and 1872 for HO (Holstein Association USA, 2022b). The dynamics of N e across time for AY can be divided into 3 notable periods, a period of very high and relatively stable N e estimates from about 100 to 85 generations ago (circa 1520-1595), a period of a slow decline in N e from about 85 to 43 generations ago (circa 1595-1805), and a period of seemingly faster decline in N e from 43 to 5 generations ago (circa 1805-1995). We speculate that the development of the original AY breed in Scotland (before 1800) and the later importation to and development of this breed in the United States in the early 1800s, could be responsible for the loss of genetic diversity seen in the second and third periods, respectively (U.S. Ayrshire Breeders' Association, 2022b). The estimates of N e for GU and JE were similar across most of the period studied; this was expected as the breeds have a shared history of development in the Channel Islands. Noticeably, our results showed a lower genetic diversity for these 2 breeds 100 generations ago (circa 1520) compared with the other 3 breeds studied; this could be indicative of a prior bottleneck event, possibly during the importation of the founder cattle to the islands where the original breeds would later develop.
In addition, we used estimates of inbreeding rate and generation interval to estimate N e of sires and dams in recent periods. In most populations, estimates of N e after GS were low (<100 individual). For BS and HO sires, these estimates were noticeably low with none surpassing 25 individuals. Our estimates of the most recent N e are noticeably lower than those recently reported in Canadian (Makanjuola et al., 2020; N e = 43-66 individuals) and Dutch (Doekes et al., 2018; N e = 69-79 individuals) HO populations. We speculate that methodological differences between studies are the main reason for differences in the magnitude of estimates, as we opted for estimating N e within discrete periods rather than obtaining a single estimate across periods. The only population to not exhibit concerning levels of N e in the years following GS were JE dams (N e = 99-276). The effective population sizes estimated in this study suggest an urgent need to implement strategies to better evaluate and manage the genetic diversity of US dairy cattle.

CONCLUSIONS
We used pedigree and genomic data to analyze genome-wide and region-specific patterns of inbreeding in US dairy cattle. We discovered that the distribution of autozygosity varied considerably between breeds, highlighting historical and present differences in selective pressures. Highly homozygous genomic regions suggested common breeding objectives such as milk production, health, and productive life. We found differences between breeds in the average level of inbreeding, especially for the genomic inbreeding coefficients, and differences in the rate of inbreeding accumulation at 3 different periods. We confirmed the effect of GS on genetic diversity, with HO providing the most compelling evidence. Together, our findings imply substantial heterogeneity between breeds in terms of the observable patterns of genomic homozygosity and trends in inbreeding accumulation. These results should moti-vate further research to explore strategies that could be implemented to manage homozygosity within each breed and across the entire US dairy population.

ACKNOWLEDGMENTS
Christian Maltecca, Francesco Tiezzi, and Emmanuel A. Lozada-Soto acknowledge funding from Holstein Association USA (Brattleboro, VT), American Jersey Cattle Association (Reynoldsburg, OH), and Select Sires Inc. (Plain City, OH). Paul M. VanRaden was supported by the appropriated project 8042-31000-002-00-D, "Improving Dairy Animals by Increasing Accuracy of Genomic Prediction, Evaluating New Traits, and Redefining Selection Goals" of the USDA Agricultural Research Service. Mention of trade names or commercial products in this article is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the USDA. The USDA is an equal opportunity provider and employer. We thank the Council on Dairy Cattle Breeding (CDCB, Bowie, MD) for providing the data in this study. The authors have not stated any conflicts of interest.