Genome-wide identification and expression profiling analysis of Wnt family genes affecting adipocyte differentiation in cattle

The Wnt family features conserved glycoproteins that play roles in tissue regeneration, animal development and cell proliferation and differentiation. For its functional diversity and importance, this family has been studied in several species, but not in the Bovinae. Herein we identified 19 Wnt genes in cattle, and seven other species of Bovinae, and described their corresponding protein properties. Phylogenetic analysis clustered the 149 Wnt proteins in Bovinae, and 38 Wnt proteins from the human and mouse into 12 major clades. Wnt genes from the same subfamilies shared similar protein motif compositions and exon–intron patterns. Chromosomal distribution and collinearity analysis revealed that they were conservative in cattle and five species of Bovinae. RNA-seq data analysis indicated that Wnt genes exhibited tissue-specific expression in cattle. qPCR analysis revealed a unique expression pattern of each gene during bovine adipocytes differentiation. Finally, the comprehensive analysis indicated that Wnt2B may regulate adipose differentiation by activating FZD5, which is worthy of further study. Our study presents the first genome-wide study of the Wnt gene family in Bovinae, and lays the foundation for further functional characterization of this family in bovine adipocytes differentiation.


Introduction
Wnt proteins, the initiators of the Wnt signaling pathway, comprise a large family of secreted glycoproteins that rich in cysteine [41]. The transduction of Wnt signaling mainly includes three pathways: the canonical Wnt (Wnt/β-catenin) signaling pathway, the non-canonical Wnt/Ca2+ pathway, and the non-canonical planar cell polarity (PCP) pathway. The three patterns all require to bind to transmembrane receptor Frizzled to regulate intracellular responses. In the canonical Wnt pathway, the activated transmembrane receptor Frizzled (Fzd) lead to a stable accumulation of β-catenin in the cytoplasm and translocation to the karyon. Then β-catenin bind to the transcription factor T-cell factor/lymphoid enhancing factor (LEF1/TCF) family, and activate the transcription of target genes to regulate embryo development, tissue regeneration, and cell proliferation and differentiation [24,36,56,64]. In the non-canonical Wnt/Ca2+ pathway, the activated receptor Fzd (mainly Fzd2), results in intracellular Ca2+ release and irritates Ca2+-calmodulin-dependent protein kinase II (CamKII) calnexin (CaN), and protein kinase C (PKC), thus regulating cell adhesion and gene transcription [14]. In the non-canonical Wnt/PCP pathway, when the receptor Fzd is stimulated by Wnt proteins (e.g., Wnt5a and Wnt11), signals transmission are through Disheveled (Dvl) to trimeric G proteins. This is followed by activating downstream target genes Rho-associated kinase (Rock) and Jun N-terminal serine/threonine kinase (JNK), thereby regulating cytoskeletal actin and cell polarity [2,63].
Furthermore, studies have shown that Wnt signaling plays an important role in maintaining the undifferentiated state of precursor adipocytes and inhibiting adipogenesis. In 3T3-L1 preadipocytes, Wnt1 ectopic expression stabilizes β-catenin, thereby activating TCF-dependent gene transcription and blocking adipogenesis [7,52]. Glucagon Like Peptide-1 (GLP-1) promotes adipogenesis through the up-regulation of the Wnt4 and β-catenin in the canonical Wnt signaling pathway [38]. KDM5A interacts with C/EBPβ and cooperatively inhibits the transcription of Wnt6, leading to the inhibition of Wnt/β-catenin pathway and promotion of 3T3-L1 preadipocyte differentiation [22]. In human bone marrow stromal (mesenchymal) stem cell (hMSC), the treatment with Wnt3A activates Wnt/β-catenin signaling pathway, thus decreasing in adipogenesis and increasing in osteogenesis [50,51]. An activation of Wnt10B activates the Wnt signaling cascade and prevents the 3T3-L1 preadipocytes normal differentiation by inhibiting expression of C/EBPα and PPARγ [23,26,52]. In porcine adipose-derived mesenchymal stem cells (AMSCs), Wnt3A inhibits the potential of adipogenic differentiation and alters the cell fate from adipocytes to osteoblasts through Wnt/β-catenin signaling pathway [34]. In murine embryonic mesenchymal cell line C3H10T1/2, Wnt3-Fz1 chimera is an inhibitor of differentiation into the adipocyte lineage and a potent activator of differentiation into osteoblasts [8]. It has been reported Wnt gene family also functions via non-canonical Wnt signaling pathway in addition to its role as a canonical Wnt ligand. In hMSC, an inhibition of Wnt3A can suppress the non-canonical Wnt/JNK pathway and enhance adipocyte differentiation whereas its activation enhance osteoblast differentiation [51]. In 3T3-L1 preadipocytes, Wnt4 and Wnt5A positively regulate adipogenesis at the initial stage of the differentiation process by activation of PKC and calcium/calmodulindependent kinase II [44]. Wnt5B partially represses canonical Wnt/β-catenin signaling pathway and enhances adipogenesis [27]. These ndings stimulated our interest and guided us to explore the function of Wnt gene family on bovine adipocyte differentiation.
So far, the Wnt family has been extensively studied in some species (Drosophila melanogaster, Tribolium castaneum, Acyrthosiphon pisum, Anopheles gambiae, and Apis mellifera, etc. ) [5,11,15,39,40,43]. However, there was limited understanding on their expression patterns and regulatory mechanisms in bovine adipocytes differentiation. To investigate the evolution of Wnt family genes in Bovinae and elucidate the function in bovine adipocytes differentiation, we carried out a genome-wide identi cation and evolutionary analysis of Wnt gene family in eight species of Bovinae. And the expression pro les in different tissues and stages during adipocytes differentiation were also analyzed in cattle based on transcriptome data and qPCR, respectively. Our study provides a basis for understanding the distribution of Wnt genes in Bovinae and is bene cial to further study of potential function of Wnt genes in bovine adipocytes differentiation.
The amino acid lengths of the 19 bovine Wnt proteins ranged from 333 (Wnt8B) to 585 (Wnt4), whereas molecular weight (Mw) ranged from 36.61 to 63.17 kDa, consistent with protein length. Except Wnt3, the rest Wnt family proteins had the isoelectric points (pI) that higher than 8.0, as they contained more basic amino acids than acidic amino acids. Wnt3 was neutral, with a pI of 7.73. All 19 Wnt proteins contained the WNT conserved domain (Supplementary info File 4).

Structural features of bovine Wnt family members
In order to learn about the structural characteristics of bovine Wnt proteins and genes, we projected the conserved motifs and gene structures based on their phylogenetic rela tionships (Fig. 1). According to the evolutionary clades, the 19 bovine Wnt family menbers clustered into six main subfamilies ( -). All the Wnt family proteins shared six conserved domains termed motifs 1, 2, Introns, coding sequences (CDS) and untrans lated regions (UTR) were variable among the Wnt gene family. For instance, the length of Wnt genes ranged from 3,084 nt (Wnt1) to 64,231 nt (Wnt7A), mainly due to intron variation. The number of CDS varied from 3 to 6, and the length and layout of the non coding areas, 3'UTR and 5'UTR, were also variable. Despite this variability in CDS, introns and UTRs, the Wnt members in the same evolutionary subfamily tend to possess similar patterns in gene structures and conserved motifs.
Phylogenetic relationship of Wnt proteins in different organisms Phylogenetic analysis can provide a reference for understanding functional diversi cation of the Wnt family in Bovinae. So, phylo genetic analysis was conducted, which included eight Bovinae species (Bos taurus, Bos indicus, Hybrid-Bos taurus, Hybrid-Bos Indicus, Bos grunniens, Bubalus bubalis, Bos mutus, and Bison bison bison) and Homo sapiens and Mus musculus. Since Wnt proteins in model organisms (human and mouse) have been studied extensively, they were also included in this study. Finally, 186 amino acid sequences from these ten species were aligned to generate a nonrooted Neighbor-Joining (NJ) tree (Fig.  2). It revealed that Wnt family proteins could be subdivided into 12 proposed subfamilies, including Wnt1-11 and Wnt16 subfamily. There were seven subfamilies that contain two Wnt members, they are subfamily (Wnt7A and Wnt7B), (Wnt3 and Wnt3A), (Wnt2 and Wnt2B), (Wnt5A and Wnt5B), (Wnt10A and Wnt10B), (Wnt9A and Wnt9B), (Wnt8A and Wnt8B), respectively.
Chromosomal distribution and collinearity analysis of Wnt genes Wnt genes were mapped on nine chromo somes in organisms of Bovinae (Fig. 3). The bovine Wnt genes showed a similar distribution with the other ve species. However, the order of Wnt1 (30832513-30835596 bp) and Wnt10B (30841487-30847512 bp) in Chr 5 and Wnt3A (3035810-3087361 bp) and Wnt9A (3155441-3163239 Mb) in Chr 7 of Bos taurus was reversed from that in Bos grunniens. In addition, compared with Bos taurus, Wnt9B and Wnt16 were lacking in Bos Indicus.
Genome collinearity analysis revealed a satisfactory corresponding relationship between chromosomes of Bos taurus and Hybrid-Bos Indicus, Hybrid-Bos taurus, Bos indicus and Bos grunniens, respectively (Fig. 4a). Although the chromosome number is different between cattle (2 N = 60) and buffalo (2 N = 50), there is also large chromosome homology between these two species. Also, collinearity modules explain the difference in the position of Wnt gene family between cattle and the other ve species in Bovinae. For instance, the position variation of Wnt2B, Wnt11, Wnt1 and Wnt10B between Bos taurus and Bos grunniens may be caused by complex intra-chromosomal translocation events (Fig. 4b), whereas Wnt3 and Wnt9B are distributed on different chromosomes between cattle and buffalo (bovine Chr 19 and buffalo Chr 3), which may be caused by inter-chromosomal rupture or fusion during the evolutionary process (Supplementary info File 6).

Expression analysis of Wnt genes in different tissue
The functionally related genes tend to show a co-expression patterns and these genes often regulate biological processes collaboratively. To explore the expression patterns of the Wnt gene family during adipogenic differentiation, we investigated their expression levels in 163 samples of 60 tissue types. The Wnt genes along with other 13 closely related genes can be classi ed into four groups ( to ) (Fig. 5a) according to their differential expression patterns in diverse tissues. Accordingly, the 60 bovine tissue types also clustered into four main clades (a-d) based on the expression patterns of all the 31 genes including Wnt family. The members of Wnt family and its receptors, FZD gene family, displayed similar and overlapping expression patterns in the 60 tissues, suggesting their broad and coordinated regulatory role in life activities.
PPARγ, a marker gene for adipocyte differentiation, showed high expression in Group a which included omental, intramuscula, subcutaneous and mammary gland fats. And CTNNB1, FZD1, FZD5, FZD6 and Wnt2B clustered and showed the most similar expression pattern with PPARγ. Further analysis of the ve different fat tissues revealed that CTNNB1, a core gene of Wnt signaling pathway, showed high expression and similar expression patterns with PPARγ (Fig. 5b).

Isolation and induced differentiation of bovine primary adipocytes
To explore the expression patterns of the Wnt gene family during adipocyte differentiation, primary adipocytes collected from subcutaneous adipose tissue of calves were induced. The results of oil red O staining showed that lipid droplet gathered together in adipocytes induced for 10 days compared to preadipocytes (Fig. 6a). Further analysis revealed that the absorbance of the differentiated adipocytes at 260nm was signi cantly higher than that of the preadipocytes (Fig. 6b). Also, The adipogenic marker genes (PPARγ, CEBPβ, FABP4 and LPL) were up-regulated (Fig. 6c). These results indicate that an induced differentiation system of primary adipocytes was established successfully, which can be used in the subsequent gene expression analysis.
Expression analysis of Wnt genes during adipocyte differentiation A qPCR analysis was conducted to detect the expression of Wnt genes and their receptors (FZD genes) at 0, 2, 6, and 10 days during adipocytes differentiation (Fig. 7). The members of Wnt8B, Wnt11, Wnt16 and their receptors Fzd1, Fzd2, Fzd3, Fzd4, Fzd6 showed a relatively high expression in preadipocytes and then reduced with the process of induced differentiation, suggesting that they may collectively play a role in keeping adipocytes undifferentiated. The members of Wnt2, Wnt6, Wnt9B, Wnt10A and their receptors Fzd9, Fzd10 showed an up-regulated tendency, indicating that they may play a regulatory role during adipocyte differentiation. Furthermore, Wnt2B, Wnt4, Wnt8A and Fzd5, Fzd8 reached the lowest expression on the second day and displayed a similar overall trend of expression.

Discussion
Structural features of bovine Wnt family proteins and genes The proteins' core motifs and domains determine their function and activity [57]. Gene families usually encode a group of proteins that possess similar active sites and play a regulatory role synergistically, and they tend to have similar conserved motifs [19]. There are six conserved amino acid sequences (Motifs 1, 2, 4, 5, 6 and 7) in all the 19 bovine Wnt family members, pointing to a common functional site. Four Wnt family members (Wnt2, Wnt3, Wnt5A, and Wnt5A) have all of the ten motifs, whereas the other 15 members lack of 1 to 4 motifs. It is likely that these four motifs (Motifs 3, 8, 9 or 10) are not located at the core of the Wnt proteins domain.
As the introns and UTRs vary in length and lay out, the CDSs distribution of Wnt genes is also variable. Whereas their amino acid sequences and motif patterns, especially in conserved motifs, were similar, and they all have the WNT conserved domain, which play essential roles in keeping their three-dimensional structure and binding function.

Phylogenetic relationship of Wnt family proteins
Phylogenetic analysis provides a credible way to explore the relationship between acid sequence similarity and function of proteins in the same family [54]. By focusing on the phylogenetic relationship of Wnt family proteins in previous studies of multicellular eukaryotes, it showed that they can divided into 13 subfamilies. For instance, a total of 11 (Zhikong scallop), 12 (Yesso scallop, Paci c oyster), 13 (Lingula anatine, Plathynereis dumerlii, Lottia gigantean, Crassostrea gigas, etc.) subfamilies were identi ed, respectively [39]. In Bovinae, Wnt family proteins were classi ed into 12 sunfamilies lacking of WntA. It is consistent with the result that vertebrates have reserved all subfamilies except WntA [21,49]. Although the Wnt family is highly conserved in molecular function, large-scale analysis showed that several Wnt members have been lost in many species after the complete set of Wnt genes emerged in cnidarians [31,33]. Among the eight species in Bovinae, Wnt7B is missing in Bison, and Wnt9B and Wnt16 is missing in Bos indicus. The Wnt9 subfamily speci c to Bilateria was also found missing in Chlamys farreri [39]. It remains unclear whether these genes were failed to identify due to genome assembly limitation or they have been lost during species evolution. In the phylogenetic relationship, two genes from distinct species located in the same clade were de ned as orthologs [29]. The orthologous gene pairs among cattle and the other ve bovine species were identi ed based on the homologous relationships. Results showed that the orthologous Wnt members rst clustered in one clade, indicating that they were conserved among different species.

Collinearity analysis of Wnts in Bovidae
The gene duplication events including tandem duplication and segmental duplication can cause gene family expansion in genome evolution [32]. Chromosome distribution and collinearity analysis indicated that both tandem duplication and segmental duplication were the reasons for the expansion of Wnt family in Bovinae (Fig. 3-4). The members of the Wnt gene family distributed in nine chromosomes in the six selected species. Due to the different starting point of chromosome annotation among species, the arrangement of genes may be totally reversed. For instance, the order of Wnt6, Wnt10A and Wnt4 in Chr2, Wnt2 and Wnt16 in Chr4, and Wnt3A, Wnt9A and Wnt8A in Chr 7 was totally opposite between Bos taurus and Hybrid-Bos taurus. Besides, intra-chromosomal translocation and rearrangement during species evolution also lead to the changes in gene arrangement [45,61]. For example, the position of Wnt2B between Bos taurus and Bos grunniens altered due to the inversion of large segments within the chromosome (Fig. 4B). It lead to the positions where Wnt genes located change from collinear (conserved in the same order) to syntenic (not neces sarily in the same order) between two species [37]. In addition, buffalo Chr 1-5 are collinear with bovine Chr 27 and Chr 1, Chr 23 and Chr 2, Chr 19 and Chr 8, Chr 5 and Chr 28, and Chr 16 and Chr 29, respectively (Supplementary info File 6). It is speculated to be caused by inter-chromosomal rupture or fusion during the evolutionary process. In a word, genome-wide collinearity analysis of Wnts provided rich perspectives for studying the function and evolution of genes in Bovidae.

Wnt genes affecting adipocyte differentiation
Wnt family proteins can activate their receptors, Fzd family proteins, and regulate adipocyte differentiation through the canonical and/or non-canonical Wnt signaling pathway. As it's known to all, functionally related genes usually exhibit similar expression patterns and gene expression clustering analysis could group together genes of similar function e ciently [17]. Previous studies have shown that Wnt is selective in recognizing its Fzd receptors [58]. For instance, Wnt3 formed a chimera with Fzd1 to regulate the canonical Wnt signaling [8]. To dissect the expression pattern of this family, the expression of 19 Wnt members and 10 Fzd members in 60 tissue types was analyzed. The overlapping in expression patterns may suggest their coordinated and selective regulatory role (Fig. 5). Group (PPARG, CTNNB1, FZD1, FZD5, FZD6 and WNT2B) were highly expressed in four fat tissues (omental, intramuscular, subcutaneous and mammary gland). Meanwhile, the studies in human showed that Wnt2B and FZD5 have physical interaction and co-expression relationship (Table 2). Furthermore, they displayed similar expression tendency during the differentiation of bovine adipocytes (Fig. 7). As PPARG is a marker gene for adipocyte differentiation and CTNNB1 is a core gene for canonical Wnt signaling pathway, it is hypothesized that WNT2B may bind to its receptor FZD5 and regulate adipogenic differentiation through canonical Wnt signaling pathway.
Adipogenic differentiation is a well-organized and complicated process regulated by various genes. Analysis of the interaction relationship of Wnt and Fzd is essential to explore their roles. An integrated network of the Wnt and Fzd gene family and their interacting genes was constructed by STRING (https://string-db.org/) and visualized by Cytoscape (Supplementary info File 7) [53]. In order to ensure the accuracy of the interaction network, only the evidences of literature mining and experimental veri cation are selected. Results showed there are extensive and complex direct or indirect relationships between Wnt and Fzd gene family. Since there are limited study on the interaction between Wnt gene and FZD gene to regulate Wnt signaling pathway in cattle, GeneMANIA (https://genemania.org/) was used to mine their relationship in human. Results displayed the bias of Wnt family members in recognizing their transmembrane receptors, Fzd family proteins (Table 2) [9,58], which provided reference on the study of cattle. Overall, these results revealed that the Wnt genes interact with Fzd genes, activating the canonical and/or non-canonical Wnt signaling pathway, and collectively regu lating adipocyte differentiation. It laid the foundation for further research on Wnt genes regu lating adipocyte differentiation in cattle.

Ethics statement
Animal experiments were conducted in accordance with the requirement of the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, 2004). It is authorized by the Animal Ethics Committee of Ningxia University (permit number NXUC20200521). The calf used in the experiment was electric shocked before released and the primary adipocytes was isolated immediately, making all efforts to minimize its suffering. It also conforms to the requirements of American Veterinary Medical Association (AVMA) Guidelines. The reporting in the study follows the recommendations in the ARRIVE guidelines [46].

Analysis of Wnt Gene Characteristics
Protein properties of identi ed Wnt genes, including protein length, molecular weights, and isoelectric points, were estimated using ExPASy (https://web.expasy.org/protparam/) [4]. The conserved motifs were detected in MEME 5.0 [6]. The minimum and maximum number of amino acids in each motif was 6 and 50, respectively. Exon-intron structures and motif patterns of Wnt family were displayed using TBtools software [12]. NCBI-CDD [42] was used to identify the conservative domain compositions of 19 Wnt proteins, and the results were visualized by the TBtools software [12]. Multiple sequence alignments of Wnts and a phylogenetic Neighbor-Joining tree was constructed using the MEGA 7.0 software [30].

Phylogenetic Tree Construction, Chromosomal distribution, and Collinearity analysis
To investigate the phylogenetic relationship of Wnt genes, the protein sequences were aligned using ClustalW. And a Neighbor-Joining tree was constructed in MEGA 7.0 [30] using the following parameters: bootstrap method (1000 replicates), Poisson model, and complete deletion. FigTree software (version 1.4.3) was used to adjust and beautify the evolutionary tree. Chromosomal locations of Wnt genes in cattle and other ve species in Bovinae were obtained from general feature format (GFF3) les. Gene Location Visualize from GFF [12]was used to map the distribution of Wnt genes. Collinearity analysis of ortholog genes between Bos taurus and the other ve bovine species was performed using the MCScanX toolkit [59]. Subsequently, genome collinearity results and orthologous Wnts were visualized by Dual Systeny Plot for MCscanX [12].
Expression Pro les of Wnt Genes in RNA-Seq RNA-Seq data of 163 tissue samples were collected from the Ruminant Genome Database (http://animal.nwsuaf.edu.cn/code/index.php/Ruminantia) [13] and all the raw data was deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (Supplementary info File 8). The sequencing quality was checked using FastQC [62]. Quality control of raw sequence data was performed using the trimmomatic-0.36 [10]. Clean reads were then mapped to the Bos taurus genome reference using STAR [16] and Hisat2 [28]. FPKMs (Fragments Per Kilobase per Million mapped reads) of the genes in each sample were computed by Ballgown (Version 2.2.0) [20,47]. The heat map was constructed using theTBtools software [12].
Isolation, culture and induction differentiation of bovine primary adipocytes Primary adipocytes were isolated and cultured from subcutaneous adipose tissue of the calf in the Zerui ecological breeding farm. The primary adipocytes were isolated and cultivated by Type collagenase digestion method. The induction of primary adipocytes differentiation [25] and oil red O staining [60] were performed as previous described. The substance extracted from adipocytes induced at day 0 day and 10 day were also measured at 510nm absorbance with isopropanol as control.

RNA extraction and quantitative RT-PCR (qPCR)
The primers of the Wnt genes were designed using Primer Premier 5.0 software (Supplementary info File 9). Total RNA was extracted by the phenol-chloroform method using TRIzol (9109, Takara) and samples with an OD260/OD280 absorbance ratio between 1.8 and 2.0 were used in the subsequent experiments. Then, a total of 1000 ng RNA was reverse transcribed using random primers with Moloney murine leukemia virus reverse transcriptase (Takara Bio, Kyoto, Japan). Realtime PCR was carried out in a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) with SYBR Green Master Mix (Takara Bio, Kyoto, Japan). The qPCR reaction procedure was 40 cycles of pre-denaturation for 3 min (95 ℃), denaturation for 10 s (95 ℃), annealing for 20 s (60 ℃), extension for 30 s (72 ℃). Relative expression was calculated using the 2−∆∆Ct method [1,3] with βactin as the reference. Three replicates were performed for each test. Statistical signi cance determined using Graphpad Prism 7.0 software.    Chromosomal distribution of Wnt genes. The black font on the left represents chromosome numbers and the red font on the right represents Wnt genes.

Figure 4
Collinearity analysis of Wnt genes between cattle and other organisms. Syntenic genes pairs are linked by grey lines whereas syntenic Wnt genes are shown as red lines.