Genome-Wide Association Study Based on Random Regression Model Reveals Candidate Genes Associated with Longitudinal Data in Chinese Simmental Beef Cattle

Simple Summary Genome-wide association study (GWAS) has become the main approach for detecting functional genes that affects complex traits. For growth traits, the conventional GWAS method can only deal with the single-record traits observed at specific time points, rather than the longitudinal traits measured at multiple time points. Previous studies have reported the random regression model (RRM) for longitudinal data could overcome the limitation of the traditional GWAS model. Here, we present an association analysis based on RRM (GWAS-RRM) for 808 Chinese Simmental beef cattle at four stages of age. Ultimately, 37 significant single-nucleotide polymorphisms (SNPs) and several important candidate genes were screened to be associated with the body weight. Enrichment analysis showed these genes were significantly enriched in the signaling transduction pathway and lipid metabolism. This study not only offers a further understanding of the genetic basis for growth traits in beef cattle, but also provides a robust analytics tool for longitudinal traits in various species. Abstract Body weight (BW) is an important longitudinal trait that directly described the growth gain of bovine in production. However, previous genome-wide association study (GWAS) mainly focused on the single-record traits, with less attention paid to longitudinal traits. Compared with traditional GWAS models, the association studies based on the random regression model (GWAS-RRM) have better performance in the control of the false positive rate through considering time-stage effects. In this study, the BW trait data were collected from 808 Chinese Simmental beef cattle aged 0, 6, 12, and 18 months, then we performed a GWAS-RRM to fit the time-varied SNP effect. The results showed a total of 37 significant SNPs were associated with BW. Gene functional annotation and enrichment analysis indicated FGF4, ANGPT4, PLA2G4A, and ITGA5 were promising candidate genes for BW. Moreover, these genes were significantly enriched in the signaling transduction pathway and lipid metabolism. These findings will provide prior molecular information for bovine gene-based selection, as well as facilitate the extensive application of GWAS-RRM in domestic animals.


Introduction
With the development of the high-throughput chip technologies and the completion of whole-genome sequencing of swine [1], cattle [2], sheep [3], chicken [4], and other domestic animals [5], genome-wide association study (GWAS) has become an indispensable statistical method that can detect significant single-nucleotide polymorphisms (SNPs) and study were derived from the Chinese Simmental beef cattle resource population established in Ulgai, Xilingol League, Inner Mongolia of China from 2008 to 2014. After weaning, the cattle were moved to the Beijing Jinweifuren fattening farm for fattening in the same feeding strategies and management conditions. Body weight was measured for each individual at 0, 6, 12, and 18 months after birth, respectively. Here, the body weight data were consistent with the data in the study by Duan et al. [24].

Genotyping and Quality Control
Blood samples for the experimental population were collected along with the periodic quarantine inspection of the farm. Genomic DNA was isolated from blood samples using the TIANamp Blood DNA Kit (Tiangen Biotech Co.Ltd., Beijing, China), and the highquality DNAs with the A260/280 ratio ranging 1.8-2.0 were considered for further analysis. In this study, the Illumina BovineHD Beadchip with 774,660 SNPs (Illumina Inc., San Diego, CA, USA) was used for qualified DNAs genotyping and Illumina's Infinium II Assay was selected as the genotyping platform. The SNPs were uniformly distributed on the whole bovine genome with a mean inter-marker space of 3.43 kb. SNP chips were scanned and analyzed using the Infinium GenomeStudio software (Illumina Inc., San Diego, CA, USA). PLINK v1.9 (http://zzz.bwh.harvard.edu/plink/, accessed on 1 July 2021) was used for quality control of SNPs according to the following empirical excluded criteria: (1) minor allele frequency (MAF) < 0.01; (2) SNP call rate (CR) < 95%; (3) Hardy-Weinberg equilibrium value p < 1 × 10 −6 ; (4) Mendelian error of SNP genotype above 2%; (5) Individuals with more than 10% SNPs deletion; (6) SNP marker sites with missing chromosomal location information. All the misplaced and duplicated SNPs were also excluded from the analysis. Ultimately, 671,192 SNPs with an average marker interval of 3 kb on 29 autosomal chromosomes remained for subsequent analysis.

Population Stratification
Population stratification usually caused serious FPR in GWAS analysis. Here, we performed a principal component analysis (PCA) by PLINK v1.9 [25]. Our previous work demonstrated the first two principal components had been selected as covariances to eliminate the influence of population stratification [24].

Genome-Wide Association Study Based on the Random Regression Model
The general expression of the random regression model is as follows: y tijk = F i + f (t) j + r(a, x, m 1 ) k + r(p, x, m 2 ) k + e ijkt (1) where y tijk is the measured value of individual k at time t; F i is the time-independent fixed environmental effect; f (t) j is fixed regression function, reflecting the average change trend of phenotypic values of animals in group j with time t; r(a, x, m 1 ) k and r(p, x, m 2 ) k are random regression functions, which represent time-varied additive genetic effect and permanent environmental effect for individual k, respectively; a and p are the random regression coefficients of additive genetic effect and permanent environmental effect, respectively; m 1 and m 2 are the orders of the corresponding regression function; x is a covariable; e ijkt is the time-independent random residual for each measurement of individual k at time t. Here, f (t) j , r(a, x, m 1 ) k and r(p, x, m 2 ) k can be described as the Legendre polynomial regression for a set of basis functions, specific form as follows: where b jl is the lth fixed regression coefficient; a jkl and p jkl are the lth random regression coefficients for additive genetic effect and permanent environmental effect of the kth individual, respectively; m f , m r1 , and m r2 are the orders of corresponding basis functions; The orders of different basis functions can be determined by model selection criteria proposed by Das et al. [26]. For instance, in the present study, Akaike information criterion (AIC) and Bayesian Information Criterion (BIC) were used to determine a fifth-order basis function for the population mean, a third order for additive genetic effects and a fifth order for permanent environmental effects were best fit to the data for BW trait. φ tjkl . represents the value of the lth basis function at time t, i.e., the value of the Legendre polynomial (covariable). The Equation (1) can be denoted as: We assume that there are m out of n individuals for which phenotypic values are measured; m k represents the phenotypic record for individual k and the total number of records for all individuals is m = ∑ n k=1 m k . Therefore, y is the vector (m × 1) of phenotypic values of all individuals. The parameters in the equation are defined as follows: Parameter f is a vector of fixed environmental effect and parameter X 1 is the corresponding incidence matrix; Parameter b is the vector (m f + 1) of fixed regression coefficients; Parameter a (n × (m r + 1)) and p (m × (m r + 1) ) are the vector of random regression coefficients for additive genetic effect and permanent environmental effect respectively for each individual; Parameter X 2 , Z 1 , and Z 2 are the corresponding covariance matrix; e is the vector of random residuals.
For matrix form (3), we have the (co) variance matrices of all random effects: here, A is the numerator relationship matrix based on pedigree information; D and P are the variance-covariance matrix of random regression coefficients for additive genetic effect and permanent environmental effect, respectively; I is the identity matrix; σ 2 e is the residual variance; The symbol "⊗" represents the Kronecker product. Therefore, the mixed model equations can be represented as: Based on the general random regression model, the functional GWAS model (fGWAS-C) for the association analysis of longitudinal traits has been proposed by Ning, C. [18], which adds an additional fixed regression term to Equation (1) to account for the effect of the SNP. Th fGWAS-C can be expressed as: here, x i is a genotype indicator variable that is coded as 0, 1, or 2 for the three genotypes, aa, Aa, and AA, respectively. SNP(t) represents the time-varied additive effect for each SNP at time t. The function expression of SNP(t) is: where m f is the order of basis functions for the time-varied SNP effect; η l is the lth fixed regression coefficient of SNP additive effect; φ tjkl is the value of the lth basis function at time t. The threshold p-value for GWAS analysis was calculated as follows: where FDR was usually set at 0.05; n is the number of SNPs with p < 0.05; m is the total number of SNPs after quality control [27].

Detection and Functional Enrichment of Candidate Gene
The significant SNPs associated with BW were screened according to the threshold p-value and then the Ensembl-BioMart was used to match these SNPs with the bovine reference genome UMD 3.1 (http://www.ensembl.org/Biomart, accessed on 1 July 2021). Candidate genes in the target region were screened through the National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/, accessed on 1 July 2021), cattle QTLdb [13], and previous relevant studies. Then Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed the function of candidate genes, and protein-protein interaction (PPI) network explored the interaction between node proteins encoded by these genes [28]. The p-value adjusted using the Benjamini-Hochberg approach (p-value < 0.05) was considered to be the threshold value for significantly enriched GO terms and pathways. Using ToppCluster and Cytoscape v3.8.0 to plot the network diagram between candidate genes and their belonged GO terms/pathways. The STRING (https://string-db.org/, accessed on 1 July 2021) was used to perform PPI network analysis.

Statistical Analysis
Using SPSS v20.0 to calculate the measurement values of body weight at four growth stages. All data were expressed as means ± standard deviation (M ± SD). Microsoft Excel 2010 software was used to check up the data with normality test.

Data Statistics of Body Weight
The BW values of Chinese Simmental beef cattle at 0, 6, 12, and 18 months of age were tested for normality. As shown in Figure 1, the phenotypic values of the individuals at different months of age presented normal distributions, which was in line with the hypothesis of the model and could be used for subsequent association analysis. The descriptive statistics of the BW trait at different months of age were presented in Supplementary  Table S1.

Genome-Wide Association Study Based on the Random Regression Model
The quantile-quantile (Q-Q) and Manhattan plots of GWAS-RRM analysis are shown in Figures 2 and 3, respectively. Distributions of the observed −log 10 (p) versus expected −log 10 (p) in the Q-Q plot represented most points revolving around the 45 • line, indicating there was no inflation or systematic bias in this study, as well as population stratification was well controlled. The Manhattan plots showed that a total of 37 significant SNPs associated with BW trait were identified, most of which were located on Bos taurus autosome (BAT) 1 (five SNPs), BAT 2 (three SNPs), BAT 11 (four SNPs), and BAT 12 (three SNPs). Two SNPs were found on each of the seven chromosomes ( BAT 4,5,7,8,14,16,and 29), and only one SNP was distributed on each of the eight chromosomes (BAT 3, 6, 9, 10, 13, 19, 21, and 25). The SNP with the smallest p-value (p = 7.55 × 10 −8 ) was located on BAT 10: 101,577,026 bp (BovineHD1000029459). However, some SNPs had the lowest significance levels, i.e., the p-value of BovineHD1100004962 and BovineHD1100011885 located on BAT 11 was 2.54 × 10 −6 and 3.42 × 10 −6 , respectively. The detailed information of all significant SNPs was displayed in Table 1 and Supplementary Table S2. A shows the normal distribution of body weight at 0 months of age; B shows the normal distribution of body weight at 6 months of age; C shows the normal distribution of body weight at 12 months of age; D shows the normal distribution of body weight at 18 months of age. In the figure, the abscissa represents weight, the left ordinate represents frequency, and the right ordinate represents normal function values.

Genome-Wide Association Study Based on the Random Regression Model
The quantile-quantile (Q-Q) and Manhattan plots of GWAS-RRM analysis are shown in Figure 2 and Figure 3, respectively. Distributions of the observed −log10(p) versus expected −log10(p) in the Q-Q plot represented most points revolving around the 45° line, indicating there was no inflation or systematic bias in this study, as well as population stratification was well controlled. The Manhattan plots showed that a total of 37 significant SNPs associated with BW trait were identified, most of which were located on Bos taurus autosome (BAT) 1 (five SNPs), BAT 2 (three SNPs), BAT 11 (four SNPs), and BAT 12 (three SNPs). Two SNPs were found on each of the seven chromosomes ( BAT 4,5,7,8,14,16,and 29), and only one SNP was distributed on each of the eight chromosomes (BAT 3, 6, 9, 10, 13, 19, 21, and 25). The SNP with the smallest p-value (p = 7.55 × 10 −8 ) was located on BAT 10: 101,577,026 bp (BovineHD1000029459). However, some SNPs had the lowest significance levels, i.e., the p-value of BovineHD1100004962 and BovineHD1100011885 located on BAT 11 was 2.54 × 10 −6 and 3.42 × 10 −6 , respectively. The detailed information of all significant SNPs was displayed in Table 1 and Supplementary Table S2.

Genes Detection
In this study, 37 significant SNPs were identified for BW trait, of which 33 were located within or near candidate genes through searching Ensembl, NCBI, and QTL databases. As shown in Table 1, significant SNPs on BAT5, BAT13, and BAT29 were within or close to genes ITGA5, ANGPTL4, and FGF4, respectively. Three significant SNPs were identified at the region of 26.11-33.58 Mb on BAT12, of which BvineHD1200027798 and BovineHD1200026844 were close to genes SHISA2 and CHMP3, respectively. No candidate genes were detected in four significant SNPs including BovineHD0200018897, Bo-vineHD0700012290, BovineHD1100030552, and BovineHD2900008350, located on BAT2, BAT7, BAT11, BAT12, and BAT29, respectively. The detailed information on these genes was listed in Supplementary Table S2.

Genes Detection
In this study, 37 significant SNPs were identified for BW trait, of which 33 were located within or near candidate genes through searching Ensembl, NCBI, and QTL databases. As shown in Table 1, significant SNPs on BAT5, BAT13, and BAT29 were within or close to genes ITGA5, ANGPTL4, and FGF4, respectively. Three significant SNPs were identified at the region of 26.11-33.58 Mb on BAT12, of which BvineHD1200027798 and BovineHD1200026844 were close to genes SHISA2 and CHMP3, respectively. No candidate genes were detected in four significant SNPs including BovineHD0200018897, BovineHD0700012290, BovineHD1100030552, and BovineHD2900008350, located on BAT2, BAT7, BAT11, BAT12, and BAT29, respectively. The detailed information on these genes was listed in Supplementary Table S2.

Functional Annotation of Candidate Genes
To further understand the function of candidate genes, GO and KEGG enrichment was performed using the R package "clusterProfiler" [29]. As shown in Figure 4, seven pathways were significantly enriched, of which five pathways were implicated in signal transduction, including PI3K-Akt signaling pathway (bta04151), Ras signaling pathway (bta04014), MAPK signaling pathway (bta04010), Rap1 signaling pathway (bta04015), and Calcium signaling pathway (bta04020); one pathway, regulation of actin cytoskeleton (bta04810), was involved in cell motility; one pathway was associated with lipid metabolism, namely alpha-Linolenic acid metabolism (bta00592). The information about these significant pathways is shown in Supplementary Table S3. No GO term was significantly enriched. Notably, FGF4, ANGPT4, ITGA5, and PLA2G4A involved in no less than two KEGG pathways deserved further attention and discussion. Additionally, based on animal QTLdb, the previous reports and gene function analysis, several candidate genes including TBC1D5, SHISA2, PDE1C, and COCX6 were also identified to affect the BW trait. The information of these candidate genes was listed in Table 2. Numerous candidate genes known to be related to BW in domestic animals were presented in Supplementary Table S4. PPI network shown in Figure 5 visualized the interaction between node proteins encoded by corresponding candidate genes. based on animal QTLdb, the previous reports and gene function analysis, several candidate genes including TBC1D5, SHISA2, PDE1C, and COCX6 were also identified to affect the BW trait. The information of these candidate genes was listed in Table 2. Numerous candidate genes known to be related to BW in domestic animals were presented in Supplementary Table S4. PPI network shown in Figure 5 visualized the interaction between node proteins encoded by corresponding candidate genes.   based on animal QTLdb, the previous reports and gene function analysis, several candidate genes including TBC1D5, SHISA2, PDE1C, and COCX6 were also identified to affect the BW trait. The information of these candidate genes was listed in Table 2. Numerous candidate genes known to be related to BW in domestic animals were presented in Supplementary Table S4. PPI network shown in Figure 5 visualized the interaction between node proteins encoded by corresponding candidate genes.  A represents the network diagram of the candidate genes affecting BW and their belonged pathways; B represents the PPI network of the candidate genes affecting BW. The genes marked in red represent potential candidate genes influencing BW in this study. Figure 5. A represents the network diagram of the candidate genes affecting BW and their belonged pathways; B represents the PPI network of the candidate genes affecting BW. The genes marked in red represent potential candidate genes influencing BW in this study.

Discussion
Longitudinal traits can better reflect the growth and development patterns of livestock and poultry. Therefore, mapping and analyzing the significant SNPs and functional genes affecting longitudinal traits have important economic value for beef cattle breeding. It is worth noting that GWAS-RRM is the main approach to analysis longitudinal traits, which could better control FPR and improve the accuracy of estimated breeding values (EBVs) [17,18], thus improving the efficiency of GWAS analysis [30]. Body weight is an important longitudinal trait that greatly reflects bovine growth performance. Previous studies have shown skeletal muscle is involved in the structure and metabolic regulation of the body and its mass accounts for 40% of total body weight in animals [31], indicating the growth and development of animals are inseparable from muscle development. The actin cytoskeleton is an important muscle structure that regulates cell adhesion, cell proliferation, cell motility, and muscle contraction via the signals transduction from the extracellular matrix to the nucleus [32]. Ito et al. elucidated that the MAPK signaling pathway could stimulate the growth of skeletal muscle and cell proliferation [33]. Therefore, it could be concluded that the pathway of regulation of actin cytoskeleton and the MAPK signaling pathway might be potential candidate pathways affecting the BW trait via affecting the development of skeletal muscle. In the present study, several candidate genes (FGF4, ANGPTL4, PLA2G4A, and ITGA5) were significantly enriched in the above pathways, which implied these genes might play special roles in regulating the BW trait.
Fibroblast growth factors (FGFs) are important growth factors that participate in many developmental and physiological processes [34]. FGF signaling could be disrupted due to the mutations of FGFs, thus causing developmental disorders, i.e., skeletal diseases, infertility, and cancer [35]. Relevant studies have reported some members of the FGF family, FGF4, and FGF14, regulate fibroblasts formation and the development of growth traits such as BW trait [36,37]. The complete nucleotide sequence of the bovine FGF4 was identified in three cattle breeds (panese Black, Japanese Shorthorn, and Holstein cattle) [38], and its coding exons encode 206 amino acid residues, perhaps including a signal peptide at the amino terminus [39]. Sato et al. reported FGF4 was related to the regulation of bovine embryo development [38]. Consistent with these findings, Feldman et al. also found FGF4 was associated with trophoblast proliferation in mice, and FGF4 null mice showed a peri-implantation lethal phenotype [40]. This evidence supported FGF4 was an important indicator in the growth and development of animals.
Angiopoietin-like 4 (ANGPTL4), known as a novel peroxisome proliferator-activated receptor target gene, is a key regulator of triglyceride, non-esterified fatty acid (NEFA) concentrations, and plasma cholesterol [41]. Just as in the mouse study, ANGPTL4 is also widely expressed in many bovine tissues, such as liver, subcutaneous adipose tissue, rumen, omasum, abomasum, etc., among which the first two are important tissues for ANGPTL4 synthesis [42]. Under fasting conditions, its expression level is strongly up-regulated in the liver and adipose tissue [43], which plays an important role in lipid metabolism via inhibition of the lipoprotein lipase (LPL) and stimulated lipolysis [44]. ANGPTL4 has been recognized as an adipokine in bovine adipose tissue and its expression could affect bovine body fat [45]. Notably, intramuscular fat becomes an important component of maturing muscles, and the mass of skeletal muscle accounts for 40% of total body weight in animals [32], thus it could be speculated that ANGPTL4 might regulate BW trait by influencing the formation of adipose tissue. Additionally, previous studies have reported ANGPTL2, a homologous family gene of ANGPTL4, influences the development of the bovine BW trait [24]. Taken together, ANGPTL4 could be recognized as the candidate gene regulating the BW trait for further research.
Phospholipase A2 (PLA2) is classified into three groups according to their chemical properties and molecular structure, namely cytosolic (cPLA2), secretory (sPLA2), and Ca 2+ -independent PLA2s [46]. Previous researchers found that activated cPLA2 could stimulate the release of arachidonic acid [46], which directly suppressed the growth and development of tumor cells [47]. Phosphorylation of cPLA2 induced by Temozolomide could also cause the suppression of cell growth [48]. More importantly, cPLA2 is an important regulator in various muscle development. Gluck et al. proposed cPLA2 activation was essential for the proliferation of bovine aortic smooth muscle cells [49]. In the work by Hirabayashi et al., cPLA2 alpha was identified to be responsible for striated muscle growth and fertility in mice [50]. In the present study, cPLA2, known as PLA2G4A, was significantly enriched in two signal transduction pathways, including the MAPK signaling pathway (bta04010) and Ras signaling pathway (bta04014). As mentioned above, the MAPK signaling pathway could stimulate the growth of skeletal muscle [33]. Hence, cPLA2 was forecasted to be the promising gene affecting the bovine BW trait via the involvement in muscle development.
Integrins are a family of heterodimeric cell-surface adhesion receptors that affect cell-matrix interaction. Some integrins encoding genes, including integrin alpha-2 (ITGA2) and integrin alpha-11 (ITGA11), have been proven to regulate the BW trait of swine and sheep, respectively [51,52]. In this study, integrin alpha-5 (ITGA5) was speculated to be associated with bovine BW trait. ITGA5 participated in various cellular processes, such as cell adhesion, survival, proliferation, differentiation, and migration of myoblasts [53], adipocytes [54], and cardiac neural crest [55]. Its differential expression was correlated with the organ specificity of tumor metastasis [56], thus ITGA5 was recognized as a potential biomarker for cancer treatment. Previous studies have shown ITGA5 could promote the proliferation, migration, and invasion of oral squamous cell carcinoma through activating the PI3K/AKT signaling pathway [57]. Chen et al. demonstrated that ITGA5 was a mediator for the proliferation and migration of retinal pigment epithelial cells [58]. ITGA5 knockdown or overexpression could inhibit or accelerate cell growth, respectively. Fang et al. reported ITGA5 participated in integrin β1 overexpression, which caused growth arrest of breast cancer cell [59]. Larzabal et al. revealed that suppressed ITGA5 could reduce adherence capacity to fibronectin and inhibit tumor growth in lung cancer cells [60]. In addition to cell growth, ITGA5 has been proposed to regulate porcine drip loss by mediating cell adhesion and extracellular matrix [61]. However, at present, there is no supporting evidence for ITGA5 on the weight of beef cattle, thus ITGA5 as a molecular marker for bovine growth traits needs further investigation.
Combined with previous studies and gene function analysis, except for candidate genes listed above, several known or potential candidate genes were also identified to affect bovine BW trait in this study. TBC1 domain family member 5 (TBC1D5), encoding GTPaseactivating protein (GAP) for Rab7, is a high-affinity ligand of the retromer cargo selective complex VPS26/VPS29/VPS35. Previous studies showed TBC1D5 was an important novel regulator that rerouted ATG9-containing vesicular carriers toward sites of autophagosome formation [62]. Bärlocher et al. illustrated TBC1D5 could promote the intracellular growth of L. pneumophila [63]. Notably, Zhuang et al. reported TBC1D5 might play a special role in BW at 18 months of age in Chinese Simmental beef cattle [11]. Consistent with these findings, TBC1D5 was also identified as the functional gene for the bovine BW trait in the present study. However, the molecular mechanism by which TBC1D5 influences bovine BW remains to be elucidated.
Previous studies have demonstrated that SHISA9 affected the growth and development traits such as pre-weaning gain in sheep and BW in beef cattle [24,64]. Protein shisa-2 homolog 2 (SHISA2) identified in this study belongs to the same family as SHISA9, which encodes an endoplasmic reticulum protein against both Wnt and FGF signaling to affect the development of Chicken and Xenopus embryos [65,66]. Liu et al. reported SHISA2 not only regulated F-actin distribution but also directly mediated the maturation of membrane protein for myoblast fusion. Its overexpression could inhibit myoblasts' proliferation but promote premature fusion [67]. Human SHISA2 overexpression led to increased cell growth and invasion [68]. In the work by Hu et al., SHISA2 was identified to be involved in growth and development in duck skeletal muscle [69]. As mentioned previously, skeletal muscle mass accounts for 40% of total body weight in animals [31], thus implying SHISA2 might regulate BW trait via affecting the development of bovine skeletal muscle.
Phosphodiesterase (PDE1C) regulates the stability of growth factor receptors such as PDGFRβ [70]. which is highly expressed in the human heart, cardiac myocytes, and mouse heart, but rarely expressed in mouse cardiac myocytes [71,72]. Cai et al. demonstrated that PDE1C positively regulated smooth muscle cells (SMCs) growth, proliferation, migration, and neointimal hyperplasia [73]. In agreement, the high expression of PDE1C was screened in proliferating human arterial SMCs in primary culture, but not in the quiescent SMCs [74], which indicated this gene was an indicator of cell proliferation. Notably, the previous study by Duan et al. found PDE1C might be the candidate gene affecting bovine BW trait through signal-trait GWAS [24]. However, many studies of PDE1C mainly rely on human and mouse SMCs, with less research on beef cattle, thus the functions of PDE1C in beef cattle should be further investigated.
Cytochrome c oxidase subunit 6C (COX6C) is eventually transported to mitochondria to form the cytochromec oxidase (COX) complex [75]. Duggan et al. suggested that COX participated in the remodeling of skeletal muscle [76]. The expression levels of COX subunits are different in vertebrate muscle [77]. COX6C overexpression could induce cell growth retardation [78]. Therefore, it was speculated that COX6C might be a valuable gene for bovine growth and development.

Conclusions
In conclusion, GWAS-RRM has been recognized as the main analysis model for longitudinal traits as it could decrease FPR and increase statistical powers. Based on this method, the present study mainly revealed four most promising candidate genes (FGF4, ANGPTL4, PLA2G4A, and ITGA5) and two significantly enriched pathways regulated bovine BW trait by affecting the growth and development of skeletal muscle and adipose tissue. The function of these candidate genes and pathways have been analyzed and discussed in detail. Moreover, further studies will be necessary to clarify their molecular mechanisms and physiological implications in regulating BW trait. This study not only offers molecular information for genomic selection of bovine growth and development traits but also provides the reference for the large-scale application of GWAS-RRM analysis of longitudinal traits in other livestock and poultry.

Informed Consent Statement: Not applicable.
Data Availability Statement: The raw data underlying our findings has been submitted to Dryad Digital Repository, which is available from the following link: https://doi.org/10.5061/dryad.4qc06, accessed on 1 July 2021.