Establishment of Bactrian Camel Induced Pluripotent Stem Cells and Prediction of Their Unique Pluripotency Genes

Induced pluripotent stem cells (iPSCs) can differentiate into all types of cells and can be used in livestock for research on biological development, genetic breeding, and in vitro genetic resource conservation. The Bactrian camel is a large domestic animal that inhabits extreme environments and holds value in the treatment of various diseases and the development of the local economy. Therefore, we transferred four mouse genes (Oct4, Sox2, Klf4, and c-Myc) into Bactrian camel fetal fibroblasts (BCFFs) using retroviruses with a large host range to obtain Bactrian camel induced pluripotent stem cells (bciPSCs). They were comprehensively identified based on cell morphology, pluripotency gene and marker expression, chromosome number, transcriptome sequencing, and differentiation potential. The results showed the pluripotency of bciPSCs. However, unlike stem cells of other species, late formation of stem cell clones was observed; moreover, the immunofluorescence of SSEA1, SSEA3, and SSEA4 were positive, and teratoma formation took four months. These findings may be related to the extremely long gestation period and species specificity of Bactrian camels. By mining RNA sequence data, 85 potential unique pluripotent genes of Bactrian camels were predicted, which could be used as candidate genes for the production of bciPSC in the future. Among them, ASF1B, DTL, CDCA5, PROM1, CYTL1, NUP210, Epha3, and SYT13 are more attractive. In conclusion, we generated bciPSCs for the first time and obtained their transcriptome information, expanding the iPSC genetic information database and exploring the applicability of iPSCs in livestock. Our results can provide an experimental basis for Bactrian camel ESC establishment, developmental research, and genetic resource conservation.


Introduction
The first embryonic stem cells (ESCs) were generated about 40 years ago by culturing mouse inner cell mass (ICM) on feeder layers [1]. ESCs are non-specialized cells and can differentiate into any type [2]. During the next 25 years, these ESCs could only be obtained by separating ICMs from early embryos, but this method is not applicable to large livestock [3,4]. In 2006, researchers first used retroviruses to transfect four mouse genes (Oct4, Sox2, Klf4, and c-Myc) into mouse embryonic fibroblasts and obtained cells similar to ESCs, which were termed induced pluripotent stem cells (iPSCs). This technology provides a practical pathway to acquire pluripotent stem cells from large livestock. With the development of this technology, iPSCs have been successfully generated from humans, mice, sheep, goats, cattle, canines, rats, rabbits, pigs, horses, and monkeys [3][4][5][6][7][8][9][10][11][12][13]. With indepth research on iPSCs, other inducible genes were discovered, including L-Myc, N-Myc,

Acquisition and Passage of bciPSCs
In this study, only AP staining was performed before determining the optimal induction culture protocol. The general flow of our experiments is shown in the schematic diagram in Figure 1C. BCFFs were transduced with retroviruses bearing the pluripotency and the reporter genes ( Figure 1C), and the cells showed morphological changes after 7 days. The cells were cultured until iPSC clones appeared, which were identified using AP staining (Figures 1Da and S1D). After statistical analysis, the best scheme for obtaining and culturing bciPSCs was found to be introducing the four factors of OSKM simultaneously, laying the feeder layer, and adding Lif to E8 as the stem cell culture medium ( Figure  1E).
When the colonies grew sufficiently, they were dissociated and transferred to new plates coated with mitomycin C-treated MEFs (Figure 1 (Db, Dc)). The cells were passaged when they reached 80% confluency.

Acquisition and Passage of bciPSCs
In this study, only AP staining was performed before determining the optimal induction culture protocol. The general flow of our experiments is shown in the schematic diagram in Figure 1C. BCFFs were transduced with retroviruses bearing the pluripotency and the reporter genes ( Figure 1C), and the cells showed morphological changes after 7 days. The cells were cultured until iPSC clones appeared, which were identified using AP staining (Figure 1Da and Figure S1D). After statistical analysis, the best scheme for obtaining and culturing bciPSCs was found to be introducing the four factors of OSKM simultaneously, laying the feeder layer, and adding Lif to E8 as the stem cell culture medium ( Figure 1E).
When the colonies grew sufficiently, they were dissociated and transferred to new plates coated with mitomycin C-treated MEFs (Figure 1(Db,Dc)). The cells were passaged when they reached 80% confluency.

Pluripotency Gene Detection
SOX2 and NANOG expression was detected in seven bciPSC cell lines ( Figure 1F)

Telomerase Gene Expression in bciPSCs
Telomerase (TERT) gene expression was significantly increased in all three bciPSC lines ( Figure 1G).

Telomerase Gene Expression in bciPSCs
Telomerase (TERT) gene expression was significantly increased in all three bciPSC lines ( Figure 1G).

Immunocytochemistry
In the 11th-generation bciPSC-A3 cells, nine pluripotency-related proteins were selected for immunofluorescence detection, and all marker proteins were expressed ( Figure 2B).

EB Formation
We inoculated the three lines of bciPSCs ( Figure 3A) into Petri dishes and cultured them in differentiation medium for 7 days to obtain EBs (Figures Figure 3A and S1E). We then detected germ layer marker genes by PCR. We found that endoderm (FOXA2), mesoderm (MYOD1, DES, and MYF5), and ectoderm (GFAP, and PAX6) genes were expressed in the differentiated cells but not in undifferentiated cells ( Figure 3B). After differentiation and additional culture, immunofluorescence analysis revealed that EBs were positive for FOXA2 (endoderm), DES (mesoderm), GFAP (ectoderm), and β-tubulin (reference gene) ( Figure 3C). The results show that bciPSCs can differentiate into three-germ layer cells in vitro. detected germ layer marker genes by PCR. We found that endoderm (FOXA2), mesoderm (MYOD1, DES, and MYF5), and ectoderm (GFAP, and PAX6) genes were expressed in the differentiated cells but not in undifferentiated cells ( Figure 3B). After differentiation and additional culture, immunofluorescence analysis revealed that EBs were positive for FOXA2 (endoderm), DES (mesoderm), GFAP (ectoderm), and β-tubulin (reference gene) ( Figure 3C). The results show that bciPSCs can differentiate into three-germ layer cells in vitro.

RNA-Seq analysis
After receiving transcriptome data, 12 randomly selected genes were tested by qRT-PCR, and their log values were plotted. The results showed consistent trends between the qRT-PCR and the transcriptome, indicating that the sequencing results are accurate and reliable ( Figure S2A).
The correlation heat map and sample clustering diagram revealed high correlation and close affinity between the seven bciPSC lines ( Figure S2B,C). The violin chart demonstrated that the sequencing data had only a few abnormal values and was highly reliable ( Figure S2D). The upregulated and downregulated genes in each transcriptome dataset were mapped ( Figure 4A). The results indicated that the expression of a large number of genes changed during the reprogramming of BCFFs and that the reprogramming direction was similar.
We selected 7561 significantly (|log 2 FC| ≥ 1, p < 0.05) upregulated and downregulated genes from seven bciPSC lines to create a heatmap ( Figure 4B). To create a non-clustering heatmap, these differentially expressed genes were mapped to the microarray data of humans, mice, and pigs ( Figure 4B). A Venn diagram was constructed using the significantly differentially expressed genes of the four species ( Figure 4C). The results indicated that bciPSCs had 2056 (human), 408 (mouse), and 1317 (pig) genes in common, indicating that bciPSCs were more similar to human ESC and pig iPSCs. Therefore, during the follow-up analysis, we focused mainly on humans and pigs as supplemental datasets.
The top 20 GO term bubble diagrams of bciPSCs indicated that they mainly belong to cellular components ( Figure S3A). The enrichment circle maps of the signal pathway with a q value < 0.05 indicated that the pathway related to genetic information processing was activated ( Figure S3B). Next, we selected the signal pathways related to stem cell formation (apoptosis, oxidative phosphorylation, signaling pathways regulating pluripotency of stem cells, and PI3K, MAPK, Wnt, JAK-STAT, and TGF-β signal pathways [39][40][41][42][43][44][45]) in addition to the signal pathways with a q value < 0.05, and assembled them into a network diagram. The results showed that MAPK and cell cycle signal pathways were at the core ( Figure S4A). In the aforesaid studies, seven bciPSCs were associated with 52 signaling pathways, from which six studied pathways were selected (cell cycle, p53, apoptosis, cellular senescence, cancer, DNA replication, and oxidative phosphorylation signaling pathways [46,47]) to draw a heat map ( Figure 4D). The heat maps of DNA replication, apoptosis, and cell cycle signaling pathways [46] have been displayed separately ( Figure S4B-D). Overall, bciPSC and human ESC had remarkably similar gene expression patterns for these pathways.

Preliminary Mining of Potentially Pluripotency Genes of bciPSCs
In the direct analysis of sequencing data by screening the top 20 GO terms ( Figure S3A) and 52 signal pathways implicated in the pluripotency verification (2.10. RNA-Seq analysis) ( Figures S3B and S4A), 511 genes were assigned to group 1.
Based on the analysis of the reported inducible genes, 43 inducible genes were identified (Table S5), and 34 were detected in the sequencing data of bciPSCs. The distribution of 34 genes in the GO and KEGG databases is displayed in an alluvial plot (Figure 5Aa, Ab) to observe their distribution more intuitively. The primary enriched pathways in-cluded KEGG-A-class (cellular processes), KEGG-B-class (cellular community-eukaryotes), and pathway ko04550 (Table S6) (Table S6). After summarization, 2934 upregulated genes were obtained, of which 401 belonged to group 2.

Discussion
Although the production of and research on iPSCs is growing, it remains limited to human and model organisms. Establishing iPSC lines and genetic information databases for various non-model organisms is significant for studying iPSC mechanisms and the The venn diagram of the two data groups ( Figure 5B) revealed an insignificant difference between the groups, which therefore could be combined for analysis. According to the analysis of five network diagrams using the STRING database with a total number of proteins ≥6 ( Figure 5C), 46 proteins had over three protein connections ( Figure 5D). Eventually, 46 core genes were obtained. The 42 complementary genes were obtained by summarizing the first 30 genes of each bciPSC line ( Figure 5D). After comparing the core and complementary genes, a total of 85 genes were obtained, three of which were shared by both ( Figure 5D).

Discussion
Although the production of and research on iPSCs is growing, it remains limited to human and model organisms. Establishing iPSC lines and genetic information databases for various non-model organisms is significant for studying iPSC mechanisms and the resource development of non-model organisms. Based on previous reports on the establishment of iPSCs in other species [3][4][5][6][7][8][9][10], we selected the retroviral system with the broadest host range to import inducible genes. The comparison of different culturing methods revealed that trophoblasts and Lif can promote the generation of bciPSCs and maintain their stable passage [42,43]. In an optimal cultural environment, the reprogramming efficiency is approximately 0.036%. Compared with that of iPSCs of other species, the induction system of bciPSCs needs to be improved [6]. Some cells undergo morphological alterations on the seventh day after exogenous gene insertion; however, stem cell clones arise later, which may be species-specific.
We verified the pluripotency of bciPSCs using various methods. bciPSCs have a high nucleocytoplasmic ratio, show epithelioid cell morphology, express pluripotency genes, are AP-positive, and have normal diploid chromosomes (2n = 74) [38]. Moreover, TERT, a marker of ES cells in several species, was active [6,58,59]. The immunofluorescence results of SSEA1, SSEA3, SSEA4, TRA-1-60, and TRA-1-81 were positive, which were different from the results of other species, perhaps due to the difference in reprogramming in different species [4,60]. bciPSCs can differentiate into the three germ layers in vitro and in vivo, and only bciPSC-A3 can form teratomas, possibly related to the test procedure or differentiation ability [6,60,61]. The teratoma formation time was 4 months or even longer, and this phenomenon has also been observed during the identification of horse iPSCs [7], which may be related to the extremely long gestation period of the Bactrian camels.
The RNA-Seq results showed that bciPSCs and BCFFs had many distinct genes, and the overall modifications were similar to those of human ESC. We further analyzed the pluripotency-related signaling pathway and pluripotency-related genes and DNA methylases, and the results indicated that the bciPSCs are pluripotent. We also found significant differences among the stem cells of Bactrian camels, humans [52], pigs [6], mice [62], bovines [63], and other species, suggesting that the Bactrian camel may have unique pluripotency-related genes. Therefore, establishing iPSCs in non-model organisms is crucial for studying the induction mechanism of iPSCs.

Plasmid Preparation
The mouse inducible pluripotency factor-carrying plasmids pMXS-Oct4, pMXS-Sox2, pMXS-Klf4, and pMXS-cMyc and the retroviral packaging plasmids pCMV-VSV-G and pUMVC were stored in our laboratory. Using EGFP as a template, we designed primers (Table S1) that were cloned and ligated into the HindIII/BamHI site of the pMXS plasmid to generate the pMXS-EGFP reporter plasmid. We introduced a pair of bases (TA) into the original pMXS-Sox2 plasmid (Table S1), complemented the three missing amino acids, and constructed a redesigned pMXS-Sox2 plasmid for our investigations. All plasmids were re-sequenced to evaluate their integrity before use.

Cell Lines and Cell Culture Conditions
BCFFs were derived from 3-month-old Bactrian camel fetuses and grown in a tissue block adherent culture after digestion and separation [90]. The cells were maintained at 37 • C and 5% CO 2 in DMEM/F12 glucose (Hyclone, Logan, UT, USA) supplemented with 10% fetal bovine serum (FBS) (BI, Kibbutz Beit Haemek, Israel). The cells were subcultured once they reached 90% confluence. The first three generations of BCFFs were preserved in liquid nitrogen.
Mouse embryonic fibroblasts (MEFs) were derived from 13.5-day-old mouse embryos. The culture method and number of cryopreserved generations were the same as those of All experimental procedures were approved by the Animal Care and Use Committee of the College of Veterinary Medicine at Gansu Agricultural University.

bciPSC Culture Method Analysis
We used 293T cells to package five retroviruses, four of which contained inducible pluripotency genes (Oct4, Sox2, Klf4, and c-Myc), and one contained a reporter gene (EGFP), and transduced them into laboratory-preserved BCFFs. The culture medium was replaced with an iPSC induction medium composed of KnockOut DMEM/F12 medium (Gibco, Grand Island, NY, USA) supplemented with 20% KnockOut Serum Replacement (Gibco, Grand Island, NY, USA), 2 mmol non-essential amino acids (NEAAs) (Gibco, Grand Island, NY, USA), 0.1 mmol β-mercaptoethanol (Sigma, St. Louis, MO, USA), and 1 mM valproic acid (Sigma, St. Louis, MO, USA). The cells were subcultured after seven days. However, because of the lack of certainty in producing these cells, we evaluated six different subsequent culture protocols (Table S2) by assessing the number of ESC-like colonies obtained from retrovirus-transduced BCFFs on day 34. Based on the results, we selected the optimal induction culture protocol, which was used for all future experiments.

bciPSC Generation
We produced bciPSCs according to the previously obtained optimal culture protocol. On day 34, the clones of bciPSCs were selected using a mechanical method for subculture. The obtained bciPSCs were divided into two groups: group A (originally preserved in the laboratory) and group B (supplemented with a missing base pair) according to the different pMXS-Sox2 plasmids in the induction combination.

Pluripotency Gene Detection
The detection primers for pluripotency genes were designed using the NCBI Primer-BLAST tool (Table S3), and the transcripts of bciPSC lines were analyzed using quantitative real-time polymerase chain reaction (qRT-PCR). Based on the results, BCFFs and seven bciPSC lines were used for transcriptome sequencing. bciPSC-A2, bciPSC-A3, and bciPSC-B13 were randomly selected from seven bciPSC lines for the subsequent culture and verification.

qRT-PCR
Total RNA was extracted from cells using TRIzol reagent (Sigma, St. Louis, MO, USA) according to the manufacturer's protocol and stored at −80 • C. cDNA was prepared using the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer's protocol. Quanti Nova SYBR Green PCR Kit (Qiagen, Dusseldorf, Germany) was used for qRT-PCR assays. All qRT-PCR reactions were performed on a Veriti 96-well thermal cycler (Bio-Rad, Hercules, CA, USA). This method was used to detect endogenous genes (such as nSOX2 and nNANOG) (Table S3), the telomerase gene (TERT) (Table S3), and RNA-Seq results (Table S4). Specific primers were designed using the NCBI Primer-BLAST tool to amplify fragments corresponding to the selected genes. All samples were amplified in triplicate, and the mean and standard error values were calculated. Expression levels of all genes were calculated using the 2 −∆∆CT method and were normalized to β-actin.

Alkaline Phosphatase Staining
Alkaline phosphatase (AP) staining was performed using the Alkaline Phosphatase Detection Kit (Beyotime, Shanghai, China).

Karyotype Analysis
The cells were transferred to a T25 bottle. When their density reached approximately 70%, colchicine (Sigma, St. Louis, MO, USA) was added to the final concentration of 0.25 µg/mL; the cells were digested and collected with trypsin after a 2.5 h treatment. The cells were then resuspended in 8 mL of potassium chloride solution (37 • C; 0.075 M), allowed to stand at 37 • C for 30 min, and 1 mL of freshly prepared fixative (glacial acetic acid: methanol = 1:3) was added to collect the cells. A few drops of the fixative were slowly added and mixed well, and another 8 mL was added to collect the cells. The cells were resuspended in 0.3 mL of the fixative, and 2-3 drops of the suspension were dropped onto the 4 • C glass slide. The slide was blown onto and baked several times with an alcohol lamp. The slide was then treated with 1% Giemsa solution (Beyotime, Shanghai, China) dye solution for 10-15 min; the seal film was washed off, observed, and micrographed with an oil lens. As there was no standard banding mode map for Bactrian camels, only a statistical analysis of the number of chromosomes was performed. The number of three bciPSCs for mitotic-phase spreads assessment was 22 (bciPSCs-A2), 25 (bciPSCs-A3), and 20 (bciPSCs-B13), respectively.

Immunofluorescence
Immunofluorescence was performed using the same steps as in fibroblast identification, outlined above. Primary antibodies included anti-mouseOCT4 (

Teratoma Formation
To evaluate pluripotency in vivo, we used matrix glue (Thermo Fisher Scientific, Waltham, MA, USA) to attach BCFFs and three bciPSC lines (1 × 10 7 cells/site), which were then resuspended and transplanted into the dorsal flank of immune-deficient BALB/cNude mice. Three mice were used for each cell, for a total of twelve mice. Six-to eight-week-old nude mice were housed under pathogen-free conditions in a temperature-controlled room on a 12/12 h light/dark schedule with food and water ad libitum. After four months, the mice were sacrificed, and the tumors were dissected. Teratomas were then analyzed by hematoxylin and eosin (HE) staining (Solarbio, Beijing, China), a reticular fiber staining kit (Solarbio, Beijing, China), and collagen fiber staining (Solarbio, Beijing, China). Simultaneously, the three-germ layer identification antibody was used for the immunohistochemical identification of teratomas.

Heatmap Generation and Data Visualization
Since there is a dearth of information and research on Bactrian camel stem cells, we retrieved databases of microarray data from human [52], mouse [62], and pig [6] stem cells microarray data from Gene Expression Omnibus (GEO) DataSets (Human, GSE33298; Mouse, GSE51483; Pig, GSE15472) as references and used these for subsequent data analysis. We screened the microarray data of these three species according to the conditions of |log 2 FC| ≥ 1, p < 0.05.
The RNA-Seq sequencing results were verified using qRT-PCR and analyzed. We summarized the signal pathways with q < 0.05 and the first 20 items of Gene Ontology (GO) analysis used to observe the main changes in the genetic information. We compared the differentially expressed genes in the RNA-Seq data and the GEO dataset. We also analyzed the enriched signaling pathways, pluripotency-related genes, and important genes in the reprogramming process in RNA-Seq data. The online software of the sequencing company was used to visualize the sequencing data.

Preliminary Mining of Potentially Inducible Genes in bciPSCs
First, we used RNA-Seq sequencing data from bciPSCs to perform direct screening. We summarized upregulated genes involved in the signaling pathways identified during data visualization and the first 20 terms of the GO analysis in the seven bciPSCs. The upregulated genes were screened using counts ≤50 and ≥20 in BCFFs and bciPSCs, respectively, and labeled as group 1. Next, we performed screening based on known induced genes. Based on the relevant literature review, we compiled a comprehensive list of inducible genes across different species. Each bciPSC gene in the RNA-Seq data was annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (KEGG-A, KEGG-B, and pathways) and the GO databases (GO components, GO functions, and GO processes). Subsequently, we selected the main subdatasets containing the summarized inducible genes at the six annotation levels and screened the upregulated genes. The upregulated genes in the new data set were screened using counts of ≤50 in BCFFs and ≥20 in bciPSCs and designated as group 2.
Ven map analysis for groups 1 and 2 determined whether to combine the analyses based on the degree of dissimilarity between the gene sets. After confirmation, the screened genes were imported to the Search Tool for the Retrieval of Interacting Genes database(STRING, http://string.embl.de/ (accessed on 28 November 2021)), to both predict gene interactions and select core genes. The Bactrian camels were not included in the reference species list. Accordingly, the most closely related Bos taurus was selected as the reference species.
To compensate for the disadvantages of selecting groups 1 and 2, we screened the upregulated genes without any screening treatment using counts of ≤50 in BCFFs and ≥20 in bciPSCs and designated the gene set as group 3. The genes of each of the seven bciPSCs in group 3 were sorted in descending order of copy number. The first 30 genes were summarized and used as supplementary genes. Finally, the core genes and complementary genes were combined to form the final candidate gene set.

Conclusions
This is the first study to obtain bciPSCs with stable cell morphology, expression of pluripotency genes and markers, and the ability to differentiate in vitro and in vivo. Furthermore, we predicted 85 candidate genes that can be used for reprogramming Bactrian camel cells. Our efforts in transcriptome data mining and basic research support the de-velopment of the iPSC genetic information database and the use of iPSC in livestock by providing biological materials and data for establishing the Bactrian camel ESC, developmental research, genetic breeding, and the conservation of genetic resources. Considering the development of research on the disease treatment of Bactrian camels, this research is expected to contribute to the transformation of regenerative medicine in the future.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author on reasonable request.