A Comprehensive Analysis of the Myocardial Transcriptome in ZBED6-Knockout Bama Xiang Pigs

The ZBED6 gene is a transcription factor that regulates the expression of IGF2 and affects muscle growth and development. However, its effect on the growth and development of the heart is still unknown. Emerging evidence suggests that long noncoding RNAs (lncRNAs) can regulate genes at the epigenetic, transcriptional, and posttranscriptional levels and play an important role in the development of eukaryotes. To investigate the function of ZBED6 in the cardiac development of pigs, we constructed the expression profiles of mRNAs and lncRNAs in myocardial tissue obtained from Bama Xiang pigs in the ZBED6 knockout group (ZBED6-KO) and the wild-type group (ZBED6-WT). A total of 248 differentially expressed genes (DEGs) and 209 differentially expressed lncRNAs (DELs) were detected, and 105 potential cis target genes of DELs were identified. The functional annotation analysis based on the Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) databases revealed two GO items related to muscle development by the cis target genes of DELs. Moreover, IGF2 was the direct target gene of ZBED6 by ChIP-PCR experiment. Our results explored the mechanism and expression profile of mRNAs and lncRNAs of ZBED6 gene knockout on myocardium tissue development, mining the key candidate genes in that process like IGF2.


Introduction
Pork is the third primary source of animal protein for human consumption and is still the most important meat production and marketing variety in the world (FAO). China is one of the leading countries in pork consumption. According to FAO statistics, China's pork production ranked first in the world in 2020. Pork plays an important role in China's consumption and daily diet and a vital role in the national economy. With the development of society and improvements in people's economic levels, consumers are paying more attention to the quality of pork. The world's average supply of protein of animal origin (g/cap/day) has increased from 26.8 g in 2000 to 31.6 g today (FAO). At present, the fat content of pork in China is high, but the amount of muscle is insufficient. Certain problems have been found in Chinese pork breeds, such as a low lean meat rate and a lack of lean pig breeds to conduct breed improvements [1]. Moreover, domestic pig breeding started

Library Construction and Sequencing
According to the quality of the RNA, a total of 6 samples were sequenced. Using NEBNext Poly (A) mRNA Magnetic Isolation Module to separate poly(A) mRNA from total RNA. NEBNext ® Ultra™ RNA Library Prep Kit (Illumina Inc., San Diego, CA, USA) was used to prepare the Sequencing library of 2 µg total RNA, the kit produced a 275 bp fragment including the average size of the adapters. In accordance with the instructions for the TruSeq PE Cluster Kit v3-cBot-HS kit (PE-401-3001, Illumina Inc., San Diego, CA, USA), the cBot Cluster Generation system was used to cluster samples with index joints. After the cluster was formed, the above samples were sequenced on the Illumina Hiseq 2500 platform for 150 cycles in paired-end mode, and 150 bp paired-end reads were generated. The Fastq files were uploaded on the Sequence Read Archive database: https://www.ncbi.nlm.nih.gov/sra, accessed on 1 June 2021 with the accession number PRJNA663759.

RNA-Seq Analysis
The raw data in Fastq format were first processed by inhouse scripts. The original reads of Illumina sequencing were obtained by removing the reads containing adapters, the reads containing ployN and the low-quality reads from the original data. FastQC v0.11.4 software was used to filter the original reads, and the low-quality bases (Phred value < 20) at the end of the reads and the reads with lengths less than 25 bp were removed. The Q20, Q30 and GC contents of clean reads were calculated. HISAT2 (v2.0.4 by Daehwan Kim) software [27] with default parameters was used to compare clean data after quality control with the Ensemble Sscrofa11.1 reference genome: ftp://ftp.ensembl.org/pub/ release-101/fasta/sus_scrofa/DNA/Sus_scrofa.Sscrofa11.1.DNA.toplevel.fa.gz, accessed on 10 April 2022. StringTie software (v1.3.4d) carried out the reconstruction of transcripts. Together with HISAT2, it allows biologists to identify new genes and new spice variants of known ones [28,29].

Quantification of mRNA and lncRNA Expression
The quantification of the transcripts and genes was performed using RSEM [35] software, and reads per kilobase of transcript per million mapped reads (RPKM) were obtained. The EdgeR [36] R package was used for differential expression analysis.
The differentially expressed genes (DEGs) were scanned with p-value ≤ 0.05, fold change (FC) ≥ 1.5 and value of each sample ≥1. The differentially expressed lncRNAs (DELs) between ZBED6-KO group and ZBED6-WT group were identified with p-value ≤ 0.05, fold change (FC) ≥ 1.5 and value of each sample ≥0. The pipeline used to find DEGs and DELs between two groups from the RNA-seq data is shown in Supplementary Figure S1B.

lncRNA Target Genes Prediction
Among its roles, lncRNA can regulate the expression of mRNA by cis or trans. When the role of lncRNA is limited to the transcriptional chromosomes (especially neighboring genes), it is mainly regulated by cis, and when it affects the expression of genes on other chromosomes (remote), it plays the role of trans [37]. Based on the theory of the cis-acting regulatory element, the protein coding genes located within 100 kb from lncRNA were selected as potential cis-acting target; for trans-acting target prediction, the correlation coefficients between the coding genes and lncRNAs were calculated, which required the sample size to be more than five. Therefore, we did not predict the trans target genes of DELs in this study.

Cluster Analysis
The RPKM values for all of the annotated transcripts from the six samples transcriptomes were used to perform PCA, which was implemented with the "ggfortify" package in R (version 3.1.3, http://cran.r-project.org/, accessed on 18 April 2022). The hierarchical cluster analysis of DEGs was carried out using the cluster analysis function in the "pheatmap" package of R language.

Enrichment Analysis
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway Analyses of DEGs and DELs were carried out with the g:Profiler web server: http://biit.cs.ut.ee/gprofiler/, accessed on 23 April 2022 [38]. Significantly enriched GO terms or KEGG pathways were identified by Benjamini-Hochberg FDR, and p.adj < 0.05 was used as the threshold. The results were visualized using the R language package.

Chromatin Immunoprecipitation (ChIP)-PCR Analysis
ChIP was performed in native conditions. Briefly, freshly snap-frozen tissues were treated with 1% formaldehyde in medium for 18 min and neutralized with glycine (AM-RESCO, Radnor, PA, USA) for 5 min at room temperature. The tissues were smashed and resuspended in SDS lysis buffer (Beyotime, Shanghai, China). After incubation for 30 min at 4 • C, the lysates were sonicated 24 times (30 s each) (Bioruptor ® Sonication System, Denville, NJ, USA). An equal amount of chromatin was immunoprecipitated at 4 • C overnight with 8 ug of the following antibodies: anti-ZBED6 (Atlas Antibodies, Stockholm, Sweden, HPA068807) and normal mouse IgG (CST, 2729S). Immunoprecipitate were collected after incubation with Protein G agarose beads (ThermoFisher Scientific, Waltham, MA, USA). The beads were washed, and bound chromatin was eluted in ChIP Elution Buffer. Then, 20 microliters 5 M NaCl was added to the combined eluates and reverse histone-DNA crosslinks by heating at 65 • C for 4 h, after which RNase A (TIANGEN, Beijing, China) was added for 30 min at 37 • C; then the proteins were digested with Proteinase K (TIANGEN, Beijing, China) for 2 h at 65 • C. The coprecipitated DNAs were purified using a ChIP DNA Clean & Concentrator purification spin column (ZYMO, Irvine, CA, USA) and eluted in 20 µL dilution buffer. The ZBED6 binding site was evaluated using PCR and normalized with total chromatin (input). Normal mouse IgG was used as the negative control, and the primers are described in Supplementary Table S1. The PCR conditions were: 1 cycle at 95 • C for 5 min followed by 32 cycles at 95 • C for 30 s, Tm • C for 30 s, and 72 • C for 30 s. The PCR products were then electrophoresed on 1.5% agarose gels stained with GelGreen (TIANGEN, Beijing, China).
It has been shown that ZBED6 can bind to IGF2 through the "GCTCG" motif and then inhibit its expression in mammals [4,5]. Previous ChIP-seq using mouse C2C12 cells showed that the ZBED6 binding site of IGF2 corresponding to the porcine QTN site was among the most highly enriched regions [5]. Therefore, we blasted the IGF2 binding sequence in C2C12 with Sus scrofa IGF2 located on chromosome 2 and designed primers in the blast region for ChIP-PCR to examine the binding of IGF2 and ZBED6.

Quantitative Real-Time PCR (RT-qPCR) Analysis
Differentially expressed genes were randomly selected to identify the gene expression. In order to verify the stability of candidate genes in offspring, myocardial tissue samples of six three-month-old individuals of F3 generation pigs were collected, three in each group.
For the RT-qPCR analysis of genes, reverse transcription was performed using a PrimerScriptTM RT regent Kit (TaKaRa, Dalian, China) according to the manufacturer's instructions. RT-qPCR with TB Green ® Premix Ex Taq™ II (TaKaRa, Dalian, China) was performed with a QuantStudio TM 7 Flex 96-well System (Thermo Fisher Scientific, Waltham, MA, USA) as follows: 95 • C for 30 s; 40 cycles of 95 • C for 5 s, 60 • C for 34 s; and 95 • C for 15 s, 60 • C for 1 min, and 95 • C for 15 s. The data were analyzed with the 2 −∆∆Ct method [39] (Livak & Schmittgen 2001). Relative gene expression data were analyzed using real-time quantitative PCR and the 2(−Delta Delta C(T)) Method. The pig β-actin gene and GAPAH gene were used as reference genes for the normalization of the target gene data. The sequences of RT-qPCR primers are listed in Supplementary Table S1.

Statistical Analysis
Statistical analyses of the RT-qPCR results and graphs were carried out in GraphPad Prism (version 8.3) software (San Diego, CA, USA). The statistical significance of the data was tested by performing paired t-tests. The results are presented as means ± SEM of three replicates, and the statistical significance was represented as p-values < 0.05 or 0.01. Error bars show standard error of the mean, *** p < 0.001, ** p < 0.01, * p < 0.05.

Overview of RNA Sequencing
The number of raw reads generated from each library exceeded 100 million reads. After filtering low-quality reads, each sample was sequenced to obtain at least 3G data, the clean reads ranged from 3.9G to 5.2G, and the clean read quality scores of Q20 and Q30 were above 94.72% and 88.46%, respectively. The GC content of the clean data ranged from 47.77% to 50.70% (Supplementary Table S2). More than 81.2% of the clean reads mapped to the genome using Hisat (v2.0.4) software (Supplementary Table S3). Therefore, the reliability and quality of the sequencing data were adequate for further analysis.

Features of mRNAs and lncRNAs
The RPKM values of six samples in the two groups are shown in Supplementary  Table S4. There were marked differences between coding and noncoding at the gene structure and expression levels. The distribution of RPKM values of lncRNAs was more concentrated, while that of mRNAs was more extensive. Moreover, the expression of mRNA was generally higher than the value of lncRNA in each sample ( Figure 1A). After counting the numbers of mRNAs and lncRNAs in six samples, we found that the number of mRNAs with RPKM values ≥1 was significantly greater than those lncRNAs ( Figure 1B). We also counted the ORF length and exon number of the coding and noncoding genes, finding that the ORF of mRNA was longer than noncoding and the exon number of mRNA was much more than the lncRNAs (Supplementary Figure S2). centrated, while that of mRNAs was more extensive. Moreover, the expression of mRNA was generally higher than the value of lncRNA in each sample ( Figure 1A). After counting the numbers of mRNAs and lncRNAs in six samples, we found that the number of mRNAs with RPKM values ≥1 was significantly greater than those lncRNAs ( Figure 1B). We also counted the ORF length and exon number of the coding and noncoding genes, finding that the ORF of mRNA was longer than noncoding and the exon number of mRNA was much more than the lncRNAs (Supplementary Figure S2). Principal component analysis of all genes and lncRNAs was conducted using the R language package. Figure 2A shows that the mRNAs expressed in the three samples of WT and the samples of the ZBED6-KO group were clustered together, with PC1 accounting for 62.18%, indicating that the repeatability of each group was good. Similarly, the PCA of lncRNAs showed the same result, which significantly separated the two groups, with PC1 accounting for 100% ( Figure 3B). Although individual "Heart 4" deviated in the Principal component analysis of all genes and lncRNAs was conducted using the R language package. Figure 2A shows that the mRNAs expressed in the three samples of WT and the samples of the ZBED6-KO group were clustered together, with PC1 accounting for 62.18%, indicating that the repeatability of each group was good. Similarly, the PCA of lncRNAs showed the same result, which significantly separated the two groups, with PC1 accounting for 100% ( Figure 3B). Although individual "Heart 4" deviated in the group, the cluster analysis of the differential genes deviated in the group, the cluster analysis of the differential genes indicated that the data were accurate and dependable and could satisfy the subsequent difference analysis. group, the cluster analysis of the differential genes deviated in the group, the cluster analysis of the differential genes indicated that the data were accurate and dependable and could satisfy the subsequent difference analysis.   The heavier the red is, the higher the difference in gene expression is, and the heavier the blue is, the lower the difference in gene expression is.

Identification of lncRNAs
The protein-coding capability of all noncoding transcripts with CPC2, CNCI, Pfam and PLEK was predicted, and then the 13667 lncRNA was obtained by taking the intersection of four kinds of software results ( Figure 3A).
After predicting the coding functions of the lncRNAs, the identification rate of the novel lncRNAs was shown in Supplementary Figure S3. It was divided into three categories, among which "Antisense lncRNA" accounted for 73.73%, "Inter lncRNA" accounted for 26.72% and "Intro lncRNA" accounted for 0%.

Differentially Expressed mRNAs and lncRNAs
According to the screening criteria, a total of 248 DEGs and 209 DELs were screened (Supplementary Table S5). The DEGs included 163 upregulated and 85 downregulated DEGs ( Figure 2B). Among the 209 DELs, 66 were upregulated and 143 were downregulated in the ZBED6-KO group compared with the WT group, and the log2 ratio ranged from −6.3478 to 3.6507 ( Figure 3C). The heavier the red is, the higher the difference in gene expression is, and the heavier the blue is, the lower the difference in gene expression is.

Identification of lncRNAs
The protein-coding capability of all noncoding transcripts with CPC2, CNCI, Pfam and PLEK was predicted, and then the 13667 lncRNA was obtained by taking the intersection of four kinds of software results ( Figure 3A).
After predicting the coding functions of the lncRNAs, the identification rate of the novel lncRNAs was shown in Supplementary Figure S3. It was divided into three categories, among which "Antisense lncRNA" accounted for 73.73%, "Inter lncRNA" accounted for 26.72% and "Intro lncRNA" accounted for 0%.

Differentially Expressed mRNAs and lncRNAs
According to the screening criteria, a total of 248 DEGs and 209 DELs were screened (Supplementary Table S5). The DEGs included 163 upregulated and 85 downregulated DEGs ( Figure 2B). Among the 209 DELs, 66 were upregulated and 143 were downregulated in the ZBED6-KO group compared with the WT group, and the log2 ratio ranged from −6.3478 to 3.6507 ( Figure 3C).

Hierarchical Cluster Analysis
Using ggplot2 in the R software package and the RPKM value of six samples to perform the hierarchical cluster analysis of the DEGs and DELs, the results, presented using a direct heatmap plot ( Figures 2C and 3D), showed that all gene expression patterns including mRNA and lncRNA patterns were more similar for three individuals (×1, ×3, ×6) in the ZBED6-KO group and for three individuals (×2, ×4, ×5) in the ZBED6-WT group. It was further demonstrated that the data obtained by RNA-seq were accurate and dependable.

Prediction of lncRNA Target Genes
To investigate the function of lncRNAs, we predicted the potential target genes of lncRNAs in cis actions by searching for protein-coding genes both 100 kb upstream and downstream of the lncRNAs. A total of 105 potential related cis-acting target genes were obtained (Supplementary Table S6), among which ACTB, ACTC1, CDH15, HAND2, MYO16 and MYOCD genes are related to muscle growth and heart development and are located near TCONS_00033689, TCONS_00091379, TCONS_00016689, TCONS_00021032, TCONS_00005021 and TCONS_00068622, respectively, thereby suggesting that the cisacting genes and the neighboring lncRNAs could simultaneously play a role in growth and development of myocardium.

GO and KEGG Enrichment Analysis
Using a p.adj ≤ 0.05 as the screening criteria, significant enrichment was found in 112 GO terms and 4 KEGG pathways from DEGs with 90 biological process, 10 cellular component terms and 12 molecular function as well as in 3 biological processes from cis target genes of DELs ( Figure 4, Supplementary Table S7). Enrichment analysis of DEGs demonstrated that most of these genes were clustered into the immune process, such as three GO terms in different classes with lower p-values: defense response to other organisms, immunoglobulin complex, circulating and immunoglobulin receptor activity. Although no GO terms of DEGs enriched in muscle growth and development, there were several related genes, like MYO19 and TNNT1. Among the GO terms enriched by cis target genes of DELs, there were two items related to muscle development, skeletal muscle thin filament assembly and skeletal myofibril assembly, and ACTC1 was significantly enriched in these two terms. These findings demonstrate one of the roles of lncRNAs: the regulation of development through their cis actions on their neighboring protein-coding genes, e.g., regulating the ACTC1 protein during muscle development.

ZBED6 Target Genes Analysis by ChIP-PCR
We examined the binding of IGF2 and ZBED6 to verify the DNA in ChIP experiments. The blast rate of IGF2 binding sequence in C2C12 with Sus scrofa IGF2 located on chromosome 2 was as high as 86%, and the ChIP-PCR result showed that IGF2 could bind

ZBED6 Target Genes Analysis by ChIP-PCR
We examined the binding of IGF2 and ZBED6 to verify the DNA in ChIP experiments. The blast rate of IGF2 binding sequence in C2C12 with Sus scrofa IGF2 located on chromosome 2 was as high as 86%, and the ChIP-PCR result showed that IGF2 could bind with ZBED6 in the myocardium, which indicated that the DNA by ChIP was reliable and usable ( Figure 5C).

ZBED6 Target Genes Analysis by ChIP-PCR
We examined the binding of IGF2 and ZBED6 to verify the DNA in ChIP experiments. The blast rate of IGF2 binding sequence in C2C12 with Sus scrofa IGF2 located on chromosome 2 was as high as 86%, and the ChIP-PCR result showed that IGF2 could bind with ZBED6 in the myocardium, which indicated that the DNA by ChIP was reliable and usable ( Figure 5C).

Verification of Gene Expression Profiles with RT-qPCR
To confirm the expression patterns of the mRNA and lncRNA, ten DEGs and six DELs were selected for RT-qPCR validation. We selected b-actin and GAPDH as the reference genes and assumed that the expression of the reference genes was constant for all the samples. The DEG expression levels were determined and compared with the RNA-seq data, which showed similar patterns to the RNA-seq data ( Figure 5A), and five DELs provided the same result ( Figure 5B).
To investigate the expression of candidate genes in offspring, we also quantitatively verified IGF2 and TCONS_00005021 in the myocardium tissue of F3 generation pigs. The result showed that the expression of IGF2 was significantly upregulated and the TCONS_00005021 was downregulated after ZBED6 knockout with GAPDH as reference gene ( Figure S4).
Based on the above results, we hypothesized a working model in which ZBED6-IGF2 axis was the major influence on myocardial development by direct interaction with IGF2. Concurrently, ZBED6 also regulated the expression of TCONS_00005021 that regulate the expression of ACTC1 in cis-mode. It was another factor to promote the growth of the myocardium ( Figure 6). TCONS_00005021 was downregulated after ZBED6 knockout with GAPDH as reference gene ( Figure S4). Based on the above results, we hypothesized a working model in which ZBED6-IGF2 axis was the major influence on myocardial development by direct interaction with IGF2. Concurrently, ZBED6 also regulated the expression of TCONS_00005021 that regulate the expression of ACTC1 in cis-mode. It was another factor to promote the growth of the myocardium (Figure 6).

Discussion
ZBED6 can inhibit the expression of IGF2, reduce the deposition of subcutaneous fat, promote the growth of skeletal muscle [26] and also affect the growth and development of the myocardium [40]. By knocking in the same IGF2 intron 3-3072 site, ZBED6 improved meat production in Chinese Bama Xiang pigs [41]. We identified the core lncRNAs and mRNAs by performing comprehensive analysis of RNA-seq data from ZBED6 knockout and wild-type Bama Xiang pigs and validated their expression with RT-qPCR. Overall, our work uncovered some candidate genes that are involved in myocardial development in ZBED6-KO pigs.
By comparing the characteristics of mRNA and lncRNA, we found that the expression and exon number of lncRNA were less than that of mRNA, which was consistent with the results of Lagarde et al. [42]. In addition, after screening and statistics with "value one" as the threshold, we found that the gene number of lncRNA in six individuals was less than that of mRNA. This was because the low expression of lncRNA and the proportion of lncRNA in mammalian cells was significantly less than that of mRNA [43]. Working model: ZBED6 could directly influence the expression of IGF2 and TCONS_00005021. Also, TCONS_00005021 regulates the expression of ACTC1 by cis-mode.

Discussion
ZBED6 can inhibit the expression of IGF2, reduce the deposition of subcutaneous fat, promote the growth of skeletal muscle [26] and also affect the growth and development of the myocardium [40]. By knocking in the same IGF2 intron 3-3072 site, ZBED6 improved meat production in Chinese Bama Xiang pigs [41]. We identified the core lncRNAs and mRNAs by performing comprehensive analysis of RNA-seq data from ZBED6 knockout and wild-type Bama Xiang pigs and validated their expression with RT-qPCR. Overall, our work uncovered some candidate genes that are involved in myocardial development in ZBED6-KO pigs.
By comparing the characteristics of mRNA and lncRNA, we found that the expression and exon number of lncRNA were less than that of mRNA, which was consistent with the results of Lagarde et al. [42]. In addition, after screening and statistics with "value one" as the threshold, we found that the gene number of lncRNA in six individuals was less than that of mRNA. This was because the low expression of lncRNA and the proportion of lncRNA in mammalian cells was significantly less than that of mRNA [43].
In this paper, 248 DEGs and 209 DELs were screened by transcriptome analysis, which should have a close relationship with the knocked-out ZBED6 gene. After ZBED6 knockout, the IGV of RNA-seq results showed that IGF2 was significantly upregulated in myocardial tissue. Previous studies have confirmed with RT-qPCR that IGF2 expression is significantly higher in ZBED6-KO group than in the ZBED6-WT group [26], which is consistent with the results in mice and cattle [4,6,7]. We again confirmed the reliability of our data.
As the first organ formed in the process of vertebrate growth and development, the heart is mostly composed of myocardial tissue. Previous histological analysis of myocardial HE-stained sections showed that knockout of ZBED6 resulted in thicker myocardial muscle fibers and a reduction in connective tissue. These findings are complementary to the heart weight results, indicating that ZBED6 knockout could promote the growth and development of pig hearts to a certain extent [26].
In the enrichment analysis of cis target genes of DELs, we found that there were two items, skeletal myofibril assembly and skeletal muscle thin filament assembly, related to muscle development, and ACTC1 was significantly enriched in these two GO terms. ACTC1 (actin α cardiac muscle1) is an α-actin protein found in cardiac muscle and is the main component of thin-filament genes in mature cardiac myocytes, and it is completely conserved in mammals. Mutations in ACTC1 have been phenotypically related to various cardiac anomalies including dilated cardiomyopathy (DCM), restrictive cardiomyopathy (RCM), left ventricular non-compaction (LVNC), HCM, and congenital heart disease (CHD) [44][45][46]. Mice lacking ACTC1 do not survive more than two weeks, and knockdown of ACTC1 in chick embryos causes ASD [44,47]. Findings confirmed that ACTC1 is critical for normal cardiac morphogenesis and muscle contraction [48]. Based on our results, we believe that the expression of TCONS_00005021 was significantly downregulated after ZBED6 knockout and that this affects the expression of its cis target gene ACTC1. Although ACTC1 is not a differentially expressed gene, it changes its original transcriptional expression due to the regulation of lncRNA, which has a certain effect on the myocardium of Bama Xiang pig.
Although there were no significantly enriched GO terms or pathways of DEGs related to muscle or heart development, TNNT1 and MYO19 were significantly expressed after ZBED6 knockout.
The coordination of the TPM (tropomyosin)-TNN (troponin)-actin complex and myosin controls muscle contraction and relaxation [49,50]. TNNT (mammalian troponin T), TNNI and TNNC are troponin members. These proteins are important for regulating the contraction of myofibrils and maintaining the structural integrity of the myofibril ganglion [51]. Three homologous genes have evolved in vertebrates to encode three muscle type-specific TnT isoforms: TNNT1 for slow skeletal muscle TnT. Troponin T (TnT) is a central player in the calcium regulation of actin thin filament function and is essential for the contraction of striated muscles. Mutations in TNNT1 resulting in a complete loss of slow TnT in slow skeletal muscle also causes severe nemaline myopathy with childhood lethality [52,53]. The alternative splicing regulation of TNNT1 expression may play an important role in modulating muscle contractility in physiological and pathophysiological adaptations. Typical histological features of TNNT1 NM (Nemaline) include rods in type I and II fibers, increased fiber size variation, fiber type disproportion and selective type II fiber hypertrophy [54,55]. We think TNNT1 could play a certain role in muscle contraction: after ZBED6 knockout, the expression of TNNT1 was significantly upregulated and then promoted muscle contraction, thus affecting the growth and development of the myocardium.
Myosins are actin-based motors that power processes such as muscle contraction, cytokinesis and organelle transport [56]. Myosin-XIX (Myo19), a novel myosin, functions as an actin-based motor for mitochondrial movement in vertebrate cells [57]. MYO19 is broadly expressed in vertebrate cells, tissues and tumors. However, the mechanism of actin-based mitochondrial movements in vertebrates remains unclear [58]. The downregulation of MYO19 expression mediated by ZBED6 may provide a theoretical basis for the myocardial development and movement mechanism of actin in mitochondria.
ZBED6 could inhibit the expression of IGF2 through the combination of "GCTCG" motif. Using ChIP-PCR experiment, we first demonstrated that ZBED6 could bind with IGF2 in myocardial tissue through this binding motif, and we proved this point again in skeletal muscle and liver tissues of ZBED6-KO pigs [26], which is consistent with previous research results in long-white pigs and mice [5,59]. These findings indicated that in the myocardial tissue of Bama Xiang pig, the ZBED6-IGF2 axis still plays an important role in the growth and development of the myocardium. IGF2, as the main target of ZBED6, not only plays a role in skeletal muscle tissue but also has a significant effect on myocardial tissue. Therefore, we believe that knockout of ZBED6 increases heart weight by promoting the growth and development of the myocardium. At present, porcine heart has been successfully transplanted into baboon [60], and the studies have found the structural characteristics of the porcine heart valve [61]. In this study, we clarified the transcriptional expression profile in the myocardium of Bama Xiang pig, which could provide a theoretical basis and practical model for the study of heart organ transplantation and cardiovascular disease in pigs.

Conclusions
Through high-throughput RNA-seq, a total of 248 differentially expressed genes (DEGs) and 209 differentially expressed lncRNAs (DELs) were obtained between ZBED6-KO and ZBED6-WT pigs. The TNNT1 and MYO19 were both significantly expressed between groups, which had a certain effect on the development of muscle and myocardium. Based on the target gene prediction of DELs, a total of 105 potential related cis target genes were found. The DEGs were significantly enriched in 112 GO terms and 4 KEGG pathways. Additionally, there was significant enrichment in three GO terms of the cis target genes of DELs, among which two GO terms were related to muscle growth and development, and ACTC1 was significantly enriched in both two items. IGF2 was confirmed to be the target gene of ZBED6 in Bama Xiang pig by direct binding in ChIP-PCR experiments. We regard IGF2, TCONS_00005021 and its target gene-ACTC1 as key candidate genes in promoting myocardial growth and development in Bama Xiang pigs after ZBED6 gene knockout.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/genes13081382/s1, Figure S1: Overview of the experiment design, Figure S2: Features of mRNA and lncRNA, Figure S3: Classification of the novel lncRNAs, Figure S4: Verification of candidate genes in F3 generation pigs. Table S1: A list of primers, Table S2: The quality analysis of sequencing data, Table S3: The genome mapping of sequencing data, Table S4: A expression of all lncRNAs, Table S5: A list of all differentially expressed genes, Table S6: Potential related genes of differentially expressed lncRNA, Table S7: A enrichment analysis of DEGs.

Conflicts of Interest:
The authors declare no conflict of interest.