The MyoD1 Promoted Muscle Differentiation and Generation by Activating CCND2 in Guanling Cattle

Simple Summary Muscle proliferation and differentiation is a complex process, which is mainly regulated by the myogenic regulatory factors (MRFs) gene family. As one of its family members, myogenic differentiation 1 (MyoD1) is an important regulator of myoblast differentiation and fibrogenesis, and its main function is to regulate the proliferation of myoblasts and the directed development of myogenic cell lines. In this study, CRISPR/CAS9 technology was used to construct the MyoD1 gene knockout MDBK cell lines, and then the transcriptome of MDBK cells was analyzed based on high-throughput sequencing technology. The differentially expressed genes in the transcriptome and the functional enrichment of differentially expressed genes were analyzed. The differentially expressed genes in the transcriptome and the functional enrichment were analyzed to provide experimental data for further studies on the mechanism of MyoD1 regulating CCND2 in myocyte differentiation. Abstract The purpose of this study was to analyze the transcriptome of MyoD1 gene knockout MDBK cells (bovine kidney cells) using high-throughput sequencing. For the first time, CRISPR/CAS9 technology was used to construct a MyoD1 knockout in MDBK cells and transcriptome sequence analysis was used to examine MyoD1-related target gene expression. Transcriptome sequencing indicated the presence of 723 differentially expressed genes (DEGs) by comparing wild type and MyoD1 knockout MDBK cells and included 178 upregulated and 72 downregulated genes. The DEGs are mainly enriched in Pl-3-kinase and AKT, p53 signaling pathways. Quantitative RT-PCR confirmed that PDE1B, ADAMTS1, DPT, and CCND2 were highly expressed in the leg muscle, longissimus dorsi, and shoulder of Guanling cattle, and CCND2 was inhibited after MyoD1 knockout, suggesting it may be a key downstream gene of MyoD1 and associated with muscle formation and differentiation in Guanling cattle. This provides experimental data for subsequent studies on the regulatory mechanisms of muscle differentiation in Guanling cattle.


Introduction
Muscle development is regulated by members of the myogenic regulatory factor (MRF) gene family that regulate the proliferation and differentiation of muscle cells and fibers during embryonic as well as postnatal development. Myogenic differentiation 1 (MyoD1), an MRF member, is a key regulator of myoblast differentiation and muscle fiber formation [1]. MyoD1 is a basic helix-loop-helix (bHLH) protein that acts in the nucleus and possesses an amino-terminal transcriptional activation domain and a carboxyl-terminal α-helical domain [2,3]. The protein is activated by Ser, Thr, and Tyr phosphorylation. The bHLH domain specifically binds and interacts with proteins possessing helix-loop-helix domains, and the adjacent basic region is necessary to bind to the promoters or enhancers of muscle-specific genes such as CK and myogenin. The presence of MyoD1 can promote the specific transcription of target genes by forming an active heterodimer with proteins of the E protein family containing the bHLH domain and then combining with the conserved sequence E-box of the promoter or enhancer region of the target gene [4,5].
MyoD1 functions as the primary switch for muscle-specific gene activation and transcription and can recruit the histone lysine methyl-and acetyltransferases Set7 and p300 that perform H3K4me1 and H3K27ac modifications necessary for stem cell differentiation into myoblasts. In addition, MyoD1 also regulates muscle development through cooperation with the transcription factors c-Jun, Jdp2, Meis, and Runx1 that then enhance muscle-related gene expression [6]. For example, when GFP was fused with the bFGF gene, it was found that GFP positive myoblasts existed in the recipient muscle, which promoted the ability of muscle fiber regeneration and significantly increased the expression of myofibril protein [7]. MyoD1 also induces expression of miR-494-3p that inhibits p300 gene expression and myosin (MYH2) expression in multifunctional stem cells during skeletal muscle development and inhibits skeletal muscle differentiation [8]. Currently, significant progress has been made in the study of the MyoD1 gene in pigs, mice, and chickens [9][10][11], but it is difficult to figure out how to promote bovine muscle cell production.
Guanling cattle are a local fine breed in Guizhou Province, one of the 78 local livestock breeds under national key protection, ranking first among the four major yellow cattle in Guizhou with a high reproduction rate, high slaughter rate, high meat yield, high amino acid content, high protein content, low fat content, and other characteristics, with high food value, economic value, and development value. Previous studies of Guanling cattle have utilized promoter and transcription factor analyses combined with yeast one-hybrid systems to demonstrate that Myf6 strongly activated the MyoD1 promoter [12]. However, the role of MyoD1 in muscle differentiation is unclear. In our pre-experiments, we found that knockdown of the MyoD1 gene inhibited the growth of muscle cells in Guanyang cattle, so we developed an experimental system employing bovine kidney cells (MDBK) to examine the effects of a MyoD1 knockout more closely. We employed the CRISPR/CAS9 technique and used transcriptome analysis of MyoD1 knockout MDBK cells to identify differentially expressed genes (DEG). Our research focus was to understand the role MyoD1 plays in muscle cell differentiation.

Cultivation of MDBK
MDBK was purchased from the Shanghai Institute of Cell Biology, Chinese Academy of Sciences, and was cultured in high-sugar DMEM medium.

sgRNA Construction
The bovine MyoD1 gene sequence (Acc. No. NM_001040478.2) was used to design CRISPR gene knockouts using the online resource http://www.e-crisp.org/ (accessed on 10 April 2021). The 4 sgRNAs possessing the highest specificity and annotation scores were synthesized by Jin Kairui Bioengineering (Wuhan, China) (see Table 1). Primers for MyoD1 sgRNA were combined with DNA oligonucleotide annealing buffer and the LentiCrispr v2 vector (YaJi Biological of Shanghai, China) at 16 • C for 12 h using directions provided by the kit manufacturer (Thermo Fisher Scientific, Waltham, MA, USA). The linked products were transformed into Escherichia coli DH5α competent cells (Sangon Biotech, Shanghai, China) and single colonies were picked from LB ampicillin plates. Recombinant plasmids were extracted using an endotoxin-free extraction kit (Omega Bio-Tek, Dallas, TX, USA) and commercially sequenced to confirm their identity (Wuhan Bio-Engineering).

MyoD1 Gene Knockouts and Validation
We utilized lentivirus packaging of the plasmid constructs since preliminary experiments with MDBK cells (Institute of Cell Biology, Chinese Academy of Sciences, Shanghai, China) [13] had indicated that Lipofectamine 2000 (Invitrogen) transfection efficiencies were <5%. We therefore co-transfected 293T cells (Shanghai Biological Technology Co., Ltd. Enzymere, Shanghai, China) with the lentivirus packaging plasmid psPAX2 and lentiviral vector pMD2.G at 2:1:1) that were then incubated for 5 h according to the manufacturer's specifications (Omega Bio-Tek, Dallas, Texas, USA). The medium was exchanged, and incubation was continued for 48 h. The lysed cells were filtered (0.45 µm) to obtain a 6 mL virus solution, and 3 mL was added to MDBK cells in a 6-well plate (4 sg types in total) and incubated for 8 h. The medium was changed, and the cells were incubated for an additional 24 h. The cells were then re-infected with the remaining 3 mL of virus solution and incubated as per above. The cells were added to fresh medium containing 2 µg/mL puromycin (ChemeGen, Los Angeles, CA, USA) and incubated for 60 h. Samples were taken for genomic DNA extraction and the remainder of the cells were expanded to 24-well plates and genomic DNA was again extracted to identify the knockout effect using PCR. PCR reactions utilized 2 × Taq PCR Master Mix (Sangon Biotech, Shanghai, China) using 10 pmol of primer and 0.5 µL of template DNA. The amplification conditions were as follows: 95 • C for 5 min, followed by 30 cycles of 95 • C for 30 s, 58 • C for 30 s, and 72 • C for 20 s, with a final incubation at 72 • C for 5 min. The amplicons were visualized by electrophoresis in 2% agarose gels and stained with SYBR green. Further verification using qRT-PCR utilized a commercial SYBR green amplification kit (Bio-Rad Laboratories, Hercules, CA, USA) containing 5 pmol of each primer ( Table 2) and 1 µL template DNA. The reactions were carried out at 95 • C for 1 min, then 40 cycles at 95 • C for 15 s, and annealing at 55 • C for 30 s.

cDNA Library Construction and Sequencing
Trizol was used for RNA extraction using the protocol of the manufacturer (Invitrogen, Carlsbad, CA, USA). Total RNA was extracted from wild type and MyoD1 knockout MDBK cells and mRNA was isolated using oligo dT magnetic beads (BioMag, Wuxi, China). The mRNA was fragmented into 200-300 nucleotide (nt) lengths using ion disruption (Shanghai Xinfan Technology Co., Ltd., Shanghai, China) and used as cDNA templates that were amplified using 6 nt random primers and reverse transcription (AXYGEN, Silicon Valley, USA). Second strand cDNA was synthesized following a similar procedure using the first-strand cDNA as template. These amplified fragments were used to construct the library that was sequenced by Nanjing Paterno Gene Technology (Nanjing, China) using the Illumina HiSeq platform (Illumina, San Diego, CA, USA) and tested for quality by an Aglent 2100 Bioanalyzer.

Enrichment of Differentially Expressed Genes
Clean reads were obtained by filtering the original sequencing data (Cutadapt) compared with the bovine genome (Bos_taurus.UMD3.1.dna.toplevel.fa) to obtain the gene sequence (Trimmomatic 0.39). Differentially expressed genes (DEG) were analyzed by the number of reads per 1000 bases of a gene per million reads (FPKM) and defined with DESeq software (false discovery rate (FDR) ≤ 0.01, absolute p value of Log2 ratio of ≥1). All genes were mapped to each term in the Gene Ontology (GO) database, and the number of differential genes in each term was calculated. Based on the whole genome, a hypergeometric distribution was used to calculate the terms and pathways of significantly enriched genes.

Validation of Differential Gene Expression
Heart, liver, spleen, lung, kidney, leg muscle, latissimus dorsi, and shoulder muscle tissues were collected from the Guanling cattle breeding farm in Guizhou province, and the amount of each tissue was about 0.25 g. The collected tissues were treated with PBS solution (10,000 U/mL penicillin and 10 mg/mL streptomycin) and put into 1.5 mL EP tubes, which were quickly transferred to the laboratory after being put into liquid nitrogen to extract tissue RNA using Trizol (Invitrogen, Carlsbad, CA, USA). After RT-PCR validation of the most significantly differential genes enriched, qPCR validation was then performed on the Guanling cattle tissue.
All animal experiments conformed to the Guizhou University Institutional Animal Care and Use Committee guidelines (Guizhou, China). The β-actin gene was used as a loading control reference gene (Table 3). All experiments were performed with 3 biological replicates and 3 technical replicates. The relative expression levels of the selected genes were calculated using the 2 −∆∆Ct method [14].  PDE1B is phosphodiesterase 1B, ADAMTS1 is a disintegrin and metalloproteinase with thrombospondin motifs 1, DPT is dermatopontin, CYP4F2 is cytochrome P450 family 4 subfamily F member 2, WNT7A is Wnt family member 7A, SEMA3A is semaphorin 3A, HOXD10 is homeobox D10, CCND2 is cyclin D2, and HR is Hairless.

Statistical Analysis
The data were analyzed using SSPS 18.0 software (IBM, Chicago, Ill, USA). Significant differences were determined using a one-way ANOVA F-test at p < 0.05. DESeq2 was used to analyze the difference in gene expression and the condition of screening DEGs was multiple of differential expression where |log Fold change| > 1, at a significance of p < 0.05. Volcano plots for DEGs were drawn using the R language ggplot2 software package.

MyoD1 Gene Knockouts
Our procedure for MyoD1 gene knockouts utilized a plasmid vector encoding puromycin resistance. Therefore, we conducted preliminary experiments of the selection procedure and added puromycin to wild type MDBK cells to evaluate the lowest concentration that could completely kill the cells. We found that 2 µg/mL resulted in almost complete cell death by 48 h ( Figure S1). We also examined whether Lipofectamine 2000 could be used for direct transfection of the CRISPR plasmids, and we utilized a commercial plasmid GFP vector to test this. However, the lipofectamine transfection efficiency of MDBK cells was <5%, so this direct avenue of gene introduction into these cells was abandoned ( Figure S2).
We subsequently utilized the lentivirus infection method and infected MDBK cells with the lentivirus constructs for 48 h and then analyzed the cells for evidence of MyoD1 gene interruption. We found that PCR amplification of genomic DNA samples from cells using four different targets generated~700 bp amplicons that indicated an overall knockout effect ( Figure S3). These amplicons were subsequently sequenced and the target for sgRNA6 generated the most unambiguous interruption of the MyoD1 gene ( Figure S4). We therefore expanded the gRNA6 cells and generated clones by plating at 1 cell per well in 96-well plates (limiting dilution analysis). We acquired a total of 20 clones after 2 weeks of culture. Target gene mutations were present in eight of the new cell lines, and we selected four obvious mutations to confirm the addition/deletion/mutation of each allele ( Figure S5). The deletions in the MyoD1 coding regions occurred in clones 3 and 22 (Figure 1). These results indicated that the MyoD1 gene was successfully knocked out in these clones obtained from limiting dilution analysis. These results were corroborated by measurements of MyoD1 steady state mRNA levels, and MyoD1 was significantly downregulated in three mutants compared with wild type cells (Figure 2). Animals 2021, 11, x FOR PEER REVIEW 6 of 12 gene interruption. We found that PCR amplification of genomic DNA samples from cells using four different targets generated ~700 bp amplicons that indicated an overall knockout effect ( Figure S3). These amplicons were subsequently sequenced and the target for sgRNA6 generated the most unambiguous interruption of the MyoD1 gene ( Figure S4). We therefore expanded the gRNA6 cells and generated clones by plating at 1 cell per well in 96-well plates (limiting dilution analysis). We acquired a total of 20 clones after 2 weeks of culture. Target gene mutations were present in eight of the new cell lines, and we selected four obvious mutations to confirm the addition/deletion/mutation of each allele ( Figure S5). The deletions in the MyoD1 coding regions occurred in clones 3 and 22 ( Figure  1). These results indicated that the MyoD1 gene was successfully knocked out in these clones obtained from limiting dilution analysis. These results were corroborated by measurements of MyoD1 steady state mRNA levels, and MyoD1 was significantly downregulated in three mutants compared with wild type cells (Figure 2).

Transcriptome Sequencing
We utilized the clone gRNA6 cells for genomic analysis and procured >40 million clean reads for both wild type and MyoD1 knockout MDBK cells, and high-quality read percentages were accessed by a comparison to the reference genome. The total number of sequences was >36 million and the alignment rate was 89% (Table 4). In addition, transcriptome analysis indicated >4000 new transcripts, 4600 generated by variable splicing and >36,000 SNPs in wild type MDBK cells and 6300 new transcripts, 4500 variable splicing, and 38,000 SNPs in MyoD1 knockout MDBK cells (Table 5). By using clustering analysis of all samples (Figure 3), we can infer that MyoD1 knockout was effectual, and some genes were used to perform difference and significance analysis (Pheatmap).

Transcriptome Sequencing
We utilized the clone gRNA6 cells for genomic analysis and procured >40 million clean reads for both wild type and MyoD1 knockout MDBK cells, and high-quality read percentages were accessed by a comparison to the reference genome. The total number of sequences was >36 million and the alignment rate was 89% (Table 4). In addition, transcriptome analysis indicated >4000 new transcripts, 4600 generated by variable splicing and >36,000 SNPs in wild type MDBK cells and 6300 new transcripts, 4500 variable splicing, and 38,000 SNPs in MyoD1 knockout MDBK cells (Table 5). By using clustering analysis of all samples (Figure 3), we can infer that MyoD1 knockout was effectual, and some genes were used to perform difference and significance analysis (Pheatmap).   The data were utilized to generate DEG populations, and we identified 723 total DEGs. These included 370 upregulated and 353 downregulated genes in the experimental group compared with the control group (Figure 4).

GO Functional Annotation of DEGs
The DEGs were mapped using GO analysis by calculating the number of DEGs in each term. Employing the whole genome as the background, a hypergeometric distribution was used to calculate the term with significant enriched DEGs. The enrichment results indicated 668 DEGs that were involved in cellular components (CC), molecular function (MF), and biological processes (BP). The BP group consisted primarily of cell differentiation throughout biological processes; it was exclusively focused on muscle cell differentiation and muscle structure development ( Figure 5 and Table S1).

GO Functional Annotation of DEGs
The DEGs were mapped using GO analysis by calculating the number of DEGs in each term. Employing the whole genome as the background, a hypergeometric distribution was used to calculate the term with significant enriched DEGs. The enrichment results indicated 668 DEGs that were involved in cellular components (CC), molecular function (MF), and biological processes (BP). The BP group consisted primarily of cell differentiation Animals 2022, 12, 2571 8 of 12 throughout biological processes; it was exclusively focused on muscle cell differentiation and muscle structure development ( Figure 5 and Table S1).

GO Functional Annotation of DEGs
The DEGs were mapped using GO analysis by calculating the number of DEGs in each term. Employing the whole genome as the background, a hypergeometric distribution was used to calculate the term with significant enriched DEGs. The enrichment results indicated 668 DEGs that were involved in cellular components (CC), molecular function (MF), and biological processes (BP). The BP group consisted primarily of cell differentiation throughout biological processes; it was exclusively focused on muscle cell differentiation and muscle structure development ( Figure 5 and Table S1).

KEGG Enrichment Analysis of DEGs
All transcripts from wild type and mutant were also annotated to the KEGG database and the number of DEGs at different levels of the KEGG pathway were then counted to determine the primary metabolic and signaling pathways for these DEGs. The results indicated that 25,644 genes were enriched to the KEGG pathway, and this included 1417 DEGs. An analysis of the 20 most significant pathways revealed that these differential genes were enriched in Pl-3-kinase and AKT, p53 signaling pathways ( Figure 6).

KEGG Enrichment Analysis of DEGs
All transcripts from wild type and mutant were also annotated to the KEGG database and the number of DEGs at different levels of the KEGG pathway were then counted to determine the primary metabolic and signaling pathways for these DEGs. The results indicated that 25,644 genes were enriched to the KEGG pathway, and this included 1417 DEGs. An analysis of the 20 most significant pathways revealed that these differential genes were enriched in Pl-3-kinase and AKT, p53 signaling pathways ( Figure 6).

DEG Verification
Finally, based on the results of highly abundant and significant differences in cell differentiation, ten genes that showed the most significant differences in expression were

DEG Verification
Finally, based on the results of highly abundant and significant differences in cell differentiation, ten genes that showed the most significant differences in expression were selected to confirm whether the steady state transcriptome data correlated with quantitative RT-PCR measurements. The expressions of PDE1B, Histone H3.1, ADAMTS1, DPT, and CYP4F2 were significantly upregulated ( Figure 7a) and WNT7A, SEMA3A, HOXD10, CCND2, and HR were significantly downregulated in the MyoD1 knockout group (Figure 7b). This is consistent with the cluster analysis results.

DEG Verification
Finally, based on the results of highly abundant and significant differences in cell differentiation, ten genes that showed the most significant differences in expression were selected to confirm whether the steady state transcriptome data correlated with quantitative RT-PCR measurements. The expressions of PDE1B, Histone H3.1, ADAMTS1, DPT, and CYP4F2 were significantly upregulated ( Figure 7a) and WNT7A, SEMA3A, HOXD10, CCND2, and HR were significantly downregulated in the MyoD1 knockout group ( Figure  7b). This is consistent with the cluster analysis results.

Discussion
Muscle growth and development are usually regulated by specific core genes and signal transduction pathways [15,16], and the MyoD gene family (MyoD1, MYF5, MyoG,  and MYF6) are key regulators that control the expression of specific proteins in muscle cell proliferation and differentiation [17][18][19][20]. MyoD1, as a member of its family, is a major transcriptional regulator of muscle-specific genes, and its activity is directly related to

Discussion
Muscle growth and development are usually regulated by specific core genes and signal transduction pathways [15,16], and the MyoD gene family (MyoD1, MYF5, MyoG, and MYF6) are key regulators that control the expression of specific proteins in muscle cell proliferation and differentiation [17][18][19][20]. MyoD1, as a member of its family, is a major transcriptional regulator of muscle-specific genes, and its activity is directly related to muscle growth. Although previous studies related to the promoter binding site of the MyoD1 gene in Guanling cattle have been reported [12,21,22], the mechanism of action of the bovine MyoD1 gene on muscle cells has rarely been explored. In view of the cell growth arrest caused by MyoD1 knockdown in muscle cells in previous experiments, in the present study, after knockdown of the MyoD1 gene in bovine kidney cells, we discovered that PDE1B, histone 3.1, ADAMTS1, DPT, and CYP4F2 were upregulated, implying they may be a compensatory pathway after MyoD1 knockdown. In contrast, WNT7A, SEMA3A, HOXD10, CCND2, and HR were significantly downregulated, so we inferred that they might be downstream genes of the MyoD1 gene.
To clarify the relationship between these genes and economic traits in Guanling cattle, we examined their expression in visceral tissues (heart, liver, spleen, lung, kidney) and skeletal muscles (leg muscles, dorsal longus, shoulder), respectively, and found that all 10 genes were expressed in the tissues, but PDE1B, ADAMTS1, DPT, and CCND2 were expressed at significantly higher levels in leg muscles, dorsal longus, and shoulder than several other genes. This is also consistent with previous reports that PDE1B is associated with fat thickness and muscle in the latissimus dorsi region [23], ADAMTS1 plays a regulatory role in muscle regeneration after injury [24], and DPT enhances cell adhesion, reduces cell proliferation, and promotes myogenic hyperplasia [25]. PDE1B, ADAMTS1, and DPT were elicited after knockdown of the MyoD1 gene. The reason may be that they have a compensatory relationship in the regulation of muscle formation. In contrast, CCND2 was significantly downregulated after MyoD1 knockdown, with expression at the lowest level among the 10 genes but the highest expression in leg muscle tissue, suggesting that it is a downstream gene of the MyoD1 gene that promotes muscle cell growth and differentiation. This is consistent with the report that CCND2 is a novel key regulator of terminal differentiation of muscle progenitor cells [26], and we speculate that MyoD1-CCND2 in promoting cell growth, proliferation, and metabolism through the Pl3K-AKt signaling pathway.
In summary, we suggest that CCND2 is a key downstream gene of MyoD1 and is closely associated with the involvement in muscle formation and differentiation in Guanling cattle through the Pl3K-AKt signaling pathway. This provides experimental data for subsequent studies on the regulatory mechanisms of muscle differentiation in Guanling cattle.

Conclusions
In this study, the MDBK cell line possessing a MyoD1 gene knockout was successfully constructed using the CRISPR/CAS9 technique, and differentially expressed genes were identified by transcriptome sequencing and screening. These DEGs were significantly enriched in GO terms and KEGG pathways like PDE1B, ADAMTS1, and CCND2, suggesting that MyoD1 affects muscle cell differentiation and muscle growth mainly by regulating the action of these pathways. The PCR results also confirmed this inference.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/ani12192571/s1, Figure S1: MDBK cell survival in the presence of the indicated concentrations of puromycin after 48 h, Figure S2: Lipofectamine transfection of GFP control plasmid pX458, Figure S3: Lipofectamine transfection of GFP control plasmid pX458, Figure S4: Analysis of gRNA target sequences, Figure S5: Gel electorophoresis of PCR amplicons generated from MyoD1 gene disruptions in cell lines identified using limiting dilution, Table S1: GO enrichment analysis of DEGs.