Gene co-expression network and differential expression analyses of subcutaneous white adipose tissue reveal novel insights into the pathological mechanisms underlying ketosis in dairy cows

Ketosis is a common nutritional metabolic disease during the perinatal period in dairy cows. Although various risk factors have been identified, the molecular mechanism underlying ketosis remains elusive. In this study, subcutaneous white adipose tissue (sWAT) was biopsied for transcriptome sequencing on 10 Holstein cows with type II ketosis [blood β-hydroxybutyric acid (BHB) >1.4 mmol/L; Ket group] and another 10 cows without type II ketosis (BHB ≤1.4 mmol/L; Nket group) at d 10 after calving. Serum concentrations of nonesterified fatty acids (NEFA) and BHB, as indicators of excessive fat mobilization and circulating ketone bodies, respectively, were significantly higher in the Ket group than in the Nket group. Aspartate transaminase (AST) and total bilirubin (TBIL), as indicators of liver damage, were higher in the Ket group than in the Nket group. Weighted gene co-expression network analysis (WGCNA) of the sWAT transcriptome revealed modules significantly correlated with serum BHB, NEFA, AST, TBIL, and total cholesterol. The genes in these modules were enriched in the regulation of the lipid bio-synthesis process. Neurotrophic tyrosine kinase receptor type 2 ( NTRK2 ) was identified as the key hub gene by intramodular connectivity, gene significance, and module membership. Quantitative reverse transcription PCR analyses for these samples, as well as a set of independent samples, validated the downregulation of NTRK2 expression in the sWAT of dairy cows with type II ketosis. NTRK2 encodes tyrosine protein kinase receptor B (TrkB), which is a high-affinity receptor for brain-derived neurotrophic factor, suggesting that abnormal lipid mobilization in cows with type II ketosis might be associated with impaired central nervous sys-tem regulation of adipose tissue metabolism, providing a novel insight into the pathogenesis underlying type II ketosis in dairy cows.


INTRODUCTION
Ketosis is a common metabolic disease occurring in the early lactation of dairy cows. Cows with ketosis have decreased milk production and poor reproductive performance, leading to significant economic losses. Ketosis is characterized by a high concentration of serum BHB (Duffield, 2000;Suthar et al., 2013). The highest incidence of ketosis occurs at d 5 after parturition (McArt et al., 2013). Compared with the relatively lower incidence (2%-15%) of clinical ketosis (BHB ≥3.0 mmol/L), subclinical ketosis (1.4 mmol/L ≤ BHB <3.0 mmol/L) has a significantly higher incidence (8.9%-34%) in the first 2 mo of lactation (Duffield, 2000).
Ketosis is usually associated with negative energy balance (NEB) in early lactation, when energy intake is insufficient to meet energy requirements for milk synthesis. Type II ketosis, often characterized by fatty liver, is commonly found in high-producing dairy cows and occurs frequently in the first 1 to 2 wk of lactation, due to the failure of adaptation to NEB (Herdt, 2000). The pathology of ketosis is related to intense lipid mobilization from adipose tissue, leading to increased plasma concentrations of nonesterified fatty acids (NEFA). The accumulation of NEFA impairs stimulation of gluconeogenesis, thus increasing the risk of excessive production of BHB in liver (Grummer, 1993;Steeneveld et al., 2020). Therefore, as an energy-storing reservoir, adipose tissue plays a critical role in the process of the metabolic adaptation to type II ketosis (Tamura and Shimomura, 2005;Martin and Sauvant, 2007).
Elevated NEFA concentration in serum is typically associated with subclinical ketosis, which may be due to an excessive mobilization of NEFA from adipose tissue under severe NEB. The release of NEFA from adipose tissue is controlled by the balance of lipogenesis (fatty acid synthesis) and lipolysis (hydrolysis of triglycerides; Proença et al., 2014;Dragos et al., 2017). Adipose deposition is determined by the balance of lipolysis and lipogenesis, which is regulated by both endocrine and neural mechanisms. Previous gene expression studies have shown that increased lipolysis might not be a major factor for excessive lipid mobilization, because the expression and phosphorylation of hormone-sensitive lipase (HSL), a key enzyme for lipolysis in adipose, did not change in dairy cows with high serum NEFA (Mc-Namara, 1989;van Dorland et al., 2012;Mann et al., 2016). In contrast, Alharthi et al. (2018) proposed that high levels of circulating NEFA might be attributed to reduced lipogenesis. However, another study showed that no significant change in lipogenesis was observed in dairy cows with severe lipid accumulation in the liver during the postpartum period (Ji et al., 2012). In recent years, RNA sequencing (RNA-seq) has been used to comprehensively understand mechanisms of gene expression underlying distinct biological processes. Recent transcriptomic studies revealed that the network of lipid metabolism was significantly affected by NEB in dairy cows (Mellouk et al., 2019), but the key genes associated with this influence remained unidentified.
Weighted gene co-expression network analysis (WGCNA) of the transcriptome is a powerful tool for the establishment of gene regulatory networks and the identification of key gene modules related to diseases (Langfelder and Horvath, 2008). In this study, we aimed to perform WGCNA coupled with differential expression analysis of the subcutaneous white adipose tissue (sWAT) transcriptome to explore the key genes and pathways associated with type II ketosis in dairy cattle.

Ethics Statement
All animals were treated in accordance with the Guidelines for Laboratory Animal Use and Care from the Chinese Center for Disease Control and Prevention and the Rules for Medical Laboratory Animals (1998) from the Chinese Ministry of Health, under protocol (CAU-AEC-2012-351) approved by the Animal Ethics Committee of China Agricultural University.

Animals
This longitudinal cohort study was conducted in a large-scale intensive farm (Anhui Province, China) where the incidence of subclinical ketosis was 30% to 40% for postparturient Holstein cows. The experiment was designed to collect 10 cows each for the type II ketotic group (Ket) and the nonketotic group (Nket) to meet the sample size requirement of WGCNA analysis. Cows in the dry period were monitored by assessing BCS and measuring blood BHB at 21 d before the expected calving date. The cohort was formed by choosing a total of 30 cows with 3.5 ≤ BCS < 4.0 and BHB <1.4 mmol/L. Blood BHB was measured on these cows at d 5 and 10 postpartum. The cows were assigned into the Ket group (BHB >1.4 mmol/L at either 5 d or 10 d after calving) or the Nket group (BHB ≤1.4 mmol/L at both 5 d and 10 d after calving). Cows with any postpartum clinical diseases other than ketosis were excluded from the study. Once the number of cows reached 10 in either group, the rest of cows in that group were not included. Eventually, 10 cows were retained for each group subject to adipose tissue biopsy (Supplemental Table S1, https: / / doi .org/ 10 .17632/ yxg6zvgjnm .1; Zhao, 2023).
The cows were fed a TMR 3 times a day and had free access to water throughout the day. Blood samples were collected from the coccygeal vessels using BD Vacutainer Rapid Serum Tube (Becton, Dickinson and Co.). The concentration of BHB in the blood samples was determined using the FreeStyle Optium Neo H Ketone Meter (Abbott Diabetes Care Ltd.). This device is reported to have a sensitivity of 98% and a specificity of 95% (Macmillan et al., 2017).

Adipose Tissue Biopsy
The biopsy was done within 10 min to collect adipose tissue from the tail-head region of cows at 10 d after parturition. Briefly, anesthetic management was achieved using fan-shaped anesthesia application with procaine hydrochloride (3%, 25 mL). Then, an incision (1.3-2.5 cm) parallel to the spine was made to collect the sWAT (1-2 g) using a scalpel. The accidental cut-ting of blood vessels was avoided as much as possible. After washing the sWAT using normal saline stored at 4°C, the adipose tissue was immediately cut into small pieces (about 0.7 cm in diameter) using surgical scissors and tweezers. The samples were stabilized in RNAlater (Invitrogen, Thermo Fisher Scientific) and stored in a freezer at −20°C.

RNA-Seq
Trizol reagent (Invitrogen) was used to isolate the total RNA from the sWAT samples using the standard protocol with modifications. Briefly, a piece of adipose tissue of approximately 100 to 150 mg was cleaned using sterilized distilled deionized H 2 O (ddH 2 O), placed in a 1.8-mL centrifuge tube, and mixed with 1 mL of Trizol reagent and four 5-mm stainless steel beads. The sample was homogenized at 70 Hz for 5 min and the homogenate was set for 15 min, vortexed vigorously for 30 s, and centrifuged for 10 min at 12,000 × g and 25°C.
The lipid droplets were removed thoroughly, the supernatant was transferred to a new tube, and 0.2 mL of chloroform was added. The solution was then vortexed for 30 s, incubated for 15 min at 25°C, and centrifuged at 12,000 × g for 15 min at 4°C. The upper aqueous phase was carefully transferred to a new tube, followed by mixing with 0.4 mL of isopropyl alcohol. After incubation for 10 min and centrifugation at 12,000 × g for 10 min at 4°C, the RNA precipitate formed a pellet at the bottom of the tube. The resulting pellet was washed in 1 mL of 75% ethanol and centrifuged at 7,500 × g for 5 min. The supernatant was discarded, and the RNA pellet was air-dried and dissolved in 15 μL of sterilized ddH 2 O and finally stored at −80°C for use. A NanoDrop 1000 spectrophotometer (Thermo Scientific) was used to determine the purity and concentration of extracted RNA. An Agilent Bioanalyzer 2100 system (Agilent Technologies) was used to measure the integrity of the RNA samples. RNA samples with an optical density ratio at 260 and 280 nm (OD 260/280 ) between 1.9 and 2.1 and integrity number ≥8 were selected for RNA-seq.
RNA sequencing was performed on the Illumina NovaSeq 6000 platform (read length of 150 bp). Raw reads of RNA-seq in this study were deposited in the Sequence Read Archive database of the National Center for Biotechnology Information under accession number PRJNA835668. Then, FastQC v. 0.11 (https: / / www .bioinformatics .babraham .ac .uk/ projects/ fastqc/ ; accessed October 22, 2022) and STAR v2.7.5b (Dobin et al., 2013) were applied to check the quality of the raw data and align the clean data to the ARS-UCD1.2 bosTau9 reference genome. The average alignment rate was 95.57%. The read count matrix and the transcripts per million (TPM) matrix were constructed using RSEM V1.3.1 (Wang et al., 2022).

Construction of Co-Expression Gene Network and Identification of Modules Correlated with Type II Ketosis Using WGCNA
The TPM matrix was imported into R software v. 4.1.1 (R Core Team, 2021) and the first 7,000 genes based on the median absolute deviation were selected. A topological overlap matrix (TOM) was constructed using the WGCNA R package v. 1.70-3 (Langfelder and Horvath, 2008) to describe the weighted correlation between nodes (soft threshold β = 7, scale-free R 2 = 0.8010, mean connectivity = 100.0). Then, the TOM was converted to a dissimilarity measure (dis-sTOM, equal to 1 − TOM). Hierarchical clustering with dynamic tree cutting was performed to group the genes into modules based on their expression patterns. Modules with similar co-expression were merged based on the values of height (h = 0.14). Finally, the Pearson correlation coefficient between the module eigengene, which was defined as the first principal component of each module, and the serum biochemical indicators were calculated to identify the modules significantly associated with these phenotypes (P < 0.05; Langfelder and Horvath, 2008).

Identification of the Hub Gene in Modules Correlated with Phenotypes.
To identify the hub genes in the modules, 2 complementary methods were used. First, the CytoHubba v. 0.1 plug-in of Cytoscape (v. 3.7.1; https: / / cytoscape .org) was applied to search for candidate hub genes (top 10% of genes) according to the intramodular connectivity (kIN). Second, genes with high gene significance (GS) (|GS| > 0.60) and module membership (MM) (|MM| > 0.80) were selected as candidate hub genes of the module. Gene Ontology (GO) enrichment analysis was performed using the David online database (https: / / david .ncifcrf .gov/ ; accessed October 22, 2022).

Identification of Key Genes for Type II Ketosis in Postpartum Dairy
Cows. Differential expression analysis was conducted in DESeq2 R package v. 1.32.0 (https: / / www .bioconductor .org) using the read count matrix. Differentially expressed genes (DEG) were defined as | log 2 fold change) | >1, Benjamini-Hochberg adjusted P-value <0.05, and base-mean >100 (Love et al., 2014). The GO enrichment analysis was performed for DEG as described above. Overlapping genes between DEG and candidate hub genes from WGCNA analyses were considered key candidate genes for type II ketosis.

Validation of Key Candidate Genes
Reverse transcription-quantitative PCT (RT-qP-CR) was performed to validate the key candidate genes that were identified by transcriptomic analyses. First, 3 cows were randomly selected from each group (Ket and Nket) to validate the RNA-seq results. Second, a repeat experiment was done to collect adipose samples from 3 postparturient cows with ketosis and 3 without ketosis in an independent dairy herd.

Statistical Analyses
Student's t-test was performed in SPSS 26.0 (IBM Corp.) software to statistically compare the 2 groups in terms of serum biochemical indicators and gene expression levels. GraphPad Prism v. 9.3.0 (Graph Pad Software Inc.) was applied to conduct the receiver operating characteristic (ROC) curve analysis of NTRK2. To show expression levels of NTRK2 from both RNA-seq and RT-qPCR in figures, the results were standardized using the scale function in R software v. 4.1.1 (R Core Team, 2021). GraphPad Prism was used to visualize the standardized TPM and RT-qPCR results.

Serum Biochemical Indicators
Serum of cows in the Ket group had significantly higher BHB and NEFA concentrations compared with serum of Nket cows (P < 0.01, Table 1). Serum ALB, AST, TBIL, and TCH are commonly used as indicators of liver function. No significant difference was observed in serum ALB between the Ket and Nket groups (P > 0.05). However, compared with the Nket group, the Ket group had significantly higher serum AST and TBIL (P < 0.01) but lower TCH (P < 0.01; Table 1).

Establishment of a Gene Co-Expression Network and Identification of Modules Correlated with Type II Ketosis Using WGCNA
A total of 35 gene co-expression modules were constructed from the TOM ( Figure 1A). The number of genes in each module ranged from 31 to 922 ( Figure  1B). The correlations of module eigengene and serum biochemical indicators showed that the darkolivegreen module was significantly associated with mutiple biochemical indicators including BHB (r = −0.66, P = 0.002), NEFA (r = −0.62, P = 0.003), AST (r = −0.7, P = 6E-04), TBIL (r = −0.6, P = 0.005), and TCH (r = 0.59, P = 0.007; Figure 1B). Therefore this module was selected for subsequent analysis as the key module correlated with type II ketosis. There are 89 genes in the darkolivegreen module (Supplemental Table S2, https: / / doi .org/ 10 .17632/ yxg6zvgjnm .1; Zhao, 2023). Functional enrichment analysis showed that these genes were significantly enriched in the regulation of the lipid biosynthetic process and positive regulation of synaptic development (P < 0.05; Figure 1C; Supplemental Type II ketosis was defined as a concentration of BHB >1.4 mmol/L at either 5 or 10 d after calving (Duffield, 2000).
3 P-value of the Student's t-test.

Analysis of Candidate Hub Genes of the Module Correlated with Type II Ketosis
Based on intramodular gene connectivity (kIN), 9 genes were selected as hub genes (Figure 2A; Supplemental Table S4, https: / / doi .org/ 10 .17632/ yxg6zvgjnm .1; Zhao, 2023). These genes were significantly involved in the positive regulation of endothelial cell migration (BSG, EMC10), the glutathione metabolic process (GSTM2), and the regulation of GTPase activity (NTRK2, PLXND1; P < 0.05). Furthermore, accord-ing to GS and MM measurements, the hub genes associated with BHB, NEFA, and AST were determined. The 10 BHB-associated hub genes participate in the positive regulation of synaptic development (CAMK1, CNTN1, NTRK2) and cell adhesion (BSG, LAMA4, CNTN1, AOC3; P < 0.05) ( Figure 2B). The 10 hub genes correlated with NEFA participate in the positive regulation of synaptic development and the positive regulation of peptidyl-serine phosphorylation (CAMK1 and NTRK2; P < 0.05; Figure 2C). The 25 AST-associated hub genes participate in the positive regulation of synaptic development (CAMK1, CNTN1, NTRK2) and cell adhesion (BSG, LAMA4, CNTN1, AOC3; P < 0.05; Figure 2D; Supplemental Table S4; Zhao, 2023). Remarkably, NTRK2 was shared among all 4 sets of candidate hub genes, which were associated with serum BHB, NEFA, and AST, and it was highly correlated with the remaining genes in the darkolivegreen module.

Identification and Validation of Candidate Key Genes
Differential expression analysis detected a total of 52 DEG, 7 of which were upregulated and 45 downregu-  Table 2).
Notably, NTRK2 was significantly downregulated in the Ket group compared with the Nket group. Combining the differential expression and WGCNA analysis revealed NTRK2 to be the most promising candidate gene for type II ketosis in dairy cows ( Figure 3B). Additionally, using ROC curve analysis, we observed that the area under the curve of NTRK2 was 0.990 (SE = 0.016, P = 0.0002), indicating the high predictive ability of the NTRK2 expression in sWAT for type II ketosis ( Figure 3C). Figure 3D, compared with the Nket group, samples in the Ket group showed significantly decreased expression of NTRK2 (P < 0.01), including both TPM from RNA-seq and relative expression from RT-qPCR ( Figure 3D). Additionally, the expression of NTRK2 was validated by RT-qPCR in an independent Ning et al.: PATHOLOGY OF KETOSIS Figure 3. Identification and validation of the candidate key gene. (A) Volcano map; red = downregulated differentially expressed genes (DEG); blue = upregulated DEG; black = not significant genes (|log 2 (fold change)| >1, P adj < 0.05, base-mean >100). (B) Venn diagrams of candidate hub genes and DEG; NTRK2 was the only overlapped gene. (C) The receiver operating characteristic (ROC) curve was used to analyze the predictive ability of NTRK2 expression in subcutaneous white adipose tissue for type II ketosis in dairy cows. (D) Quantitative reverse transcription PCR (RT-qPCR) was used to verify the RNA-sequencing results (n = 3 for the Ket group, n = 3 for the Nket group); then, the decreased expression level of NTRK2 in the Ket group was validated using RT-qPCR in another batch of animals (n = 3 for the Ket group, n = 3 for the Nket group). TPM = transcripts per million. To show NTRK2 expression levels from both RNA-seq and RT-qPCR in the figure, the results were standardized; ** and * indicate P < 0.01 and P < 0.05, respectively. Error bars indicate the standard error of the mean. sample set (n = 3 for the Ket group, n = 3 for the Nket group). Consistently, samples from the Ket group showed lower NTRK2 expression compared with that from the Nket group (P < 0.05; Figure 3D).

DISCUSSION
Adipose tissue plays an important role in a cow's adaptation to metabolic switching during the periparturient period. In the present study, comparative transcriptome analysis of adipose tissues revealed a gene module significantly associated with serum biomarkers (BHB, NEFA, AST, TBIL, and TCH), with gene enrichment in the regulation of the lipid biosynthetic process. The hub gene NTRK2 in the module was identified as a key gene for type II ketosis.
NTRK2 encodes the tyrosine-protein kinase receptor B (TrkB), a high-affinity receptor for brain-derived neurotrophic factor (BDNF; Marosi and Mattson, 2014), identified in the cytosol and along the plasma membrane of adipocytes (Colitti et al., 2015). The abundance of TrkB is associated with the lipid metabolism in adipose tissue (Xu et al., 2003). Therefore, TrkB is considered an important regulator for energy homeostasis and significantly participates in the intake and comsumption of energy in the central nervous system and peripheral organs (Podyma et al., 2021). Additionally, the anorexia and reduced lipid storage that were observed in the adipose tissue of NTRK2knockout female mice under a high-energy diet (Nakagomi et al., 2015) is similar to the clinical symptoms of type II ketosis in dairy cows, such as reduced feed intake (Berge and Vertenten, 2014) and elevated blood NEFA. Based on this evidence, type II ketosis might be associated with central nervous system regulation of the peripheral lipid metabolism. The reduced expression of NTRK2 in white adipose tissue might impair the nervous regulation of lipid metabolism and disrupt the balance between lipogenesis and lipolysis in type II ketosis dairy cows.
Previous studies reported that the binding of BDNF to its receptor TrkB can activate its downstream pathways in cells, such as the mitogen-activated protein kinase (MAPK), phosphoinositide 3-kinase (PI3K), and the phospholipase Cγ (PLCγ) signaling pathways (Marosi and Mattson, 2014;Arosio et al., 2021;Podyma et al., 2021). The MAPK and PI3K signaling pathways play important roles in the metabolism of lipid in adipocytes (Prusty et al., 2002;Pirola et al., 2004;Liang et al., 2020). The PI3K signaling pathway plays a positive role in the prevention of type II ketosis by promoting the synthesis of lipids and secretion of adiponectin, and increasing the sensitivity to insulin (Pirola et al., 2004;del Rincon et al., 2007;McFadden, 2020). The MAPK signaling pathway regulates insulin resistance and is negatively associated with type II ketosis (Jager et al., 2007;Sumara et al., 2009;Ozaki et al., 2016). Therefore, reduced signaling of BDNF/ TrkB might impair its downstream pathways, leading to decreased lipid synthesis and elevated release of NEFA into the blood in dairy cows with type II ketosis.
Lipogenesis is largely regulated at the transcriptional level during early lactation (Bionaz and Loor, 2012;McNamara, 2012;Minuti et al., 2020). In this study, the downregulated genes in adipose tissue of the Ket group were most significantly enriched in the fatty acid biosynthetic process, including FASN, LPL, SCD, and ACACA. These genes play important roles in regulating fatty acid synthesis and are involved in the coordinated downregulation of lipogenesis during early lactation (Sumner-Thomson et al., 2011;Ji et al., 2012;Khan et al., 2013). In addition, THRSP and MID1IP1, 2 genes from the darkolivegreen coexpression module correlated with ketosis, were also enriched in the regulation of the lipid biosynthetic process. THRSP encodes thyroid hormone responsive spot 14 protein (Spot14), a nuclear protein that regulates fatty acid synthesis in the liver, adipose, and lactating mammary gland (Wellberg et al., 2014;Cui et al., 2015). MID1IP1 encodes MID1 interacting protein 1, which promotes the synthesis of fatty acids by enhancing the activity of acetyl-CoA carboxylase α (ACACA; Kim et al., 2010;Park et al., 2013). Based on our RNA-seq results, both THRSP (P < 0.1) and MID1IP1 (P < 0.05) showed a lower level of expression in the Ket group than in the Nket group (data not shown). Collectively, these findings revealed a more pronounced reduction of lipogenesis during early lactation in cows with type II ketosis compared with cows without type II ketosis.
In this study, we investigated the sWAT transcriptome to identify the key genes associated with subclinical ketosis. However, recent studies showed that both sWAT and abdominal adipose tissue (AAT) were mobilized to offset NEB during the perinatal period (Raschka et al., 2016). Cows stored about 2 to 3 times more fat in AAT than in sWAT depots during the dry period, and more adipose tissue mass was mobilized from AAT after parturition (Szura et al., 2020). Moreover, lipolytic activity is higher in AAT than in sWAT in the first weeks after calving (Saremi et al., 2014;Szura et al., 2020). These functional differences between sWAT and AAT may have a significant effect on the development of metabolic health disorders in dairy cows. Therefore, futher transcriptome studies on AAT would provide more insights into the molecular mechanism of ketosis in dairy cattle.

CONCLUSIONS
The gene co-expression network analysis and gene differential expression analysis of the sWAT transcriptome demonstrated that NTRK2 acts a potential key gene significantly affecting lipid metabolism in pluriparous Holstein cows. Cows with type II ketosis expressed significantly lower levels of NTRK2 in sWAT, supporting the hypothesis that dysfunctional lipid metabolism in cows with type II ketosis might be associated with central nervous system regulation of peripheral lipid metabolism. These results provide a new perspective on the pathological mechanisms underlying type II ketosis in dairy cows. However, further study on the regulatory role of NTRK2 at the cellular level and using different types of adipose tissue is warranted to clarify the pathogenesis of type II ketosis.

ACKNOWLEDGMENTS
This study was supported by the earmarked fund for CARS36 (Beijing, China), Ningxia Key Research and Development Program of China (Ningxia, China; 2021BEF02019; 2022BBF02017), and the National Key Research and Development Program of China (Beijing, China; 2021YFD1200904). We thank the veterinarians of Modern Farming (group) Co. Ltd. (Anhui, China) for their assistance collecting samples and J. S. F. Barker (University of New England, Armidale, Australia) for his assistance in editing the English in this article. We gratefully acknowledge the critical review of our manuscript by the anonymous reviewers. The authors have not stated any conflicts of interest.