Comprehensive Analysis of lncRNA and mRNA Reveals the Effect of ZBED6 on Spleen Growth in Pigs

: Transcription factor zinc-ﬁnger BED domain-containing protein 6 ( ZBED6 ) is unique to placental mammals and regulates insulin-like growth factor 2 ( IGF2 ) expression, which lead to muscle growth. However, the effect of ZBED6 on the growth of spleen is still elusive. In this study, we explored the regulation of ZBED6 on spleen growth, and the results showed ZBED6 knockout ( ZBED6 KO) pigs had heavier spleens than wild-type (WT) pigs. To analyze the mechanism of increased spleen weight in ZBED6 KO pigs, long noncoding RNAs (lncRNAs) and mRNAs in the spleen samples (WT: ZBED6 KO pigs = 3:3) were analyzed to identify differentially expressed lncRNAs (DE-lncRNAs) and genes (DEGs) based on the RNA sequencing (RNA-seq) method. Then, 142 DEGs and 82 DE-lncRNAs were obtained. The qRT-PCR results were consistent with those of the RNA-seq, indicating that the data were reliable. The heavier spleen weight of ZBED6 KO pigs coincided with the signiﬁcantly upregulated IGF2 mRNA. Functional enrichment analysis of DEGs showed enrichment mainly in myoﬁbril assembly and sarcomere. In addition, 252 cis-and 109 trans-acting target genes of 82 DE-lncRNAs were predicted. By conjoint analysis of lncRNA and mRNA revealed that IGF2 , DE-lnRNAs ( XLOC_113021 , XLOC_078852 , NONSUSG004057.1 , NONSUSG014354.1 , and NONSUSG009750.1 ), and their target gene ACTN2 may be the key candidate genes in promoting spleen growth in ZBED6 KO pigs. This study provides new directions to understand the global functions of ZBED6 and lncRNAs in spleen growth in pigs.

Zinc-finger BED domain-containing protein 6 (ZBED6) is a transcription regulatory factor [14] and mainly affects muscle growth by inhibiting insulin-like growth factor 2 (IGF2) [15][16][17][18].ZBED6-IGF2 binding is hindered by a mutation (G to A) of IGF2 which results in postnatally upregulated mRNA expression in muscle tissue (three-fold) and increased muscle weight in pigs [17].Recently, studies about transgenic mice (ZBED6 knockout and IGF2 mutation knockin) showed that these two transgenic mice had higher expression of IGF2, more skeletal muscle, and a bigger heart [19].Chromatin immunoprecipitation (ChIP) sequencing showed that ZBED6 has more than 1200 target genes including IGF2 in C2C12 cells, which suggested that ZBED6 is highly significantly associated with other biological processes through changing the expression of various important function genes [14,20,21].ZBED6 knockout Chinese Bama pigs (ZBED6 KO) were produced in our previous study to explore ZBED6 function [22].Using the ZBED6 KO pigs, we focused on how ZBED6 regulates skeletal muscle and internal organ development by RNA-seq and lncRNA.Then, we found that the knockout of ZBED6 markedly results in higher lean mass and heavier heart and liver tissues by regulating the expression of IGF2, lncRNA, and other genes in pigs [22][23][24].However, the effect of ZBED6 on the spleen and the specific function of mRNA and lncRNA in the spleen of ZBED6 KO pigs are poorly understood.In summary, it is necessary to analyze the effect of ZBED6 on spleen development and its molecular mechanism using the ZBED6 knockout Chinese Bama pig model.
In this study, ZBED6-KO pigs were used to study the regulatory effect and its molecular mechanisms of ZBED6 on spleen development.Spleen traits (founder and F4 generation) were measured to illustrate the effects of the regulation of ZBED6 on spleen growth.Then, mRNAs and lncRNAs of the spleen (WT:ZBED6 KO pigs = 3:3) were integrated and analyzed to investigate the molecular mechanism.Here, we report the change in spleen phenotype and the function of mRNAs and lncRNAs in spleen growth in ZBED6 KO pigs.

Slaughtering Experiment
ZBED6 knockout Chinese Bama pigs (ZBED6 KO) that we developed in previous research [22] were used to investigate the regulatory effect of ZBED6 on the spleen.Six female founders (WT:ZBED6 KO = 3:3, 8 months), eleven female F4 (WT:ZBED6 KO = 5:6, 5 months), and six male F4 (WT:ZBED6 KO = 3:3, 5 months) pigs were slaughtered after anesthetized, and then we obtained the spleen weight and spleen weight rate (spleen weight/body weight) to evaluate the effect on spleen caused by ZBED6 knockout.Spleen samples were collected for RNA extraction and paraffin sectioning.

Histological Analysis
Histological sectioning was carried out following the instruction manual that were used in this experiment.Spleen samples were prepared for histological sectioning Briefly, paraffin waxes were used to wrap spleen tissues, which were preserved in 4% paraformaldehyde solution.The paraffin waxes wrapping tissues were sectioned into 5 µm thick slices, and we performed HE staining.We used a microscope (Olympus Corporation, Kusatsu, Japan) was used to obtain images of the HE staining to observe the morphology.

Total RNA Extraction and RNA Sequencing
Using an RNeasy Mini kit (QIAGEN, Hilden, Germany) and an Agilent 2100 Bioanalyzer (USA), we isolated total RNA from the spleen samples of female founders (WT:ZBED6 KO = 3:3), and we detected the RNA concentration.RNA (RIN value > 7.0) met the requirements for library construction and RNA sequencing (Beijing, China).

RNA-Seq Analysis
Clean reads were generated after deleting ployN, adapters, and low-quality reads from the raw reads, and were mapped to pig genome assembly Sscrofa11.1 using TopHat with default parameters, and we generated BAM files.The fragments per kilobase of transcript per million fragments mapped (FPKM) values of genes were quantified using the BAM files and Cufflinks.DEGs (adjusted p-value ≤ 0.05 and log2|fold change|> 1) between the two groups pigs were detected by Cuffdiff.

lncRNAs Analysis
To obtain clean reads, we removed low-quality raw reads (adapter sequences, N sequences > 10%, and low-quality reads).The Q30 was calculated and used to assess the quality of clean reads.TopHat software was used to map the clean reads to Sscrofa11.1 and generate the BAM files.Transcripts were constructed by Cufflinks.Then, we analyzed lncRNA using the lncRNA calling protocol and CNCI.The RNAs with lengths more than 200 nt, had more than 2 exons, and had a CNCI score more than 0 were screened as lncRNAs.We detected differentially expressed lncRNAs (DE-lncRNAs) of the spleen between the two groups pig using Cuffdiff.The screening criteria of DE-lncRNAs were adjusted p-value ≤ 0.05 and log2|fold change|>1.

lncRNA Target Gene Prediction
A previous study suggested that lncRNAs have cis-regulatory effects on their neighboring target genes [23].According to this theory, the genes that range from 100 KD upstream to 100 KD downstream of lncRNAs are selected as potential cis-acting targets.Additionally, the genes that are trans-regulated by lncRNAs on other chromosomes are considered to be trans-acting targets.The Pearson correlation coefficients of the lncRNAs and trans-acting targets were calculated, and the screening criteria of trans target genes were correlation value > 0.05 and q-value ≤ 0.5.

Clustering and Correlation and GO Enrichment Analysis
The FPKM values of all genes and gmodels were used to analyze the heatmap in R. The volcano plot was generated with the DEGs and DE-lncRNAs using the pheatmap package in R. GO enrichment analysis was carried out using the g:Profiler web server (https://biit.cs.ut.ee/gprofiler/),URL(accessed on 25 September 2022).The screening criterion of significant GO terms was adjusted p-value ≤ 0.05.

Association Analysis of lncRNA and mRNA
We overlapped the cis-acting targets, trans-acting targets, and DEGs to identify common genes.Then, we analyzed the functions of these common genes to screen for candidate genes associated with growth and development.Then, based on the expression of candidate genes and their targeted lncRNAs, we speculated the regulatory roles of lncRNA on spleen growth in ZBED6 KO pigs.

qRT-PCR
Four mRNA and five lncRNA were selected to carry out qRT-PCR.A kit (PrimeScript RT, Vazyme, Nanjing, China) was used to synthesize single-strand cDNA from RNA.An SYBR Premix Ex Taq (Takara, Kusatsu, Japan), RT PCR instrument and 384-well reaction plates (ABI 7500, Thermo Scientific, Waltham, MA, USA) were used to carried out the qPCR reactions.Premier Primer 5.0 was used to design primer sequences (Supplementary Table S1).The comparative cycle threshold (2−∆∆Ct) was used to analyze relative gene expression.

Statistical Analysis
GraphPad Prism 7 (USA) was used to draw graphs of spleen phenotypes and qRT-PCR results.We analyzed statistical significance of these data using paired t-tests.The results are presented as the mean ± SEM (***, p < 0.001; **, p < 0.01; *, p < 0.05).

ZBED6 Knockout Promotes Pig Spleen Growth
The detailed spleen indicators are recorded in Supplementary Table S2.The results of these seventeen pigs (seven founder pigs and ten F4 pigs) showed that ZBED6 KO pigs had a body weight similar to that of WT pigs (Supplementary Table S2).It is surprising that female ZBED6 KO pigs had heavier spleens (founder: 62.5% and F4: 65.7%) and a higher spleen weight rate (founder: 10.4%; F4: 11.8%), whereas male F4 ZBED6 KO pigs did not differ in spleen weight or rate from WT pigs (Figure 1A,B).The results of HE staining showed that the morphological structure and inflammation of the spleens of ZBED6 KO pigs were no different from those of WT pigs (Figure 1C).These result suggest that ZBED6 KO female pigs have heavier spleens.
had a body weight similar to that of WT pigs (Supplementary Table S2).It is surprising that female ZBED6 KO pigs had heavier spleens (founder: 62.5% and F4: 65.7%) and a higher spleen weight rate (founder: 10.4%; F4: 11.8%), whereas male F4 ZBED6 KO pigs did not differ in spleen weight or rate from WT pigs (Figure 1A,B).The results of HE staining showed that the morphological structure and inflammation of the spleens of ZBED6 KO pigs were no different from those of WT pigs (Figure 1C).These result suggest that ZBED6 KO female pigs have heavier spleens.

Whole-Transcriptome Sequencing Results
We performed cDNA sequencing of six spleen samples (WT:ZBED6 KO = 3:3) (Supplementary Table S3) to study the effect of ZBED6 on spleen growth.We generated 40 million raw reads from each sample.After filtering the adapters, ployN, and low-quality reads, a total of 307760381 clean reads were obtained from all libraries with a clean ratio of greater than 90%.The Q30, which is used to evaluate the quality of clean reads, was above 86.865% in each sample, suggesting the good quality of the clean reads.The

Whole-Transcriptome Sequencing Results
We performed cDNA sequencing of six spleen samples (WT:ZBED6 KO = 3:3) (Supplementary Table S3) to study the effect of ZBED6 on spleen growth.We generated 40 million raw reads from each sample.After filtering the adapters, ployN, and low-quality reads, a total of 307760381 clean reads were obtained from all libraries with a clean ratio of greater than 90%.The Q30, which is used to evaluate the quality of clean reads, was above 86.865% in each sample, suggesting the good quality of the clean reads.The mapping rate of the clean reads was 70%~88.3%.This suggested that our RNA-seq data were appropriate for further analysis.

Identification and Characteristics of lncRNAs
We generated 210,858 transcripts in the spleen to identify lncRNAs using Cufflinks software.Ultimately, 12,405 lncRNAs containing 10,378 known lncRNAs and 2027 potentially novel lncRNAs were identified (Figure 2A).We found the lncRNAs were mainly longer than 2000 bp (Figure 2B), and most of the lncRNAs had more than two exons (Figure 2C).In addition, we also obtained 25,878 expressed mRNAs.Base on the FPKM of lncRNAs and mRNAs that were counted, we compared results of the two groups and found that lncRNAs transcript expressions were relatively lower than those of mRNAs (Figure 3A).These features of lncRNAs, in agreement with those of previous studies, suggested the reliability of our results.
Agriculture 2023, 13, x FOR PEER REVIEW 5 of 13 mapping rate of the clean reads was 70%~88.3%.This suggested that our RNA-seq data were appropriate for further analysis.

Identification and Characteristics of lncRNAs
We generated 210,858 transcripts in the spleen to identify lncRNAs using Cufflinks software.Ultimately, 12,405 lncRNAs containing 10,378 known lncRNAs and 2027 potentially novel lncRNAs were identified (Figure 2A).We found the lncRNAs were mainly longer than 2000 bp (Figure 2B), and most of the lncRNAs had more than two exons (Figure 2C).In addition, we also obtained 25,878 expressed mRNAs.Base on the FPKM of lncRNAs and mRNAs that were counted, we compared results of the two groups and found that lncRNAs transcript expressions were relatively lower than those of mRNAs (Figure 3A).These features of lncRNAs, in agreement with those of previous studies, suggested the reliability of our results.

Analysis of DEGs and DE-lncRNAs
Cufflinks was applied to measure the expression level of lncRNAs and mRNAs, and Cuffdiff was applied to evaluate differences in the lncRNAs and mRNAs between ZBED6 KO and WT pigs.Hierarchical cluster analysis was performed with the FPKM values of the expressed mRNAs and lncRNAs in spleen tissues.The results of the heatmap plot showed that all mRNA and lncRNA expression patterns of the individuals in the two groups were more similar (Figure 3B,C), which further suggested the accuracy and dependability of our data.The screening threshold of DEGs and DE-lncRNAs was set to an adjusted p-value ≤ 0.05 and fold change > 2, and the results are shown in Supplementary Tables S4 and S5.In total, 142 DEGs were identified (upregulated: 78; downregulated: 64) (Figure 3D), and 82 DE-lncRNAs (upregulated: 42, downregulated: 40) were found (Figure 3E).IGF2 is a crucial regulator of development.The spleen of ZBED6 KO pigs showed higher IGF2 levels (2.8-fold; WT:ZBED6 KO = 43.6:124)than WT pigs (Figure 3F), suggesting that ZBED6-IGF2 mainly plays a role in spleen development in pigs.

Analysis of DEGs and DE-lncRNAs
Cufflinks was applied to measure the expression level of lncRNAs and mRNAs, and Cuffdiff was applied to evaluate differences in the lncRNAs and mRNAs between ZBED6 KO and WT pigs.Hierarchical cluster analysis was performed with the FPKM values of the expressed mRNAs and lncRNAs in spleen tissues.The results of the heatmap plot showed that all mRNA and lncRNA expression patterns of the individuals in the two groups were more similar (Figure 3B  The g:Profiler web server was used to perform GO enrichment analysis of 142 DEGs.We found 27 GO terms (Supplementary Table S6, adjusted p-value ≤ 0.05); most of the DEGs were enriched in the growth pathways, such as myofibril assembly and sarcomere (Figure 4A).From the GO terms, we found that growth factors, actin alpha cardiac muscle 1 (ACTC1), actinin alpha 2 (ACTN2), and myosin heavy chain 7 (MYH7), were significantly enriched in the myofibril assembly and sarcomere pathways.These results demonstrated that these three growth factors may play important roles in spleen development in ZBED6 KO pigs (Figure 4 and Supplementary Table S6).
The g:Profiler web server was used to perform GO enrichment analysis of 142 DEGs.We found 27 GO terms (Supplementary Table S6, adjusted p-value ≤ 0.05); most the DEGs were enriched in the growth pathways, such as myofibril assembly and sarcomere (Figure 4A).From the GO terms, we found that growth factors, actin alpha cardiac muscle 1 (ACTC1), actinin alpha 2 (ACTN2), and myosin heavy chain 7 (MYH7), were significantly enriched in the myofibril assembly and sarcomere pathways.These results demonstrated that these three growth factors may play important roles in spleen development in ZBED6 KO pigs (Figure 4 and Supplementary Table S6).

Prediction of DE-lncRNA Target Genes
The regulatory networks of lncRNAs and their target genes (mRNAs) in spleen growth were constructed to explore the potential key molecules.The potential cis and trans targets of DE-lncRNAs were predicted (Figure 5A).We obtained 252 cis-acting target genes and 109 trans-acting targets (Supplementary Tables S7 and S8).Five GO pathways were enriched with the 361 potential target genes, including the negative regulation of cell activation and response to external stimulus (Figure 4B and Supplementary Table S9).

Prediction of DE-lncRNA Target Genes
The regulatory networks of lncRNAs and their target genes (mRNAs) in spleen growth were constructed to explore the potential key molecules.The potential cis and trans targets of DE-lncRNAs were predicted (Figure 5A).We obtained 252 cis-acting target genes and 109 trans-acting targets (Supplementary Tables S7 and S8).Five GO pathways were enriched with the 361 potential target genes, including the negative regulation of cell activation and response to external stimulus (Figure 4B and Supplementary Table S9).

Association Analysis of lncRNA and mRNA
To study the effect of lncRNAs on spleen growth, the DEGs and potential target genes of DE-lncRNAs were overlapped.Five common genes were identified between the DEGs and cis-acting target genes, and thirty-five common genes were identified between the DEGs and trans-acting target genes (Figure 5B).These 40 genes were enriched in five immune process pathways (Figure 5C and Supplementary Table S10).In the present study, the network of DE-lncRNAs and their targets were performed by Cytoscape.We found no growth GO pathways, but there were several growth factors in the potential target genes of DE-lncRNAs.Among these 40 common genes, ACTN2 is related to growth and development and was upregulated in the spleen of ZBED6 KO pigs (Figure 4A).ACTN2 was the trans-acting target gene of five DE-lncRNAs (two novel lncRNAs, XLOC_113021 and XLOC_078852; three known lncRNAs, NONSUSG004057.1, NONSUSG014354.1, and NONSUSG009750.1) (Figure 5D), which were all upregulated in the spleen of ZBED6 KO pigs.In summary, XLOC_113021, XLOC_078852, NONSUSG004057.1, NONSUSG014354.1, and NONSUSG009750.1 may be potential regulators of ACTN2 expression during spleen development in ZBED6 KO pigs.

Association Analysis of lncRNA and mRNA
To study the effect of lncRNAs on spleen growth, the DEGs and potential target genes of DE-lncRNAs were overlapped.Five common genes were identified between the DEGs and cis-acting target genes, and thirty-five common genes were identified between the DEGs and trans-acting target genes (Figure 5B).These 40 genes were enriched in five immune process pathways (Figure 5C and Supplementary Table S10).In the present study, the network of DE-lncRNAs and their targets were performed by Cytoscape.We found no growth GO pathways, but there were several growth factors in the potential target genes of DE-lncRNAs.Among these 40 common genes, ACTN2 is related to growth and development and was upregulated in the spleen of ZBED6 KO pigs (Figure 4A).ACTN2 was the trans-acting target gene of five DE-lncRNAs (two novel lncRNAs, XLOC_113021 and XLOC_078852; three known lncRNAs, NONSUSG004057.1, NONSUSG014354.1, and NONSUSG009750.1) (Figure 5D), which were all upregulated in the spleen of ZBED6 KO pigs.In summary, XLOC_113021, XLOC_078852, NONSUSG004057.1, NON-SUSG014354.1, and NONSUSG009750.1 may be potential regulators of ACTN2 expression during spleen development in ZBED6 KO pigs.

RT-qPCR Verification
We performed RT-qPCR to confirm the reliability of transcriptomic data.B2M and HMBS, the expressions of which were constant for all samples, were selected as the reference genes.The primers are listed in Supplementary Table S1.Four DEGs and five DE-lncRNAs were verified.For lncRNAs, the expressions of XLOC_113021, XLOC_078852, NONSUSG004057.1, NONSUSG014354.1, and NONSUSG009750.1 were higher in ZBED6 KO pigs than in WT pigs (Figure 6).Among the target genes regulated by DE-lncRNAs, the expression of ACTN2 was increased in ZBED6 KO pigs.The RT-qPCR results of the other seven DEGs, including growth factors IGF2, ACTC1, and MYH7, had similar expression trends as in the RNA-seq data (Figure 6), indicating our results were reliable.lncRNAs were verified.For lncRNAs, the expressions of XLOC_113021, XLOC_078852, NONSUSG004057.1, NONSUSG014354.1, and NONSUSG009750.1 were higher in ZBED6 KO pigs than in WT pigs (Figure 6).Among the target genes regulated by DE-lncRNAs, the expression of ACTN2 was increased in ZBED6 KO pigs.The RT-qPCR results of the other seven DEGs, including growth factors IGF2, ACTC1, and MYH7, had similar expression trends as in the RNA-seq data (Figure 6), indicating our results were reliable.

Discussion
ZBED6, as a transcriptional repressor that is highly conserved in placental mammals, affects skeletal muscle growth by regulating IGF2 expression [17,19,25,26].Recently, ZBED6 KO pigs were produced through the CRISPR/Cas9 technique and showed increased IGF2 expression, higher muscle mass, and heavier heart and liver [22].Surprisingly, heavier spleen tissues were found in founder ZBED6 KO pigs.The founder pigs were generated from surrogate sows that contained zygotes, whereas WT pigs were naturally mating.This may be the reason that led to the different phenotypes of WT and ZBED6 KO pigs.The spleen phenotypic data of 17 F4 pigs were considered to further verify that the differences between the two groups in spleen growth were due to the editing of ZBED6, which showed that the spleens of ZBED6 KO pigs actually grew because of the knockout of ZBED6.The results of the spleen phenotypes in female F4 pigs were similar to those of founder ZBED6 KO pigs.Interestingly, the effect on spleen growth was only reflected in female ZBED6 KO pigs, whereas we observed no difference in males, which is similar to the results for liver weight of ZBED6 KO pigs and mice [19,22].The action of testosterone of male pigs that are not castrated may mask the change resulting from ZBED6 knockout.
In this study, the upregulated IGF2 is in agreement with the overgrowth of spleen in ZBED6 KO pigs, suggesting that the ZBED6-IGF2 axis is important for spleen development.Previously, IGF2 was found to be significantly higher in most tissues except the liver after ZBED6 knockout [22], which was also reported in mice and cattle [18,19,27].These findings confirmed the reliability of our data again.Previously, ChIP-seq, mRNA-

Discussion
ZBED6, as a transcriptional repressor that is highly conserved in placental mammals, affects skeletal muscle growth by regulating IGF2 expression [17,19,25,26].Recently, ZBED6 KO pigs were produced through the CRISPR/Cas9 technique and showed increased IGF2 expression, higher muscle mass, and heavier heart and liver [22].Surprisingly, heavier spleen tissues were found in founder ZBED6 KO pigs.The founder pigs were generated from surrogate sows that contained zygotes, whereas WT pigs were naturally mating.This may be the reason that led to the different phenotypes of WT and ZBED6 KO pigs.The spleen phenotypic data of 17 F4 pigs were considered to further verify that the differences between the two groups in spleen growth were due to the editing of ZBED6, which showed that the spleens of ZBED6 KO pigs actually grew because of the knockout of ZBED6.The results of the spleen phenotypes in female F4 pigs were similar to those of founder ZBED6 KO pigs.Interestingly, the effect on spleen growth was only reflected in female ZBED6 KO pigs, whereas we observed no difference in males, which is similar to the results for liver weight of ZBED6 KO pigs and mice [19,22].The action of testosterone of male pigs that are not castrated may mask the change resulting from ZBED6 knockout.
In this study, the upregulated IGF2 is in agreement with the overgrowth of spleen in ZBED6 KO pigs, suggesting that the ZBED6-IGF2 axis is important for spleen development.Previously, IGF2 was found to be significantly higher in most tissues except the liver after ZBED6 knockout [22], which was also reported in mice and cattle [18,19,27].These findings confirmed the reliability of our data again.Previously, ChIP-seq, mRNA-seq, and lncRNA analysis of muscle, heart, and liver provided reliable evidence that the ZBED6-IGF2 axis is important for growth, and there were other ZBED6 targets related to growth [14,19,22,23,28].A large amount of lncRNAs were proved to be related to embryonic development [29,30], imprinting control [31], cell differentiation and development [32][33][34], and muscle development [35,36].This showed that lncRNAs can be found by RNA-seq analyses.With the rapid development of RNA-seq analyses, many studies on lncRNAs in humans and mice have been reported [3,32,[37][38][39], while studies targeting lncRNAs are still in their infancy in pigs.Therefore, to analyze other target genes of ZBED6 affecting spleen growth, we performed RNA-seq and lncRNA analyses.In summary, our study found some candidate genes that play important roles in spleen growth in ZBED6 KO pigs.
We generated more than 300 millions reads of RNA-seq data from the spleen.Using the data, we identified more than 25,000 mRNAs in the spleen.The lncRNAs had a length of more than 2000 bp, with two or more exons, and expression levels were lower than those of mRNAs.These characteristics are similar to those reported in humans [40], mice [41], sheep [33,35], and horses [42].The reason that 2027 potentially novel lncRNAs were identified in this study may be the differences in the types of pigs and the tissue sampled in these studies.Thus, we performed RT-qPCR validation to verify that the lncRNAs were really transcripts.In this study, with adjusted p-value ≤ 0.05 and fold change > 2, 142 DEGs including IGF2 and 82 DE-lncRNAs, were screened.These results suggested that ZBED6 can affect spleen growth of pig by regulating lncRNAs and mRNAs.The GO terms of DEGs showed that growth pathways form complex mechanisms that may regulate spleen development, such as myofibril assembly and sarcomere, which include three growth factors, ACTC1, ACTN2, and MYH7.In addition, the GO terms of the DEGs and potential target genes of DE-lncRNAs suggested that ZBED6 may affect the immune process in the spleen.Interestingly, ACTN2, which was upregulated in the spleen of ZBED6 KO pigs, was also included in the 361 potential target genes of DE-lncRNAs (252 cis-and 109 trans-acting targets).We identified five DE-lncRNAs, XLOC_113021, XLOC_078852, NONSUSG004057.1, NONSUSG014354.1, and NONSUSG009750.1 that coexpressed with ACTN2, which encodes sarcomeric a-actinin-2, which constitutes the Z-line in mammalian muscle fibers.Due to the central role of ACTN2 in development [43][44][45][46], XLOC_113021, XLOC_078852, NONSUSG004057.1, NONSUSG014354.1, and NONSUSG009750.1 may influence spleen growth in ZBED6 KO pigs by regulating the expression of ACTN2.

Conclusions
In this study, we found that founder and F4 female ZBED6 KO pigs had heavier spleen tissue than WT pigs.Using the RNA-seq data of spleen tissues (ZBED6 KO:WT = 3:3), a total of 142 DEGs and 82 DE-lncRNAs were obtained.Among these 142 DEGs, IGF2, as a ZBED6 target gene and a key molecule related to development, was upregulated 2.8-fold in ZBED6 KO pigs compared with in WT pigs.We found 27 DEGs were significantly enriched in GO terms, and two terms were related to growth and development, both of which included ACTC1.With 82 DE-lncRNAs, 252 cis-and 109 trans-acting targets were predicted, among which XLOC_113021, XLOC_078852, NONSUSG004057.1, NONSUSG014354.1, and NONSUSG009750.1 regulate trans-acting target gene ACTN2.We speculate that IGF2, the five lncRNAs (XLOC_113021, XLOC_078852, NONSUSG004057.1, NONSUSG014354.1 and NONSUSG009750.1), and their target gene ACTN2 are key candidate genes in promoting spleen growth in ZBED6 KO pigs.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/agriculture13010108/s1,Table S1: Details of primers used for qPCR analysis; Table S2: Spleen indexes between WT and ZBED6 KO pigs; Table S3: Quality and genome mapping analyses of transcriptome sequencing; Table S4: List of DEGs in spleens between WT and ZBED6 KO pigs; Table S5: List of DE-lncRNAs in spleens between WT and ZBED6 KO pigs; Table S6: GO terms of DEGs in spleens between WT and ZBED6 KO pigs; Table S7: Cis-target genes of DE-lncRNAs in spleen between WT and ZBED6 KO pigs; Table S8: Trans-target genes of DE-lncRNAs in spleens between WT and ZBED6 KO pigs; Table S9: GO terms of cis-and trans-target genes of DE-lncRNAs in spleen between WT and ZBED6 KO pigs; Table S10: Gene Ontology analysis of common genes between the DEGs and target genes of DE-lncRNAs in spleen between WT and ZBED6 KO pigs.

Figure 1 .
Figure 1.Spleen phenotypes of ZBED6 KO founder and F4 pigs.(A,B) The phenotypes of WT and ZBED6 KO in terms of spleen weight and ratio.Red points indicate actual phenotype data.(C) The results of spleen HE staining in WT and ZBED6 KO founder pigs.Scale bar = 200 μm; **, p < 0.01; *, p < 0.05.

Figure 1 .
Figure 1.Spleen phenotypes of ZBED6 KO founder and F4 pigs.(A,B) The phenotypes of WT and ZBED6 KO in terms of spleen weight and ratio.Red points indicate actual phenotype data.(C) The results of spleen HE staining in WT and ZBED6 KO founder pigs.Scale bar = 200 µm; **, p < 0.01; *, p < 0.05.
,C), which further suggested the accuracy and dependability of our data.The screening threshold of DEGs and DE-lncRNAs was set to an adjusted p-value ≤ 0.05 and fold change > 2, and the results are shown in Supplementary

Figure 4 .
Figure 4. Enrichment analysis of DEGs and DE-lncRNAs target genes.

Figure 4 .
Figure 4. Enrichment analysis of DEGs and DE-lncRNAs target genes.

Figure 5 .
Figure 5. lncRNAs in spleen development in ZBED6 KO pigs.(A) Number of cis-and trans-acting genes of lncRNA.(B) Venn diagram of DEGs, cis-and trans-acting genes.(C) Enriched gene pathway terms of 40 overlapping genes in Figure (D) Function network of ACTN2.Red indicates key gene ACTN2, green indicates known lncRNAs, and gray indicates potentially novel lncRNA.

Figure 5 .
Figure 5. lncRNAs in spleen development in ZBED6 KO pigs.(A) Number of cis-and trans-acting genes of lncRNA.(B) Venn diagram of DEGs, cis-and trans-acting genes.(C) Enriched gene pathway terms of 40 overlapping genes in Figure 5B.(D) Function network of ACTN2.Red indicates key gene ACTN2, green indicates known lncRNAs, and gray indicates potentially novel lncRNA.