Genetic background of hematological parameters in Holstein cattle based on genome-wide association and RNA sequencing analyses

Hematological parameters refer to the assessment of changes in the number and distribution of blood cells, including leukocytes ( LES ), erythrocytes ( ERS ), and platelets ( PLS ), which are essential for the early diagnosis of hematological system disorders and other systemic diseases in livestock. In this context, the primary objectives of this study were to investigate the genomic background of 19 hematological parameters in Holstein cattle, focusing on LES, ERS, and PLS blood components. Genetic and phenotypic (co)variances of hematological parameters were calculated based on the Average Information Restricted Maximum Likelihood (AIREML) method and 1,610 genotyped individuals and 5,499 hematological parameter records from 4,543 cows. Furthermore, we assessed the genetic relationship between these hematological parameters and other economically important traits in dairy cattle breeding programs. We also carried out genome-wide association studies and candidate gene analyses. Blood samples from 21 primiparous cows were used to identify candidate genes further through RNA sequencing ( RNA-seq ) analyses. Hematological parameters generally exhibited low-to-moderate heritabilities ranging from 0.01 to 0.29, with genetic correlations between them ranging from −0.88 ± 0.09 (between mononuclear cell ratio and lymphocyte cell ratio) to 0.99 ± 0.01 (between white blood cell count and granulocyte cell count). Furthermore, low to moderate approximate genetic correlations between hematological parameters with one longevity, 4 fertility, and 5 health traits were observed. One-hundred-and-99 significant single nucleotide polymorphisms ( SNP ) located primarily on the


INTRODUCTION
Blood is a body fluid in the circulatory system of vertebrates that transports essential substances such as nutrients and oxygen to the cells.Hematology refers to the study of the cellular elements of the blood -red cells (erythrocytes), white cells (leucocytes), and platelets (thrombocytes).The variability in these blood parameters can be used in the diagnosing and monitoring animals' physiological states and diseases (Salvagno et al., 2014;Saravanan et al., 2020).Erythrocytes transport essential substances between different tissues via intracellular hemoglobin to maintain acid-base balance (Barbour et al., 2008).Leukocytes are the first defense against microbial invasion and contribute critically to overall immune functions (Roland et al., 2014).In addition, thrombocytes are involved in blood coagulation, wound healing, and inflammatory responses (Hoffman, Genetic background of hematological parameters in Holstein cattle based on genome-wide association and RNA sequencing analyses Tongtong Yang, 1 Hanpeng Luo, 2 Wenqi Lou, 1 Yao Chang, 1 Luiz F. Brito, 3 Hailiang Zhang, 1 Longgang Ma, 1 Lirong Hu, 1,3 Ao Wang, 1 Shanshan Li, 1 Gang Guo, 4 and Yachun Wang 1 * 2018).Hematological parameters (HP) are relevant for diagnosing disorders of the hematologic system and many other organs and systemic diseases.
Genetic regulation (Ganesh et al., 2009) and environmental factors jointly influence the variations of HPs in multiple species, which are involved with counts, ratios, volume, and distribution characteristics of peripheral blood cells and their biological activities (Kelada et al., 2012).In human populations, HPs have been extensively studied and shown to be associated with age, gender, obesity level, and smoking (Fisch et al., 1975;Okada et al., 2012;Whitfield et al., 1985).HPs can reflect the individual's health status and provide valuable insights into health conditions.The heritability of HPs in the erythrocytes-related series (ERS) in the blood routine examination ranges from 0.4 to 0.9 in human twin studies (Lin et al., 2007), indicating that variability in HPs is under genetic control.
The recent wide availability of automated hematology analyzers has facilitated the batch determination of HPs in livestock, making the process more efficient and cost-effective.Veterinary clinicians rely on hematology reference intervals (RI) to make decisions on whether laboratory results indicate pathologic changes in their patients.Initially, a study comparing RI over 50 years (1957 to 2006) in healthy North American Holstein cattle suggested that changes in husbandry and parasite control programs could decrease eosinophils (George et al., 2010).Further research has established RI for HPs in Holstein cattle in various countries, demonstrating their association with factors such as altitude, temperature, sex, pregnancy status, parity, and heat stress (Contiero et al., 2018;Cozzi et al., 2011;Jonsson et al., 2013;Kim et al., 2021;Mohri et al., 2007;Moretti et al., 2017;Moretti et al., 2015).Moreover, HPs are valuable in detecting and managing diseases in cows, including respiratory conditions (Akter et al., 2022), clinical anemia or hemorrhage (Latimer, 2011), and chronic inflammation (Zadnik, 2003).Additionally, HPs can predict disease progression and overall prognosis in individual animals (Lin et al., 2021).However, limited information exists on the genetic parameters of HPs in livestock, particularly in dairy cows, with some commonly used RI originating from research conducted in the 1960s (Roland et al., 2014).Enhancing early disease risk detection is crucial for health monitoring in dairy cattle.There is a need to investigate the genetic background of HPs in healthy Holstein cows during lactation and their genetic relationship with other relevant traits in Holstein cattle breeding programs.Therefore, the main objectives of this study were to: 1) assess the RI and the effects of environmental and physiological factors on HPs in Chinese Holstein cattle, 2) estimate variance components and genetic parameters of HPs in lactating Holstein cows, 3) perform genome-wide association studies (GWAS) for various HPs; and, 4) investigate the biological processes shared by the candidate genes identified based on GWAS and RNA sequencing (RNA-seq) analyses.

Ethics Committee Approval
The blood samples used for genotyping and RNA sequencing were collected along with the regular quarantine inspection of the farms and breeding stations.Therefore, no ethical approval was required for this study.
Data Editing.The records with no basic cow information (e.g., fertility records) and exceeding 3 standard deviations from the mean for each trait were removed from further analyses.Finally, 19 HPs recorded on 4,543 cows were kept, with the number of records per trait ranging from 5,387 to 5,499.The descriptive statistics of the systematic effects for subsequent analyses are shown in Supplemental Table S1 (https: / / doi .org/ 10 .6084/m9 .figshare.24425602).To account for the physiological status of dairy cows and ensure a balanced data distribution, the herd and testing seasons were concatenated, resulting in 10 levels.Moreover, lactation days were categorized into 3 stages (DIMs1, DIMs2, DIMs3), which encompassed 1-100 d (DIMs1), 101-200 d (DIMs2), and ≥ 201 d (DIMs3), respectively.Parity was further divided into 3 levels: 1 (Parity1), 2 (Parity2), and ≥ 3 (Parity3).
Pedigree and Genotypic Data.The pedigree data sets were provided by the Dairy Association of China (Beijing, China) and contained 20,340 animals, including 17,997 females and 2,343 males born from October 1947 to May 2020 and spanning over 3 generations.A total of 880 cows with HPs records and 1,610 animals (1,573 cows and 37 bulls) with de-regressed breeding values (dEBV) were genotyped using the Illumina 150K Bovine BeadChip (Illumina Inc., San Diego, CA, USA), which contained 123,269 single nucleotide polymorphisms (SNPs).Animals with genotype call rate lower than 0.95 (n = 2) were removed before genotype imputation.Imputation of missing SNPs was performed using the Beagle5.1 software (Browning et al., 2018).SNPs with minor allele frequency (MAF) greater than 0.05, no significant extreme deviation from Hardy-Weinberg equilibrium (P-value greater than 10 −6 ), with known chromosomes and genome position, and located on the autosomal chromosomes were kept for further analyses.After quality control, 116,509 SNPs and 1,610 animals with dEBVs were retained for the GWAS analyses.
Calculation of Reference Intervals.Reference intervals for 19 HPs in dairy cows were computed with the Reference Value Advisor freeware (RefValAdv; v.2.1), a set of macroinstructions for Microsoft Excel (Geffre et al., 2011), according to the guidelines recommended by the American Society for Veterinary Clinical Pathology (Friedrichs et al., 2012).RefValAdv output displays descriptive statistics, histograms, and Q-Q plots for each variable and calculates 95% RI based on normality and symmetry of data distribution, outliers, and sample size.
Estimation of Genetic Parameters.The GLM procedure of the SAS software (version 9.4; SAS Institute Cary, NC, USA) was utilized to determine the fixed effects to be incorporated in the genetic evaluation models for the 19 HPs traits.The effect of various lactation stages and parities on cows' HPs are presented in Figure 1, showing the least squares means and standard errors.The phenotypic and genetic (co)variances of 19 HPs were estimated based on the Average Information Restricted Maximum Likelihood (AI-REML) method with single and bivariate linear animal models implemented using the BLUPF90+ family programs (Misztal et al., 2018).The model used to estimate (co)variance components can be described as: where y is a vector of phenotypic records for 19 HPs; b is a vector of systematic effects including herd-testing season (10 levels), parity (3 levels), and lactation stages (3 levels) as previously described; X is the design matrix that associates b with y; a is a vector of random additive genetic effect; a ~N 0 2 , , Aσ a ( ) σ a 2 is the additive genetic variance, A is the pedigree-based relationship matrix; Z is the corresponding incidence matrix; and, e is the vector of random residual effects, e ~N 0 2 , , Iσ a ( ) I is an identity matrix, and σ e 2 is the residual variance.The heritability of each trait was estimated using The phenotypic and genetic correlation coefficients were calculated as r p = cov(p 1 ,p 2 )/ ˆσ σ p p 1 2 2 2 and r a = cov(a 1 ,a 2 )/ ˆˆ, σ σ a a 1 2 2 2 respectively.The σp1 2 and σp2 2 are the phenotypic variances of traits 1 and 2 and σa1 2 and σa2 2 are the additive genetic variances of traits 1 and 2, respectively.cov(p 1 ,p 2 ) and cov(a 1 ,a 2 ) are the phenotypic and additive genetic covariances, respectively.The square of the SE for the heritability estimates was calculated as (Su et al., 2007): The square of the SE for the genetic correlation coefficient was calculated as (Su et al., 2007): The variances and covariances of the estimated parameters AE, AE σ σ i ij 2 were obtained from the inverse of the average information matrix.The reliability of the estimated breeding values (EBV) was calculated as: Reli- where SEP is the standard error of prediction and σa 2 is the direct additive genetic variance for each trait.dEBVs were calculated based on the method proposed by VanRaden (VanRaden et al., 2009) and on the BLUP solutions of single-trait animal models.
Calculation of Approximate Genetic Correlations.The approximate genetic relationship between HPs and functional traits is of great importance for designing dairy cattle breeding programs.The traits included in the analyses are LONG, AFS, ICF, SB_H, SB_C, MAST, KETO, RETP, METR, and DSAB.Due  to limited access to raw data sets, the following formula was used to calculate the approximate genetic correlation between HPs with longevity, fertility, and health traits using the method proposed by Calo et al. (1973).Approximate genetic correlations were calculated based on the EBV of individuals with reliability greater than 0.10 for both traits: where ˆ, r g1 2 is the approximate genetic correlation between traits 1 and 2; ΣRL 1 and ΣRL 2 are the sum of reliabilities of traits 1 and 2; RL 1 and RL 2 are the reliabilities of traits 1 and 2; r 1,2 is the Person correlation between EBV of traits 1 and 2. The SE of the approximate genetic correlations were estimated as: Where n is the number of animals with EBVs.The EBV for all traits included in the analyses were provided by the Beijing Dairy Cattle Center (Beijing, China) in July 2022.Single-trait and 2-trait linear animal models with AI-REML were used to estimate genetic parameters and EBV for health traits (Kai et al., 2022).Single-trait linear animal models based on the AI-REML method were used to estimate genetic parameters and EBV for fertility traits (Liu et al., 2017).The trait definition and model used for the longevity trait is described in Zhang et al. (2021).Genome-Wide Association Studies.The mlma option in the GCTA software package (Yang et al., 2011) was used to estimate SNP effects based on mixed linear models with a random polygenic effect.The statistical model for single marker regression analysis can be described as: where y is a vector of dEBVs of HPs; μ is the overall mean effect; g is the vector of SNP effect; u is the vector of random polygenic effects with the variancestructure of u ~N 0 ) where G is the genomic re- lationship matrix calculated based on high-density SNP genotypes, and σ u 2 is the estimated polygenic variance; e is the vector of random residual effects with e ~N 0 2 , , Iσ e ( ) I is an identity matrix, σ e 2 is the residual vari- ance; W and Z are incidence matrices for g and u, respectively.Three traits including WLCC, RDWCV, and RDWSD were removed from further analyses due to the low reliability of their dEBVs.Thus, GWAS were done for 16 HPs.
Multiple Testing Correction.As a large number of statistical tests were performed, the traditional Bonferroni correction would be highly conservative.To avoid excessive false-negative results, we applied a modified Bonferroni correction by using the number of independent chromosome segments (M e ) at the chromosomewide level (Li et al., 2015), instead of the total number of tests.M e was calculated as (Goddard et al., 2011): Where M e is a function of both effective population size (N e ) and the average length of a chromosome in Morgan (L).The N e used was 33, which it is the most conservative value reported for the same Holstein cattle population (Shi et al., 2020).The SNP effects were considered to be statistically significant if their −log10 (P-value) was higher than the chromosome-wide threshold, which was calculated by dividing 0.05 by M e .The corresponding significance thresholds for different chromosomes ranged from 4.77 (BTA25) to 5.29 (BTA1).The values of M e and significance threshold for all chromosomes are presented in the Supplemental Table S2 (https: / / doi .org/ 10 .6084/m9 .figshare.24807780).A genomewide M e by summing across each chromosome were obtained in the same method and the corresponding significance thresholds was 6.42.Candidate Genes Detection, Quantitative Trait Loci, and Functional Enrichment Analyses.Positional candidate genes were investigated for 200 kb windows (considering the ARS-UCD1.2assembly as the reference genome) on each side of the significant SNPs using the Ensembl R package "Biomart" (Durinck et al., 2009).We scanned the corresponding QTL regions within the 200 kb genomic region surrounding the significant SNPs based on the Animal QTL database (www .animalgenome.org/cgi -bin/ QTLdb/ index) (Hu et al., 2021;Hu et al., 2007).To better understand the biological processes and pathways shared by these candidate genes, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms were enriched based on genes belonging to different blood series via the "WebGestalt" tool (Wang et al., 2017) with a cutoff P-value of 0.05.
Library Preparation and RNA Sequencing.A total of 21 primiparous cows from one herd, with DIM ranging from 117 to 166 d, were selected out of the genotyped cows with HPs measurements simultaneously for RNA-seq.The blood sample collection and measurement flow was presented in the Supplemental Figure S1 (https: / / doi .org/ 10 .6084/m9 .figshare.24901719).RNA was isolated from leukocytes according to the manufacturer's instructions for the TRIzol Reagent method (Rio et al., 2010).RNA concentration and quality were determined using Equalbit RNA BR Assay Kit (Cat No. Q10211, Invitrogen, Carlsbad, CA, USA) and Nanodrop 2000 (Thermo, Massachusetts, USA).RNA integrity was assessed using 1% agarose gel electrophoresis and used for library construction with 28S/18S > 1, A260/280 > 1.8, the concentration > 200 ng/μL, and clear electrophoretic bands without obvious degradation.We will resend a new sample if the sample does not meet the above conditions.Highquality sequencing data were obtained.Agarose electrophoresis of partial total RNA samples are presented in the Supplemental Figure S2 (https: / / doi .org/ 10 .6084/m9 .figshare.24901722).For the RNA-Seq library, 2 μg of RNA was used for purification and fragmentation using NEBNext Poly(A) mRNA Magnetic Isolation Module (Cat No. E7490S, New England Biolabs Ltd., Hitchin, Herts, UK) then building cDNA library with NEBNext Ultra RNA Library Prep Kit for Illumina (Cat No. E7530S, New England Biolabs Ltd., Hitchin, Herts, UK).All libraries were quantitated based on the Equalbit DNA BR Assay Kit (Invitrogen, CA, USA) and pooled equimolarly.All blood samples were submitted for sequencing using the NovaSeq 6000 System (Illumina Inc., San Diego, CA, USA), which generated 150 bp paired-end reads.
Quality Control and Reads Alignment.The FastQC software (v0.11.9) was used to evaluate the quality of the sequencing reads.The Fastp software (v0.20.0)default parameters (Chen et al., 2018) were used to filter all reads with the adapter, less length, low complexity, low quality, and global trimming.Subsequently, clean reads were mapped to the bovine reference genome ARS-UCD1.2using the STAR (v2.7.5c) software (Dobin et al., 2013).Gene counts were quantified using the featureCounts package (Liao et al., 2014).
Differential Expression Gene Analyses.Further differential expression gene analyses were performed for 16 HPs with significant SNPs based on the GWAS.The cows allocated to the HIGH or LOW group (3 samples vs. 3 samples) should exhibit a maximum difference in population structure and values of HPs.A principal component analysis (PCA) of read counts by RNA-seq of all cows was performed.Based on the first principal component, there is a significant difference between the HIGH group and the LOW group, as presented in Supplemental Figure S3 (https: / / doi .org/ 10 .6084/m9 .figshare.24901728).Cows in the HIGH and LOW groups exhibited extremely high and low HPs, respectively.The differences in the mean for each HP between the HIGH and the LOW groups were greater than one standard deviation of HPs for the population (Table 2).The specific values of HPs for cows in the 2 groups are presented in Supplemental Table S3 (https: / / doi .org/ 10 .6084/m9 .figshare.24807630).For screening the differentially expressed genes (DEGs), quantile-adjusted conditional maximum likelihood (qCML) analyses were performed using the DEseq2 package in R with a significance threshold of adjusted P-value <0.05 and the absolute value of fold change (|FC|) > 2.

Descriptive Statistics
The number of records, mean, standard deviation (SD), coefficient of variation (CV), and 95% RIs for HPs of the Holstein population without any clinical symptoms after quality control are shown in Table 2.In LES, WSCC accounted for more than half of WBC (53.85%), followed by WLCC (36.17%) and WMCC (9.88%).Although the CV for the ratios (WLCR, WMCR, and WCSR) were all smaller than the corresponding counts (WLCC, WMCC, and WSCC), the large differences in WLCR and WSCR's CVs were as high as 28.70% and 27.00%, respectively.Combining the 3 series, the mean of RBC far exceeded WBC and PLT, while the CV for them was the opposite case.For both volume HPs, the MCV in ERS was 48.53 fL, which is approximately 8 times the mean value of MPV in PLS.In addition, the 3 HPs with the widest RIs were PLT (240.00 -701.00),MCHC (188.89-385.60), and HGB (66.00 -121.00), while the 3 HPs with the narrowest RIs were WMCR (0.40 -1.90), MPV (5.70 -9.00), and RBC (4.58 -7.91).

Yang et al.: GENETICS OF HEMATOLOGICAL PARAMETERS IN HOLSTEINS
Statistically significant (P < 0.05) phenotypic and genetic correlations are presented above and below the diagonal of Figure 2, respectively.All phenotypic correlations were positive, and stronger phenotypic correlations were observed between HPs within the same blood series.Within ERS, high phenotypic correlations were found between 5 HPs, including 0.98 ± 0.01 (MCH and MCHC), 0.96 ± 0.01 (RBC and HCT), and 0.76 ± 0.01 (HCT and MCH).In contrast, the phenotypic correlations between the 3 blood series ranged from 0.03 ± 0.01 (RBC and PLCR) to 0.14 ± 0.02 (WMCR and PLCR).The genetic correlations were markedly higher than phenotypic correlations and followed the same trend with phenotypic correlations being stronger between HPs within the same HP groups.However, the direction of genetic correlation between some HPs changed.There were 19 pairs of positive genetic correlations for each HPs within 3 blood series and the highest was 0.99 ± 0.01 (WBC and WSCC).Among the 3 blood series, we observed 14 pairs of positive genetic correlations, such as 0.99 ± 0.22 (WBC and PDW) and 0.92 ± 0.32 (HGB and PDW).

Genetic and Approximate Genetic Correlations Between Hematological Parameters and Other Relevant Traits
Milk traits, including MY, PP, PF, and LP, exhibited high negative genetic correlations with 5 hematological parameters in the LES and ERS, ranging from −0.34 to −0.73.Only PF and WLCR (0.54 ± 0.28) and LP and PLCR (0.41 ± 0.19) displayed positive genetic correlations.The approximate genetic correlations between HPs and other economically important traits are shown in Table 4.The approximate genetic correlations between LONG and all 5 HPs within the ERS were positive and ranged from 0.38 ± 0.03 (HCT) to 0.71 ± 0.04 (MCH), while 3 HPs within PLS were negatively correlated with LONG, with correlations ranging from −0.40 ± 0.09 (LONG and MPV) to −0.83 ± 0.01 (LONG and PDW).Moreover, the estimated approximate genetic correlations between LONG and all HPs with the LES were negative and ranged from −0.13 ± 0.04 (WSCR) to −0.56 ± 0.01 (WMCC), except for a positive approximate genetic correlation between LONG and WLCR (0.48 ± 0.02).For fertility traits, ICF had high approximate genetic correlations with all HPs ranging from −0.86 ± 0.010 (HCT) to −0.67 ± 0.01 (WLCR) (negative genetic correlations) and 0.55 ± 0.01 (PDW)  to 0.92 ± 0.01 (WSCC) (positive genetic correlations).Notably, the approximate genetic correlations between AFS and all 5 HPs within the ERS group were moderate to high and ranged from 0.30 ± 0.02 (HCT and AFS) to 0.98 ± 0.01 (MCH and AFS).Moreover, there were 7 pairs of negative correlations and only 1 pair of positive approximate genetic correlations between HPs and SB_H and SB_C.In total, 21 pairs of approximate genetic correlations between the 5 health traits and all 5 HPs from the ERS group were positive and ranged from 0.34 ± 0.02 (MAST and HCT) to 0.86 ± 0.01 (RETP and HCT).There were also 37 pairs of negative approximate genetic correlations between all 5 health traits and 8 HPs within LES and PLS.However, the approximate genetic correlations of all 5 HPs within the LES with LONG, fertility, and health traits were in the opposite direction to WLCR, as this trait was negatively genetically correlated with all other HPs from the LES group.

GWAS
GWAS were performed for each of the 16 HPs using the GCTA software.A total of 199 SNPs reached the chromosome-wide threshold in association with 16 HPs.Out of these 199 SNPs, 26 SNPs reached the genomewide threshold and were associated with 10 HPs, as shown in Figure 3.Additional details about the GWAs results are presented in the Supplemental Table S4 (https: / / doi .org/ 10 .6084/m9 .figshare.25025087).The Manhattan plots for 3 individual traits are presented in the Supplemental Figure S4 (https: / / doi .org/ 10 .6084/m9 .figshare.24901731).The lambda values ranged from 0.94 to 1.24.The number of significantly associated SNPs for each HP ranged from 2 (PLCR) to 29 (HCT), with an average of 12.44 SNPs per HP.Regarding chromosomal distribution, BTA8 exhibited the highest number of SNPs significantly associated with HPs (25 SNPs), followed by BTA6 (21 SNPs) and BTA3 (16 SNPs).Notably, the 199 significant SNPs showed no overlap, indicating that HPs are highly polygenic and had distinct and uniquely significant SNPs.Further search within a flanking region of 200 kb on either side of the significant SNPs identified 420 protein-coding candidate genes associated with 16 HPs across 27 chromosomes, as presented in Supplemental Table S5 (https: / / doi .org/ 10 .6084/m9 .figshare.24425758).Overall, as shown in Table 5, the most significant SNPs associated with HPs and overlapping with candidate genes were Hapmap51971-BTA-18711 (WSCR, 10 genes), ARS-BFGL-NGS-17861 (HGB, 9 genes), and ARS-BFGL-NGS-24347 (PLT, 8 genes).Among the 199 significant SNPs identified, 75 were located within genes, and 3 belonged to the significant SNPs mentioned above.More specifically, BTA8, BTA19, and BTA5 had the highest number of candidate genes, with 52, 40, and 39 genes, respectively.Notably, the BTA1 region stood out with 9 unique genes belonging to the Olfactory Receptor Family 5.

Functional Enrichment Analyses
Table 6 presents an overview of the top KEGG pathways and GO terms based on the protein-coding candidate genes from LES, ERS, and PLS.The KEGG pathways enriched for these candidate genes include the Estrogen Signaling pathway as the most significantly enhanced pathway (P = 3.28 × 10 −9 ).In addition, other enriched pathways include the Chemokine Signaling pathway, Inflammatory Mediator Regulation of TRP Channels, Cytokine-cytokine Receptor Interaction, and Gastric Acid Secretion.Furthermore, the most significant GO terms in each of the 3 blood series are Regulation of Cardiac Muscle Contraction (P =   3.61 × 10 −3 ), Cellular Response to Interleukin-6 (P = 4.40 × 10 −3 ), and Positive Regulation of Transferase Activity (P = 3.71 × 10 −4 ), respectively.

Identification of Differentially Expressed Genes
The DEGs were detected based on the qCML.As illustrated in Figure 4, 10,806 DEGs were identified between the HIGH and LOW groups across 16 HPs, in which 6,687 of them were downregulated and 4,119 ).When compared with the HIGH group, the LOW group exhibited upregulation of only 3 genes, namely GPAT3, WDFY3, and RTL6, while downregulation was observed in the remaining 10 genes.More information on these genes is presented in Supplemental Table S7 (https: / / doi .org/ 10 .6084/m9 .figshare.25025090).The expression of 13 candidate genes in high and low groups was shown in Supplemental Figure S5 (https: / / doi .org/ 10 .6084/m9 .figshare.24901734).The investigation of corresponding QTL regions was carried out on the 13 SNPs that identified overlapping candidate genes above in both data sets.Supplemental Table S8 (https: / / doi .org/ 10 .6084/m9 .figshare.25025096)presents the main and total QTL numbers.

Descriptive Statistics and Reference Intervals of Hematological Parameters
Improving herd health programs depends on laboratory blood analyses, which is an essential tool for dairy veterinarians to assess the health status of cows (Brscic et al., 2015).The most appropriate RIs are from healthy animals with environmental and physiological factors similar to those of the assessed individuals.Chen et al. (2019) measured HPs in Holstein cattle in Beijing during summer periods, and some of the HPs, such as PLCR and MPV, showed slightly different values compared with the present study.Similarly, HPs of dairy cows from tropical herds were higher for 6 parameters in the ERS than those in the present study, including RBC, HGB, HCT, MCV, MCH, and MCHC (Vallejo-Timarán et al., 2020).Moreover, the 95% RIs for HPs from previously published studies from different Holstein populations are presented in Supplemental Table S9 (https: / / doi .org/ 10 .6084/m9 .figshare.25025108),including RIs for neonatal calves (Panousis et al., 2018), bulls (Contiero et al., 2018), cows between 30 and 120 d in milk, and cows at 2 d after calving (Tsiamadis et al., 2022).The mean values of the HPs in this study were all within 95% RIs from previous studies, while the MCHC in the studied population was slightly below the 95% RI.Like all species, Holstein cattle exhibit some variability in their HPs.Various factors such as age, sex, stressors, diet, body condition, activity level, hydration, ambient temperature, altitude, and lactation/pregnancy stage contribute to the thresholds and width of RIs, leading to variations between the results of this study and those previously reported in the literature.This indicates that it is necessary to establish RIs according to the different characteristics of the cows.However, HPs of Holstein cattle share standard features with most mammals, such as RBC significantly exceeding WBC and PLT (De Melo et al., 2019) and MCV being approximately 6-8 times the MPV values (Zhang et al., 2014).The current study, based on a large and comprehensive data set (5,499 blood samples and 4,543 cows), provides a solid foundation for disease detection and subsequent analyses of dairy cows.

Parity and Lactation Effect HPs
In modern dairy herds, managing cows to maintain their health status, fertility, and milk yield during each lactation period poses significant challenges.Additionally, the longevity of Holstein cows in large dairy herds in China remains below 3 parities (Yan et al., 2016), far from reaching the peak productivity stage for dairy cows.Additional information about individual cows, such as blood and milk analyses, can aid in managing the physiological challenges of high-producing dairy cows at both the herd and individual levels.Therefore, this study investigated the effects of different physiological stages on dairy cows' HPs.
In total, 6 HPs within ERS of primiparous cows were higher than the values from later parity cows, and the change of WMCC reached a significant level, which is similar to previous reports (Mirzadeh et al., 2010;Ziwei et al., 2019).Many studies have shown that the DNA methylation level of white blood cells in human peripheral blood decreases with age, which may accelerate the aging or death of white blood cells, leading to decreased WBC in peripheral blood.Moreover, the previous study revealed that monocytes at birth are 3 times higher than in adults (Schefold et al., 2015), a similar trend to the decrease of WMCC with increasing parity in cows from the present study.For ERS, MCV increased with parity and reached its highest value above 3 parities.Although not statistically significant, RBC, HGB,    MCH, MCHC, and RDWCV decreased to their lowest in the second parity, which is consistent with other studies (Vallejo-Timarán et al., 2020).RBC, HGB, and HCT could reflect the body's ability to transport oxygen (Patel, 2008).Therefore, the ability of dairy cows to transport oxygen might decrease in later parities.Clinically, the decreases in RBC, HGB, and HCT and the increases in MCV, MCH, and RDWSD were all critical indicators of malnutrition, anemia, and hematopoietic decline (Nutjeera et al., 2015).Lower RBC, HGB, and HCT in the blood can result in endothelial cells having hypoxia.This could result in an abnormal increase in red blood cell volume and dispersion and a compensatory increase of hemoglobin to improve oxygen-carrying capacity, which may lead to an increase in MCH.In addition, large red blood cells could be crushed and damaged when passing through the vascular system (Hall et al., 1995).In summary, the ability of dairy cows with high parity to transport oxygen is expected to decrease, and symptoms of malnutrition and anemia tend to become clearer at first parity.Five HPs in LES and ERS in DIMs1 were significantly lower than in other lactation stages, including WMCC, WSCR, WMCR, HGB, and RDWCV.Meanwhile, MCV and MCH were found to be the lowest in DIMs2.Dairy cows mobilize fat during early lactation to alleviate the negative energy balance.At the same time, the milk yield of dairy cows increases, and highintensity metabolism promotes the body to produce oxidative stress, increasing the possibility of infection with pathogenic bacteria.Delivery stress and the substantial decrease in protein intake would make dairy cows lose physical fitness, leading to anemia in dairy cows at this stage (Blanco et al., 2009).However, PLT and PLCR within PLS at DIMs3 were significantly lower than in other lactation stages, reflecting the decreased platelet coagulation activity in dairy cows at this stage.Laterparity and early lactation cows might have lower immunity, increased susceptibility to bacteria, and increased risk of anemia.Therefore, regular testing of HPs could contribute to a better understanding of the health level of dairy cows, adjust the feeding management plan as early as possible, reduce the possibility of disease in different physiological stages, and improve management practices, thereby improving the profitability of dairy herds.

Heritability and Genetic Correlations with Hematological Parameters
Estimating genetic parameters for HPs is key to better understanding the genetic background of immune and health-related traits in dairy cattle populations.Accurate estimation of these parameters is essential for effective herd management and refining breeding objectives in dairy cattle breeding programs.Our study also revealed that most HPs are heritable, with heritability estimates ranging from 0.01 to 0.29.The heritabilities of HPs in beef cattle (Angus) have been reported to range from 0.11 to 0.52 (Chinchilla-Vargas et al., 2020).Similarly, the genomic-based heritability estimates of different HPs in yaks (Bos grunniens) were found to vary from 0.23 to 0.61.Interestingly, the heritability of WBC was also found to be moderately heritable in a Sardinian human population, with estimates ranging from 0.14 to 0.40 across WBC and other HPs in the LES group (Pilia et al., 2006), which differs from our results.The differences in the estimates might be partly attributed to variations in analytical methods, the age of the animals, environmental factors, or the lack of standardization protocols across laboratories.Furthermore, the divergences could also indicate differences in the genetic determinants of immune capacity among different species.
In line with previous findings, we also observed that the genetic correlations were higher than phenotypic correlations in beef cattle (Chinchilla-Vargas et al., 2020).Moreover, both correlations followed the same trend of being stronger within trait groups.The high genetic correlations observed between HPs indicate potential pleiotropic effects (Lukowski et al., 2017).
Although direct physiological modulation between milk components and HPs has not been conclusively proven, our study revealed statistically significant and robust genetic correlations between these 2 trait groups.Research has established that maintaining HPs within reasonable ranges is key for animal health.Therefore, it is essential to consider the antagonistic relationship between production traits and HPs when breeding cows for increased milk production.The different leukocytes in the blood of healthy cows directly reflect the immune level, which could explain the moderate and negative genetic correlations between LES HPs and the health traits evaluated in this study.Moreover, decreased leukocytes could relate to peracute inflammation, metabolic disorders, liver disease, and infectious diseases (Kraft et al., 2005;Roland et al., 2014).A platelet decrease is due to blood loss, disseminated intravascular coagulation, infections, neoplasia, or immune-mediated (Kraft et al., 2005;Roland et al., 2014).

Trait-Related Genes
Despite no overlap among the SNPs identified for different HPs, they are close to each other.For example, Hapmap38825-BTA-28013 and BTB-00501626 are associated with HGB and WSCR, respectively, and are located 56 kb apart.We found 420 protein-coding genes associated with 16 HPs within 3 blood series, and these genes are involved in various biological processes, including gastric acid secretion, endocytosis, cardiac muscle contraction, T cell differentiation, neurotrophin signaling, cortisol synthesis and secretion, and apoptotic signaling.
The most significant SNP, Hapmap51971-BTA-18711, was found to be associated with WSCR and located within ITSN2 (Intersectin 2).ITSN2 is a conserved scaffold protein involved in endocytic internalization, actin cytoskeleton regulation, and epithelial morphogenesis (Novokhatska et al., 2011;Pucharcos et al., 2000;Rodriguez-Fraticelli et al., 2010).Another significant SNP, ARS-BFGL-NGS-17861, was associated with HGB and found near 9 genes, including ASXL1 (Additional Sex Combs Like Transcriptional Regulator 1), CCM2L, HCK (Hemopoietic Cell Kinase), KIF3B (Kinesin Family Member 3B), PDRG1 (P53 And DNA Damage Regulated 1), PLAGL2, POFUT1 (Protein O-Fucosyltransferase 1), TM9SF4 (Transmembrane 9 Superfamily Member 4), and XKR7 (X Kell Blood Group Precursor-Related Family, Member 7).A key role for ASXL1 in erythropoiesis has been unveiled, and the loss of ASXL1 hinders erythroid development/ maturation (Shi et al., 2016).In cattle studies, HCK is essential in synchronizing the immune response of the parenchyma against mastitis-causing bacteria (Kosciuczuk et al., 2017).In the PLS, CD163 (CD163 Molecule) and NTRK2 (Neurotrophic Receptor Tyrosine Kinase 2) were identified as candidate genes previously associated with the innate immune response.For example, CD163 is a monocyte/macrophage-specific marker expressed predominantly on cells with strong anti-inflammatory potential (Shiraishi et al., 2018), and CD163 gene-edited-pigs are resistant to porcine reproductive and respiratory syndrome (Wei et al., 2020).Confirming the results of GWAS is crucial to increase the confidence of animal breeders in genetic improvement and validate scientific hypotheses.Various methods can be employed to identify candidate genes from GWAS results, including SNP chip analysis (e.g., Zhao et al., 2019), qRT-PCR (e.g., Luo et al., 2021), or association analysis in other populations (e.g., Dikmen et al., 2015).RNA-seq, with high-throughput sequencing, provides valuable insights into translating genetic loci into biological mechanisms underlying essential traits (Zhao et al., 2021).Integrating GWAS and RNA-seq data through the overlap analysis has improved the accuracy and precision of candidate gene selection to some extent.It has been widely used in animals (e.g., Luo et al., 2022), plants (e.g., Zhao et al., 2021), and humans (e.g., Chang et al., 2022;Collado-Torres et al., 2019;Sun et al., 2019).This RNA-seq study focused on primiparous cows with similar DIM.Among the 13 protein-coding genes identified based on the GWAS (ACRBP, ADAMTS3, CANT1, CCM2L, CNN3, GPAT3, GRIP2, PLAGL2, RTL6, SOX4, WDFY3, and ZNF614), all showed significant differences in expression between cows with high and low HPs.Previous research has found associations of ADAMTS3 with dysplasia of the lymphatic system and pulmonary hypoplasia due to a missense variant (Haefliger et al., 2020), while both WDFY3 and PLAGL2 have been linked to carcass and meat quality (Grigoletto et al., 2020).Additionally, ADAMTS3 has been shown to activate vascular endothelial growth factor and promote lymphangiogenesis (Janssen et al., 2016), and WDFY3 has been identified as a novel regulator of macrophage efferocytosis (Shi et al., 2022).Interestingly, CCM2L, a paralog of CCM2, was selectively expressed in endothelial cells during active cardiovascular growth, aligning with our study's results (Zheng et al., 2012).Based on the integration of GWAS and RNAseq analyses and functional annotations, ADAMTS3, WDFY3, CCMM2L, and PLAGL2 were identified as promising candidate genes regulating HPs.This way has improved the accuracy and precision of candidate gene selection.We will continue to improve the number of animals with hematological parameters and add whole genome resequencing to explore how SNPs specifically affect gene expression.These findings contribute to a better understanding of the genetic background of HPs in Holstein cattle.

CONCLUSIONS
The HPs of Holstein cows exhibit moderate heritability, with notable high and positive genetic correlations observed among them.Furthermore, given the genetic association of HPs with certain milk performance, fertility, and health traits, these relationships could be considered when refining selection indexes in dairy cattle breeding programs.Through GWAS and RNAseq, numerous genomic regions located across various chromosomes were identified to be associated with variability in HPs.These regions were primarily enriched in the pathways of Cardiac Muscle Contraction Regulation and Cellular Response to Interleukin-6, suggesting a polygenic yet significant inheritance of HPs.Thirteen genes (ACRBP, ADAMTS3, CANT1, CCM2L, CNN3, CPLANE1, GPAT3, GRIP2, PLAGL2, RTL6, SOX4, WDFY3, and ZNF614) emerged as the main candidate genes influencing HPs in Holstein cows.In summary, our findings shed light on the genetic architecture of HPs, unveiling new avenues for utilizing HPs as indicator traits to enhance disease resistance and improve welfare in dairy cows.
Yang et al.: GENETICS OF HEMATOLOGICAL PARAMETERS IN HOLSTEINS Yang et al.: GENETICS OF HEMATOLOGICAL PARAMETERS IN HOLSTEINS Figure 1.Effect of different parities and lactation stages on Holstein cattle's hematological parameters (least squares means ± standard error).Peer data without the same letter in the shoulder marker showed a significant difference (P < 0.05) while those with the same letter showed no significant difference (P > 0.05).LES = Leukocytes series; ERS = Erythrocytes series; PLS = Platelet series; WBC = White blood cell count; WLCC = Lymphocyte cell count; WMCC = Mononuclear cell count; WSCC = Granulocyte cell count; WLCR = Lymphocyte cell ratio; WMCR = Mononuclear cell ratio; WSCR = Granulocyte cell ratio; RBC = Red blood cells; HGB = Hemoglobin; HCT = Hematocrit; MCV = Mean corpuscular volume; MCH = Mean corpuscular hemoglobin; MCHC = Mean corpuscular hemoglobin concentration; RDWSD = Red blood cell distribution width of standard deviation; RDWCV = Red blood cell distribution width of coefficient of variation; PLT = Platelet count; PLCR = Large cell ratio of platelet; MPV = Mean platelet volume; PDW = Platelet distribution width.
Figure 2. Statistically significant phenotypic (above diagonal) and genetic (below diagonal) correlations between 16 hematological parameters and 5 economically important traits.The gradient of color from brown to green represents negative to positive correlations and their strength, respectively.White parts mean that the correlations failed to reach the significance level (P > 0.05); SCS = somatic cell score; MY = milk yield; PP = protein percentage; PF = fat percentage; LP = lactose percentage; WBC = White blood cell count; WSCC = Granulocyte cell count; WMCC = Mononuclear cell count; WSCR = Granulocyte cell ratio; WMCR = Mononuclear cell ratio; WLCR = Lymphocyte cell ratio; RBC = Red blood cells; HGB = Hemoglobin; HCT = Hematocrit; MCV = Mean corpuscular volume; MCH = Mean corpuscular hemoglobin; MCHC = Mean corpuscular hemoglobin concentration; PLT = Platelet count; PDW = Platelet distribution width; MPV = Mean platelet volume; PLCR = Large cell ratio of platelet.

Table 1 .
Information of nineteen hematological parameters and five milk traits in Holstein cattle

Table 2 .
Yang et al.: GENETICS OF HEMATOLOGICAL PARAMETERS IN HOLSTEINS Descriptive statistics and reference intervals for hematological parameters in healthy Holstein cows

Table 3 .
Estimates of variance components, heritability, and standard errors for 19 hematological parameters in Chinese Holstein cattle

Table 5 .
Single nucleotide polymorphisms (SNPs) and candidate genes significantly associated with hematological parameters in Holstein cattle HP

Table 6 .
Top KEGG pathways and gene ontology (GO) terms revealed for hematological parameters in Holstein cattle Series Bolded candidate genes indicate that the genes were present in both KEGG and GO; LES = Leukocytes-related series; ERS = Erythrocytes-related series; PLS = Plateletrelated series.