Key resources table

Reagent or resource

Source

Identifier

Antibodies

 Rabbit polyclonal anti-neuron specific beta III Tubulin

Abcam

Cat# ab18207, RRID: AB_444319

 Rabbit polyclonal anti-alpha smooth muscle Actin

Abcam

Cat# ab5694, RRID: AB_2223021

 Rabbit monoclonal anti-SOX17 (D1T8M)

CST

Cat#81,778, RRID: AB_2650582

 Rabbit monoclonal anti-GATA-6 (D61E4)

CST

Cat# 5851, RRID: AB_10705521

 Rabbit polyclonal anti-human Nanog

PeproTech

Cat# 500-P236, RRID: AB_1268274

 Mouse monoclonal anti-SSEA1 (MC480)

Abcam

Cat# ab16285, RRID: AB_870663

 Mouse monoclonal anti-SSEA4 (MC813)

Abcam

Cat# ab16287, RRID: AB_778073

 Mouse monoclonal anti-E Cadherin

Abcam

Cat# ab76055, RRID: AB_1310159

 Mouse monoclonal anti-Oct-3/4 (C-10)

Santa Cruz

Cat# sc-5279, RRID: AB_628051

 Mouse monoclonal anti-Sox-2 (E-4)

Santa Cruz

Cat# sc-365823, RRID: AB_10842165

 Monoclonal Anti-Skeletal Myosin (FAST)

Sigma-Aldrich

Cat# M4276, RRID: AB_477190

 Goat anti-Mouse IgG (H + L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 488

Thermo Fisher

Cat# A-11001, RRID: AB_2534069

 Goat anti-Mouse IgG (H + L) Cross-Adsorbed Secondary Antibody, Alexa Fluor 594

Thermo Fisher

Cat# A-11005, RRID: AB_141372

 Donkey anti-Rabbit IgG (H + L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor 594

Thermo Fisher

Cat# A-21207, RRID: AB_141637

 Donkey anti-Rabbit IgG (H + L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor 488

Thermo Fisher

Cat# A-21206, RRID: AB_2535792

Chemicals, peptides, and recombinant proteins

 CHIR99021

Selleckchem

Cat# S1263

 IWR-1-endo

Selleckchem

Cat# S7086

 Y-27632

Selleckchem

Cat# S1049

 WH-4-023

Selleckchem

Cat# S7565

 SB431542

Selleckchem

Cat# S1067

 PD0325901

Selleckchem

Cat# S1036

 XAV939

Selleckchem

Cat# S1180

 LDN193189

Stemgent

Cat# 04-0074

 Recombinant Human HGF

Peprotech

Cat# 100-39H

 Recombinant Human IGF-I

Peprotech

Cat# 100-11

 Recombinant Murine BMP-4

PeproTech

Cat# 315-27

 Recombinant human EGF

PeproTech

Cat# AF-100-15

 Recombinant human SGF

PeproTech

Cat# 300-07

 Human/Murine/Rat Activin A (E.Coli)

PeproTech

Cat# 120-14E

 Recombinant Human LIF

PeproTech

Cat# 300-05

 Recombinant Human FGF-basic (154 a.a.)

PeproTech

Cat# 100-18B

 Bovine Serum Albumin

Sigma-Aldrich

Cat# A1470

 Ascorbic Acid

Sigma-Aldrich

Cat# A4544

 KnockOut Serum Replacement

Thermo Fisher

Cat# A3181502

 Neurobasal™ Medium

Thermo Fisher

Cat# 21,103-049

 DMEM/F-12, GlutaMAX™ supplement

Thermo Fisher

Cat# 10,565-018

 N-2 Supplement (100 ×)

Thermo Fisher

Cat# 17,502-048

 B-27™ Supplement (50 ×), minus vitamin A

Thermo Fisher

Cat# 12,587-010

 GlutaMAX™ Supplement

Thermo Fisher

Cat# 35,050-061

 MEM Non-Essential Amino Acids Solution (100 ×)

Thermo Fisher

Cat# 11,140-050

 2-Mercaptoethanol

Thermo Fisher

Cat# 21,985-023

 Penicillin–Streptomycin (10,000 U/mL)

Thermo Fisher

Cat# 15,140–122

 Gelatin (0.1% in water)

Stem Cell Technologies

Cat# 07903

 Trypsin–EDTA (0.05%), phenol red

Gibco

Cat# 25,300,120

 DMEM, high glucose, no glutamine

Gibco

Cat# 11,960-044

 Fetal bovine serum (FBS)

Gibco

Cat# 16,000-044

 Accutase cell dissociation reagent

Gibco

Cat# A11105-01

 TrypLE ™ Express

Gibco

Cat# 12,605,010

 Dulbecco's phosphate-buffered saline (DPBS)

Gibco

Cat# C14190500CP

Critical commercial assays

 KAPA Hyper Prep Kits

KAPA Biosystems

Cat# KK8054

 RNeasy Mini Kit

QIAGEN

Cat# 74,106

 Globin-Zero Gold rRNA Removal Kit

Illumina

Cat# GZG1224

 NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina®

New England Biolabs

Cat# E7760

 RNAprep pure Cell/Bacteria Kit

TIANGEN

Cat# DP430

 5 × All-In-One RT Master Mix

ABM

Cat# G490

 Rapid Giemsa staining kit

BBI Life Sciences

Cat# E202FA0001

 Gel & PCR Clean Up Kit

OMEGA

Cat# D2000-02

 Lipofectamine ™ 3000 Transfection Kit

Invitrogen

Cat# 2,173,193

 Endo-free Plasmid Mini Kit

OMEGA

Cat# D6950-02B

 2 × RealStar Green Power Mixture

GenStar

Cat# A311-05

 Cells-to-cDNA™ II Kit

Invitrogen

Cat# AM8723

Deposited data

 scRNA-seq original data

This paper

GEO: GSE256201

 RNA-seq original data

This paper

GEO: GSE256201

 scRNA-seq public data

Zhi et al. [20]

GSA: CRA003960

 bEPSCs and biPSCs RNA-seq public data

Zhao et al. [18]

GEO: GSE129760

 Primed bPSCs RNA-seq public data

Bogliotti et al. [49]

GEO: GSE110036

 bEDSCs RNA-seq public data

Kinoshita et al. [16]

GEO: GSE172420

Experimental models: Cell lines

 E7-bEpiSCs

This paper

N/A

 E10-bEpiSCs

This paper

N/A

 E12-bEpiSCs

This paper

N/A

 E14-bEpiSCs

This paper

N/A

Experimental models: Organisms/strains

 CD-1® (ICR) IGS Mice

Vital River

201

 BALB/c Nude Mice

Vital River

401

 Holstein cow

Shijiazhuang Tianquan Elite Dairy Ltd and Beijing Dairy Cattle Center

N/A

Oligonucleotides

 See Table S8 for oligonucleotide sequences used in this paper

N/A

N/A

Recombinant DNA

 PB-CMV-EF1A-GFP-NLS

Our Laboratory

N/A

Software and algorithms

 kallisto v-0.46.0

Bray et al. [50]

https://pachterlab.github.io/kallisto/

 tximport v-1.18.0

Soneson et al. [51]

https://bioconductor.org/packages/release/bioc/html/tximport.html

 Genome and Genome annotation

Ensembl 105

https://dec2021.archive.ensembl.org/index.html

 Seurat v-4.1.0

Hao et al. [52]

https://satijalab.org/seurat/

 Monocle3 v-1.3.1

Cao et al. [53]

https://cole-trapnell-lab.github.io/monocle3

 destiny v-2.14.0

Angerer et al. [54]

https://bioconductor.org/packages/release/bioc/html/destiny.html

 KEGGREST v-1.30.1

Tenenbaum [55]

https://bioconductor.org/packages/release/bioc/html/KEGGREST.html

 gghalves v-0.1.1

N/A

https://github.com/erocoar/gghalves

 DESeq2 v-1.30.1

Love et al. [56]

https://bioconductor.org/packages/release/bioc/html/DESeq2.html

 ComplexHeatmap v-2.6.2

Gu et al. [57]

https://github.com/jokergoo/ComplexHeatmap

 FactoMineR v-2.4

Lê et al. [58]

http://factominer.free.fr/

 R v-4.0.5

N/A

https://www.R-project.org/

Experimental model and study participant detailsMice

CD-1® (ICR) IGS and BALB/c nude mice were purchased from Beijing Vital River Laboratory Animal Technology Co., Ltd. (Beijing, China) and were used for mouse embryonic fibroblasts (MEFs) isolation, and the bEpiSCs teratoma formation test. The MEFs were treated with mitomycin C (Selleckchem, S8146) to prepare feeder cells for bEpiSCs.

Bovine

Holstein cow embryos at E5, E6, E7, E10, E12, and E14 were used for bovine embryonic single-cell collection and single-cell transcriptome analyses. Embryonic day n (E(n)) embryos were obtained n days after mating. For bEpiSCs derivation, E7, E10, E12, and E14 embryos were used.

Collection of cow embryos and isolation of embryonic single cells

The embryonic cells utilized in this study were exclusively derived from Holstein cow embryos. E5-E7 embryos were obtained by thawing and culturing frozen cow embryos until the desired developmental stage was reached. The zona pellucida was subjected to a 15–30 s Pronase (Sigma, 10,165,921,001) treatment, followed by cleaning and removal using a solution composed of DPBS + 0.1% BSA. Subsequently, the embryonic cells were mechanically separated and transferred into the lysate. E10, E12, and E14 embryos were obtained through in vivo transplantation and subsequent flushing, with single embryonic cells being isolated and collected via a combination of enzymatic treatment and mechanical manipulation.

Bovine epiblast stem cells culture medium

The components of the bovine epiblast stem cells culture system were the same as the 3i/LAF culture system previously reported [20], and the contents of different components were optimized accordingly. Specifically, the basal medium (BM) was a 1:1 mixture of DMEM/F12 medium (Thermo Fisher Scientific, 10,565–018) and Neurobasal medium (Thermo Fisher Scientific, 21,103–049), and added 1 × N2 supplement (Thermo Fisher Scientific, 17,502–048), 1 × B27 supplement (Thermo Fisher Scientific, 12,587–010), 0.5% GlutaMAX (Thermo Fisher Scientific, 35,050–061), 1% nonessential amino acids (Thermo Fisher Scientific, 11,140–050), 0.1 mM β-mercaptoethanol (Thermo Fisher Scientific, 21,985–023), 1% penicillin–streptomycin (Thermo Fisher Scientific, 15,140–122), 5% knockout serum replacement (KOSR, Thermo Fisher Scientific, A3181502, optional), and 50 μg/mL ascorbic acid (Sigma–Aldrich, A4544). To prepare the 3i/LAF culture system, it is also necessary to add the following small molecule inhibitors or cytokines to BM: CHIR99021 (1 μM, Selleckchem, S1263), IWR-1-endo (0.5 μM, Selleckchem, S7086), WH-4–023 (1 μM, Selleckchem, S7565), recombinant human LIF (10 ng/mL, PeproTech, 300–05), recombinant human Activin A (25 ng/mL, PeproTech, 120-14E), and recombinant human FGF-basic (154 aa) (10 ng/mL, PeproTech, 100-18B). Optimized, the ROCK inhibitor Y-27632 (5 μM, Selleckchem, S1049) was added. In addition, XAV939 (1 μM, Selleckchem, S1180) can be used instead of IWR-1-endo. Unless otherwise noted, the base medium and supplement factor concentrations of the condition medium (e.g. I/F medium) described herein are the same as those of the 3i/LAF medium. The bEpiSCs were cultured with mitomycin C-treated mouse fibroblasts as feeder cells.

Derivation of bEpiSCs from bovine ICMs, epiblasts, and ectoderms

Bovine ICMs, epiblasts, and ectoderms were isolated by mechanical isolation and treated with TrypLE™ Express (Gibco, 12,605,010) for 3 min, then seeded onto feeder cells supplemented with 3i/LAF medium. Cultures were incubated at 5% O2 and 5% CO2 at 37 °C. The domed outgrowths were collected and digested using Accutase cell dissociation reagent (Gibco, A11105-01), and then passaged every 3 days at a ratio of 1:4.

Method detailsSingle-cell RNA library preparation and sequencing

Single-cell RNA-seq library was prepared using a modified Smart-seq2 protocol, as previously described in studies [59, 60]. Briefly, individual embryonic cells were transferred into a lysis buffer containing an 8 bp barcode. Subsequently, first-strand cDNA was reverse-synthesized and amplified in a reverse transcription (RT) mixture consisting of 4 U RNase inhibitor, 100 U SuperScript II reverse transcriptase (Invitrogen, 18,064,071), 1 mM dNTPs (TAKARA, 4019), 60 mM MgCl2, and 3 µM RT primer with 10 µM TSO primer. Following PCR amplification, the product was purified using 0.8 × Beckman’s AMPure XP beads (A63882). Biotin PCR enrichment was then performed to further enhance the library quality. Finally, construction of the single-cell RNA-seq library followed the instructions provided by KAPA Hyper Prep Kits with PCR Library Amplification/Illumina series (KAPA KK8054). High-quality libraries were sequenced using Illumina HiSeq Xten platform (Novogene) with paired-end reads of length 150 bp. The primers utilized in these experiments are listed in key resources table.

Cell population doubling time

The growth curve of bEpiSCs was plotted by seeding them at a density of 3 × 105 in 12-well plates. Cells were subsequently digested and counted at intervals of 12 h, 24 h, 36 h, 48 h, and 60 h using the Luna™ Automated Cell Counter. Three replicates were performed for each time point. The Doubling time was calculated using the following formula: Doubling time (DT) = 12 × [lg2/(lgNt-lgN0)], where 12 is the cell culture time (hours); Nt is the number of cells cultured at 48 h; and N0 is the number of cells recorded at 36 h.

Analysis of single-cell cloning efficiency

Bovine EpiSCs were digested with TrypLE™ Express (Gibco, 12,605,010) and filtered through a 40 μm cell filter. Cells were seeded in 6-well plates with 100, 500, and 1000 cells, respectively. After 3 days of culture, the number of colonies formed was counted, and the single-cell colony formation rate was calculated and averaged.

Karyotype analyses

Before performing karyotype analyses, 1% KaryoMAX Colcemid Solution (Gibco, 15,212,012) was added to the bEpiSCs culture medium and incubated for 1 h. The bEpiSCs were then dissociated into single cells using TrypLE™ Express (Gibco, 12,605,010) and collected through centrifugation. Subsequently, the bEpiSCs were suspended in a hypotonic solution of 0.075 M KCl (Sigma, P5405) and incubated at 37 °C for 15 min. Following this step, the bEpiSCs were fixed with methanol and acetic acid in a ratio of 3:1; this process was repeated three times. The resulting suspension of bEpiSCs was dropped onto a pre-cooled slide, thoroughly dried at room temperature, and stained with a 10% Giemsa Stain Solution (Sangon, E6073140001) for a duration of 30 min. More than 45 metaphase cells were examined for each cell line.

Alkaline phosphatase (AP) staining

For details about AP staining procedures and precautions, refer to Alkaline Phosphatase Detection Kit (Millipore, SCR004).

Immunofluorescence analyses

The cells were washed with DPBS for immunofluorescence (IF) analyses, followed by fixation in 4% paraformaldehyde (PFA) at room temperature for 30 min. Subsequently, the cells were rinsed with DPBS and permeabilized with 0.1% Triton X-100 for 20 min. After another round of DPBS washing, the cells were blocked with 3% BSA at room temperature for 1 h. The primary antibodies were incubated overnight at 4 °C and subsequently washed three times using a wash buffer (DPBS containing 0.1% Triton X-100 and 0.1% Tween 20). Secondary antibodies were incubated at room temperature for 1 h and then washed three times using the same wash buffer. Finally, DAPI staining was performed to visualize the nucleus, enabling direct observation and photography.

Embryoid body differentiation

Following digestion, the bEpiSCs were seeded at a density of 1 × 106 cells per well and cultured for 5–7 days on 35 mm low-attachment plates in MEF culture medium: DMEM (Gibco, 11,960–044) supplemented with 10% FBS (Gibco, 16,000–044), 1% penicillin–streptomycin (Thermo Fisher Scientific, 15,140–122), and 1% GlutaMAX (Thermo Fisher Scientific, 35,050–061), on a horizontal shaker at a speed of 70 rpm. Subsequently, Embryoid Bodies (EB) were transferred to a 12-well plate and incubated in the same medium for one week with bi-daily medium changes. Adherent cells were utilized for subsequent immunofluorescence staining.

Teratoma formation

The prepared 1 × 107 bEpiSCs were collected by centrifugation after digestion, resuspended in 50 μL BM, and subsequently injected subcutaneously into the neck of BALB/c nude mice. After a period of 4–5 weeks, the teratoma was collected and subjected to subsequent H&E analyses.

H&E analyses

The teratomas were subjected to two washes in DPBS and fixed with 4% PFA for a duration of 2 days at a temperature of 4 °C. Subsequently, the teratoma tissues underwent dehydration using an alcohol gradient (70%, 80%, 90%, 95%, and finally, 100% for each hour) before being transferred into xylene and embedded in paraffin. The samples were sliced to a thickness of 5 μm, followed by deparaffinization in xylene and rehydration using decreasing concentrations of ethanol. Finally, the samples were stained with haematoxylin (Sigma–Aldrich, MHS16) and eosin (Sigma–Aldrich, HT110116), after which they were observed under a microscope (Leica, DM5500B).

RT–qPCR

The total RNA was extracted from bEpiSCs using the RNA prep Pure Cell/Bacteria Kit (TIANGEN, DP430), followed by reverse transcription to cDNA using the 5 × All-In-One RT Master Mix (Abm, G490). Subsequently, PCR amplification was performed with the 2 × RealStar Green Power Mixture (GenStar, A311-05) on an Archimed Real Time System (ROCGENE). The data were analyzed utilizing the comparative CT (2−ΔΔCT) method. ΔCT values were calculated employing GAPDH as an internal control. Three independent biological replicates were conducted for all experiments. The primer sequences used in quantitative real-time PCR can be found in key resources table.

Myogenic differentiation of bEpiSCs

The protocol was performed as described for pig pgEpiSCs [19]. In brief, myogenic differentiation basic medium (MDBM) consisted of DMEM/F12 supplemented with 1% nonessential amino acids, 0.1 mM β-mercaptoethanol, 1% penicillin–streptomycin, 15% KOSR, and 200 μM ascorbic acid. For stage 1, the bEpiSCs were dissociated into small pieces, and cultured in MDBM supplemented with 1% B27 supplement, 3 μM CHIR99021, and 2 μM SB431542 (Selleckchem, S1067) maintained for 3 days. In stage 2, the culture medium was replaced with a combination of 3 µM CHIR99021, 2 μM SB431542, 0.5 μM LDN193189 (Stemgent, 04–0074), and 20 ng/mL recombinant human FGF-basic (154 a.a.) from day 4 to day 6. For stage 3, the culture was replaced with 10 ng/mL HGF (Peprotech, 100-39H), 10 ng/mL IGF-1 (Peprotech, 100–11), 20 ng/mL recombinant human FGF-basic (154 a.a.), and 0.5 μM LDN193189 for 2 days. During stage 4, the bEpiSCs undergo differentiate into myoprogenitors. The bEpiSCs-MPCs were treated with a concentration of 10 ng/mL IGF-1 for 4 days, and for stage 5, a combination of 10 ng/mL HGF and 10 ng/mL IGF-1 for 20–25 days. For skeletal muscle maturation, the cells were treated with N2 medium consisted by DMEM/F12 supplemented with 15% KOSR, 1% N2 supplement, 1% penicillin–streptomycin, and 1% nonessential amino acids.

Primordial germ cell like cells (PGCLCs) induction

The induction of PGCLCs was conducted according to a modified protocol [8, 45]. Briefly, 3000 bEpiSC cells were seeded into each well of a low-cell binding U-bottom 96-well plate (Beyotime, FULA962) and incubated for 24 h in 3i/LAF medium. Subsequently, the cells were switched to PGCLC induction medium consisting of BM medium supplemented with sodium pyruvate (1 mM, Thermo Fisher Scientific), BMP4 (200 ng/mL, PeproTech), recombinant human LIF (10 ng/mL, PeproTech), recombinant human SCF (100 ng/mL; PeproTech), recombinant human EGF (50 ng/mL; Peprotech), and Y-27632 (10 μM, Selleckchem) for a duration of three days.

rRNA-depleted RNA-seq

Total RNA was extracted separately from four bEpiSCs samples using the RNeasy Mini Kit (Qiagen, 74,106). For constructing strand-specific RNA-seq libraries, we employed a rRNA depletion protocol (Globin-Zero Gold rRNA Removal Kit, Illumina, GZG1224) in conjunction with the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, E7420S). All libraries were quantified using the Qubit dsDNA High Sensitivity Assay Kit (Invitrogen, Life Technologies, Q32851) and sequenced on an Illumina HiSeq 4000 platform.

Vector construction

The GFP plasmid was stored in our laboratory. In brief, we generated a PB-CMV-EF1A-GFP-NLS plasmid by modifying the PB-CAG-MCS vector (provided by Prof. Sen Wu). Specifically, we substituted the chicken β-actin promoter with the Homo sapiens elongation factor 1 alpha (EF1A) promoter and integrated a GFP-NLS cassette downstream of the EF1A promoter.

bEpiSCs transfection

Bovine EpiSCs were transfected using Lipofectamine™ 3000 Reagent (Invitrogen, L3000008). Specifically, transfection was performed 16 h after normal cell passage in a 24-well plate. First, 25 μL Opti-MEM™ Medium was added to a centrifuge tube followed by the addition of 0.75 μL Lipofectamine™ 3000 Reagent and mixing well. Then, another centrifuge tube was taken and filled with 25 μL Opti-MEM™ Medium to which 0.5 μg DNA and 1μL P3000™ Reagent were added and mixed well. The DNA mixture was then combined with the Lipofectamine™ reagent at a ratio of 1:1 and allowed to incubate for10-15 min before being added to the cells for culture.

Generation of bEpiSCs cloned embryos

Ovaries were obtained from cattle farms around Beijing, and oocytes in 3–8 mm follicles were extracted. The oocytes with three layers of cumulus cells were transferred to Maturation medium and matured at 38.5 °C and 5% CO2 for 18 h. Maturation medium was based on TCM199 (Gibco, 12,340–030), supplemented with 10% FBS (Gibco, 16,000–044), and 0.01 IU/mL follicle stimulating hormone (FSH, Sigma, F4021), 0.01 IU/mL luteinizing hormone (LH, Sigma, L6420), 1 μg/mL estradiol (Sigma, E2257). After 18 h of oocytes maturation, 0.1% hyaluronidase (Sigma, H4272) was used to remove the excess cumulus, the oocytes with polar bodies were placed in HM medium containing 7.5 μg/mL of cytorelaxin (Sigma, C6762) for 10 min, and then transferred to HM medium containing 10% FBS for enucleation. The enucleated oocytes were transferred into the Maturation medium until nucleus injection. bEpiSCs were differentiated in basal medium containing 10 ng/mL BMP4 (PeproTech, 315–27), 5 μM SB431542, and 10 ng/mL FGF2 for more than 1 week and then used as donor cells for nuclear transfer. Clear round donor cells were selected and injected into the periovoid space, and the cells were as close to the cytoplasm as possible to improve the fusion efficiency. The reconstructed embryos were placed between two electrodes of the fusion cell and aligned with a microneedle so that the somatic cells were oriented toward one of the two electrodes. The fusion conditions were as following: double DC pulse 2.5 kV/cm, 10 μs, 1 s interval; The fusion solution was 0.3 mM mannitol (Sigma, 1,375,105), 0.15 mM CaCl2 (Sigma, C7902), 0.15 mM MgCl2 (Sigma, M2393), and 0.05% BSA. The fusion rate was detected under stereomicroscope. The reconstructed embryos were then transferred to IVC medium contained 5 μM ionomycin (Sigma, 407,950) for 4 min and then transferred to IVC medium contained 2 mM 6-DMAP (Sigma, D2629) for 4 h of culture. The activated embryos were washed three times in IVC medium and transferred into IVC medium for culture.

Quantification and statistical analysisSingle-cell RNA-seq low level processing and filtering

For the STRT-seq dataset, Raw reads were split by 8 bp cell barcodes located on Read 2 allowing 2 mismatches. Additionally, the 8 bp unique molecular identifiers (UMIs) located on Read 2 were switched to the identifier line of paired Read 1. Then Read 1 was processed to remove the template switch oligo (TSO) primer, low quality bases and polyA sequence [59]. The trimmed reads were aligned to the respective reference genomes (bovine: Bos_taurus.ARS_UCD1.2; pig: Sscrofa11.1), respectively. The unique molecular identifier (UMI) counts using kallisto (v-0.46.0) [61].

scRNA-seq data analyses and cell-type identification

To pre-process our bovine data, we loaded the entire raw count matrix dataset with metadata using Seurat (v-4.1.0.9) [52]. The expression matrix was filtered based on the following criteria to get high quality cells and detectable genes: 1) the number of genes should be more than 2 000; 2) the proportions of UMIs for external RNA controls (ERCC) and mitochondrial genes should be below 40% and 50%, re respectively; 3) genes should be detected in at least 3 out of 728 sequenced cells. After filtering, 594 (486 embryonic cells, 108 pgEpiSCs cells) out of 728 cells were retained in the high-quality expression matrix for down-stream clustering and cell type identification. For the public datasets, the original information in the references were used.

Sequencing reads for each gene were normalized and scaled in each cell after regression on the mitochondrial gene ratio and cell cycle scores by “SCTransform” function. Recalibrated data was then used for principal components analysis (PCA) using “RunPCA” with default parameters. Dimensionality reduction was applied to the first thirty principal components, and the clusters were identified by construction of a shared nearest neighbor (SNN) graph. Finally, we visualized cell clusters after projecting them into two-dimensional latent space via t-distributed Stochastic Neighbor Embedding (t-SNE) by “RunTSNE” and Uniform Manifold Approximation and Projection (UMAP) by “RunUMAP” function. For bovine each embryonic day’s cell types, we identified the cell types with known marker genes. We filtered out some cells for which cell type could not be defined and retained 456 cells for further analyses (Table S1). The UMAP embeddings and cluster assignments computed in this section were later used in bovine embryonic dataset pseudotime analysis.

Identification of differentially expressed genes between embryonic stages in respective lineages

Based on the differentiation process during embryonic development, we grouped the bovine embryonic cells into three main lineages: embryonic lineages, covering pre-ICM at E5, ICM at E6 and E7, epiblasts at E10 and E12, and the ectoderms at E14; TE lineages, including pre-TE at E5, TE at E6, E7, E10, E12, and E14; and hypoblast lineages, including hypoblasts at E10, E12, and definitive endoderm at E14 (Fig. 1C).

Construction of expression tendencies

To trace the dynamic changes of DEGs during embryonic development, we constructed the expression tendencies of DEGs in epiblast lineages. We first calculated the average expression levels for each gene at specific embryonic development time points in each respective lineage separately. The average expression levels were rescaled for embryonic lineage and analyzed by the k-means clustering method with parameters k = 10 and iter.max = 100, grouping the DEGs with similar tendencies during embryonic development into separate clusters.

Cross-species comparative analyses of scRNAseq datasets

We downloaded orthologous genes list for three species (human, cattle and pig) via BioMart tool in the ensembl genome browser 105 (http://dec2021.archive.ensembl.org/index.html) and retained 16,841 genes, all of which were 1:1 orthologs. We then retained 1:1 orthologs in the bovine and porcine embryonic lineages datasets. For integrative analyses between porcine and bovine embryonic dataset, we used Seurat CCA approach for anchors determination and datasets alignment. Top 2000 genes that were repeatedly variable across datasets were selected and anchors were identified using “FindIntegrationAnchors” function with following parameters: “reduction = ‘cca’, k.anchor = 5, normalization.method = ‘SCT’”. The datasets were then integrated based on the identified anchors with “IntegrateData” function with following parameters: “dims = 1:30, k.weight = 50, normalization.method = ‘SCT’”. The genes average expression values in the integrated assay were obtained using the “AverageExpression” function, and the Spearman’s correlation coefficient between different cell types of two species was calculated using “cor” function. The fold change of genes expression level between naïve and formative states, formative and primed states of the two species were determined with the Wilcoxon rank-sum test by ‘‘FindMarkers’’ function, respectively. Only those with | ‘avg_logFC’|> 0.25 and ‘p_val’ < 0.05 were considered as DEGs.

Pseudotime analysis

R package Monocle3 (v-1.3.1) [53] was used to reconstruct the embryonic cells developmental trajectory. The UMI matrix was used as input and variable genes obtained by Seurat to sort cells in pseudotime. Individual cells of two species were further ordered with the destiny (v-2.14.0) R package [54] using the DEGs between porcine all embryonic day's cells that were calculated by “FindAllMarkers” function to determine the developmental order between different species.

RNA-seq data processing and analyses

Expression levels of protein coding genes in all samples (our original bulk RNAseq dataset, and public datasets of biPSCs, bEPSCs, primed bESCs and pgEpiSCs) were quantified as transcripts per million (TPM) using kallisto (v-0.46.0). DEGs between different cell types were identified using the DEseq2 tool (v-1.30.1) [56]. We used Benjamini–Hochberg adjusted false discovery rate (FDR) < 0.05 and absolute log2 (fold change) > 2 as cut-offs for statistical significance.

Transcriptome relatedness analyses of bEpiSCs to embryonic cells

For integrative analyses between bEpiSCs and bovine embryonic dataset, we used Seurat RPCA approach for anchors determination and datasets alignment. Top 2000 genes that were repeatedly variable across datasets were selected and anchors were identified using “FindIntegrationAnchors” function with following parameters: “reduction = ‘rpca’, k.anchor = 15”. The datasets were then integrated based on the identified anchors with “IntegrateData” function with following parameters: “dims = 1:30, k.weight = 50”. Scaling and PCA were applied on the combined datasets and the integrated PCA coordinates were used as the input of the clustering and t-SNE visualization workflow with the following parameter: ‘dims = 1:5’.

Functional enrichment analyses

Functional enrichment analyses of selected genes were performed using Metascape (http://metascape.org) The bovine genes were mapped to their human orthologs and human (Homo sapiens) was the target species for analyses. Enrichment analyses were performed using all genes in the genome as the background set with Gene Ontology (GO)-biological processes (GO-BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway as ontology sources. Terms with a minimum count ≥ 3, adjusted P < 0.01, and enrichment factor ≥ 1.5 were considered to be significant and similar terms were grouped into clusters. Key pathways were depicted using -log10 (P-value) bar plots.

Statistical analysis

Two-way ANOVA multiple comparisons test was used to analyses data of RT–qPCR in Figs. 5C–H, 7C, F and Figure S3B, S3D, and used to analyses data presented in Figs. 4D, E, 5B. Dunnett’s multiple comparisons test were utilized in Fig. 5B, D, E, and Figure S3B, S3D. Bonferroni’s multiple comparisons test were utilized in Fig. 5C, F–H, and Fig. 7C, F. Tukey’s multiple comparisons test were utilized to compare the Cell population doubling time, and single-cell cloning efficiency in Fig. 4D, E.