Experimental animals

Animals used in this study were regularly maintained in the Laboratory Animal Centre of Yunnan Agricultural University, Kunming, China. All animal experiments were performed with the approval of the Animal Care and Use Committee of the above institution. The female pigs overexpressing leptin were obtained by backcrossing the WT female with already generated leptin overexpressing transgenic boar [18] and offspring were screened for the presence of transgenic vector containing leptin and enhanced green fluorescent protein (EGFP) by polymerase chain reaction (PCR) using primers listed in additional file Table S1, and circulating serum leptin concentrations. Furthermore, the information regarding animal used in this study is provided in additional file Table S2.

DNA extraction and polymerase chain reaction (PCR)

The genomic DNA from blood of each newborn piglet was extracted using TIANamp genomic DNA kit (Tiangen, Beijing, China) and PCR was carried out using 2×Ex Taq Master Mix (TaKaRa CW0718). The PCR reaction was composed of 0.5 µg of template, 10 µl of Ex Taq Master Mix, 0.4 µM of each primer and distilled water to a final volume of 20 µl. The amplification of each primer pairs was performed in separate reaction with 35 cycles of 94 °C for 30 s, 55–60 °C for 30 s, and 72 °C for 30 s. PCR products were separated by 1% agarose gel electrophoresis and visualized by staining with DL 2000 DNA Marker (TSINGKE). Depending upon the presence of transgenic vector comprised of leptin gene with EGFP expression, the offspring were divided into two groups; Leptin pigs and WT pigs. Next, the circulating serum leptin concentrations were measured using Porcine Leptin Elisa Kit (SEKP-0278, Solarbio life sciences, China).

Body weight and size measurements

Body weight (BW) and size including length, height, heart girth, chest depth, chest width, and abdominal circumference of leptin (n = 5) and WT (n = 4) pigs were measured at 30d intervals from d180 ~ d390 of age. The live weight of each animal was determined by suspending the animal on an electronic scale balance and the weight of each animal was taken and recorded. The body size measurements were recorded. In brief, body length was measured using the tape rule as the distance from the occipital protuberance to the base of the tail and heart girth was determined by measuring the chest with a tape. Height was measured as a distance from the surface of the ground to the withers using a meter rule and chest depth was measured as the distance from the sternum at the region just in front of the forelimb to the withers, while chest width was recorded as distance between shoulders using meter ruler. Abdominal circumference was measured from the bottom of flank on one side to the bottom of flank on the other side of pig.

Back fat thickness

The back fat thickness of WT and leptin pigs was scanned ultrasonically. Briefly, pigs were immobilized and restrained by the head in a squeeze chute and the image site at the P2 position (at the level of head of last left rib) was determined by physical palpation. The pigs were held manually, avoiding any abnormal situation that could stress the animal, and were only scanned in a relaxed posture, permitting accurate measurements. A mix of Eco Gel and isopropyl alcohol was used as a sound-conducting material to allow a better acoustic contact surface between the probe and the skin. An Aplio 500 V real-time ultrasound machine (BSM34-0927, Toshiba, Japan) equipped with a 5.0–14.0 MHz linear array transducer (PLT-1005BT, Toshiba, Japan) was used for scanning the image site.

Reproductive performance

Reproductive parameters comprising of onset of puberty, pattern of estrous expression, and breeding performance were evaluated in leptin and WT pigs. The estrus observation in all pigs was started at d120 of age and gilts were visually observed two times (9:00 AM and 05:00 PM) a day for the expression of estrus signs in the presence of boar and stages of estrous i.e. proestrus, estrus, met-estrus, and di-estrus were ascertained according to previously described estrus scoring system [19] with slight modifications. The age at which gilts displayed the first estrus was recorded as the age of puberty and the pattern of estrous expression was recorded until 540 d. The inter-estrous interval was defined as the duration (days) from the start of one estrous to the start of the next estrous. Next, we evaluated the breeding performance of leptin and WT pigs by mating with a fertile male. Briefly, leptin females (n = 4) were used for mating, and gilts were mated 2 times in one estrous cycle using the same male. Usually, the first mating was performed 24 h after standing reflex followed by the second mating 18 h later. Pregnancy was monitored by non-return to estrous and confirmed by ultrasound scanner for around 33 d after service (HS-101 V, Honda Electronics Co., Ltd., Yamazuka, Japan).

Serum hormonal and lipid profile

Blood samples were collected at five different stages of the estrous cycle: proestrus, estrus, met-estrus, d6 post-met-estrus, and d12 post-met-estrus for hormonal analysis. Briefly, the pigs were fasted overnight and 5 ~ 10 mL blood was collected from the jugular vein using a vacutainer (Cat no. 367986, BD Vacutainer SST™, USA) and transported to the laboratory within 10 min. The serum was separated by centrifugation at 4˚C at 3000 rpm for 30 min and stored at -80˚C until assayed for hormonal analysis and cholesterol esters (CE) profile. FSH, LH, thyroid stimulating hormones (TSH), porcine prolactin (PRL), anti-müllarian hormone (AMH), and testosterone were measured by enzyme-linked immunosorbent assay (Shanghai Enzyme-linked Biotechnology Co., Ltd., China). The serum progesterone (P4) and estradiol (E2) concentrations were measured using a radioimmunoassay (Iodine [125I] Progesterone Radioimmunoassay Kit (B08B), and Iodine [125I] Estradiol Radioimmunoassay Kit (B05B), Beijing North Institute of Biotechnology, Co., Ltd. (Beijing, China). The serum leptin concentrations were measured using Porcine Leptin Elisa Kit (SEKP-0278, Solarbio life sciences, China).

Serum samples lipid extraction, separation of cholesterol esters and phospholipids via thin-layer chromatography, and determination of fatty acids via gas chromatography were conducted as described previously [20]. The fatty acid composition was analyzed with An Agilent 7890 A gas chromatographer equipped with 15 m × 0.25 mm × 0.25 μm DB-WAX columns (Agilent). While, C15:0 was used as a standard for quantitation. All the data were first normalized with the phospholipids profile and then presented as mean ± SD.

Luteal cell culture

Fresh ovaries of leptin and WT pigs were removed by laparotomy and immersed immediately in sterile phosphate buffer saline (PBS) with 5% penicillin-streptomycin solution and transported to laboratory. The luteal tissues were obtained from the ovaries and cut into small pieces using sterilized ophthalmic scissors, minced well and dispersed by pipetting in M199 (Sigma, M4530) medium containing collagenase (10 mg/mL) and the tissue was digested on a horizontal shaker in an incubator at 37 °C for 30 min. The supernatants containing luteal cells were then decanted through nylon mesh to remove debris. The cells were washed twice by centrifugation and resuspended in 10 ml Dulbecco’s modified Eagle’s medium with 10% fetal bovine serum and incubated at 37 °C, 5% CO2 in a 96-well plate for 24 h and then cells were collected and used for subsequent experiment.

Quantitative polymerase chain reaction (qPCR)

Total RNA from ovarian tissues and luteal cells were isolated using the Trizol reagent (Transgen Up, China) according to the manufacture’s instruction. cDNA was synthesized from the total RNA using a PrimeScript RT reagent Kit (TaKaRa, Japan) and was used as a template to perform qPCR in SYRB green-based qPCR instrument (CFX-96, Bio-Rad, USA). The reaction was performed in a 20 µL reaction mixtures comprising 10 µL of 2× SYBR (TaKaRa, Japan), 1 µL of cDNA, 1 µL of forward primer, 1 µL of reverse primer, and 7 µL of ddH2O. The reaction program as following: 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s and 62 °C for 45 s. The relative expression levels of target genes were quantified by  2-ΔΔCT. The primers are listed in additional file Table S3.

Ovarian histology and immunohistochemistry

For histological examination, ovaries were fixed in 4% paraformaldehyde for 48 ~ 72 h, processed by an automatic tissue processor (Yd-12p, Jinhua Yidi, medical appliance Co., Ltd, Jinhua, China) and embedded in a paraffin block (Yd-6D, Jinhua Yidi, medical appliance Co., Ltd, Jinhua, China). The paraffin blocks were cut into 5-um-thick sections using a Microm HM 325 microtome (Thermo Scientific, Waltham, MA, USA) and allowed to dry on glass slides overnight at 37 °C. Thereafter, the tissue sections were deparaffinized in xylene and rehydrated through graded ethanol dilutions. Sections were stained with hematoxylin–eosin (H&E) (G1120; Solarbio, China) according to manufacturer’s instruction.

For immunohistochemistry (IHC), after dewaxing and hydration, sections were incubated in 3% H2O2 solution for 30 min, and washed with PBS for thrice (each time 3 min). After that, the sections were blocked in PBS containing 5% BSA for 15 min at room temperature. Finally, the tissue sections were incubated with caspase-3 antibody (additional file Table S4) at 4 °C overnight. After washing with PBS for thrice, sections were incubated with specific secondary antibodies (R&D, USA) for 20 min. After washing thrice again (each time 3 min), sections were stained with fresh DAB (KIT-9901, Elivision TM plus Polyer HRP IHC Kit, Fuzhou, China) solution in dark for 5 min and then washed with PBS thrice (each time 3 min). Hematoxylin counter-staining, and neutral gum sealing of the slides was performed. Finally, the slides were imaged using OLYMPUS BX53 fluorescence microscope and analyzed using software of accessories.

Protein extraction and immunoblotting

The ovarian tissues from WT and leptin pigs were used to evaluate the expression of different protein levels (listed in additional file Table S4) using western blotting. In brief, the ovarian tissues were lysed in RIPA lysis buffer (Bestbio, China) with protease inhibitors at 4 °C. After lysis, supernatants were obtained by centrifugation at 13,000 × g for 15 min at 4 °C. Equal amounts of protein (100 µg) were run on SDS-PAGE gel, along with molecular weight marker. After electrophoresis, the proteins were transferred to PVDF membranes and reacted with primary antibodies against various antibodies and β-actin at 4 °C overnight. After incubation, the membranes were washed and incubated with specific secondary antibodies (R&D, USA). The membranes were then incubated with ECL (Easysee Western Blot Kit, China) and visualized with an Imaging System (Bio-Rad, Universal Hood II, USA).

RNA−seq analysis and quality assessment

Total RNA was extracted using TRIzol Reagent (Invitrogen, CA, United States) and purified using an RNeasy Mini Kit (Qiagen, CA, United States). The quality of RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent, Palo Alto, CA, United States). The rRNA-depleted RNA samples were further processed in accordance with the Illumina protocol (New England Biolabs, Massachusetts, United States). After cDNA synthesis, the samples were sequenced with an Illumina Novaseq6000 by Gene Denovo Biotechnology Co. (Guangzhou, China). The raw data were recorded. The overall quality of the RNA−seq data was evaluated by fastp (version 0.18.0). Clean reads were aligned to the Sus scrofa reference version 11.1 using HISAT2. 2.4 [21] with the default parameters.

Screening and clustering analysis of differentially expressed genes

Data preprocessing and follow-up analysis were performed DESeq2 software. The lists of DEGs between leptin and WT pigs were generated using the edgeR package (version 3.32.0) (Robinson et al., 2010). To normalize the raw data, the genes with the parameter of false discovery rate (FDR) below 0.05 and absolute fold change ≥ 2 were considered differentially expressed genes (DEGs) between the leptin and WT groups. Hierarchical clustering analysis was performed based on the expression levels of all transcripts and significantly differentially expressed transcripts using the pheatmap R package based on Euclidean distance.

Functional enrichment analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis in each module and network was conducted using the Database for Annotation, Visualization and Integrated Discover (DAVID). DEGs and enriched pathways were mapped using KEGG pathway annotation with KOBAS3.02. The top 5 KEGG pathways were selected and ranked by the enrichment factor. Subsequently, the Venn diagram tool was used to identify the common genes that were the focus of this work. Finally, gene set enrichment analysis (GSEA) using software GSEA and MSigDB [22] was performed to identify a set of genes in specific GO terms.

Single nucleus transcriptome sequencing (snRNA-Seq) libraries construction

Cellular suspensions were loaded on a 10X Genomics GemCode Single-cell instrument that generates single-cell Gel Bead-In-EMlusion (GEMs). Libraries were generated and sequenced from the cDNAs with Chromium Next GEM Single Cell 3’ Reagent Kits v2. following the manufacturer’s instructions.

Single-cell RNA-seq data processing

FASTQ files of 2 samples were processed with the use of Cell Ranger (v.3.1.0) count pipeline coupled with Sus scrofa reference version 11.1 to generate feature-barcode matrices. Seurat object list was then generated by Seurat package (v.4.0.4) [23], with R software (v.4.1.0) following these criteria: (1) min.cells = 5; (2) 200 < nFeature_RNA < 3600; (3) percent.mt < 0.1. In other words, genes expressed in at least 5 cells and gene number detected in cells ranging from 800 to 10,000 were kept for further analysis, and low-quality cells were filtered if R 10% UMIs derived from the mitochondrial genome.

Dimensionality reduction and clustering

To remove batch effects across samples, canonical correlation analysis method was used for data integration [24]. In detail, we initially normalized the filtered gene expression data with Normalize Data function with parameter (normalization.method = “LogNormalize”, scale.factor = 10000). Then top 2,000 variable genes were identified using the ‘vst’ method with FindVariableFeatures function in each sample, respectively. At last, we used FindIntegrationAnchors with parameters (i.e. anchor.features = 2000, k.filter = 200, dims = 1:30) and IntegrateData functions with the top 30 dimensions among 2 samples. After filtering low quality cells with 200 < nFeature_RNA < 3600, percent.HB < 5 and nCount_RNA < 10,000, we ran ScaleData function. To perform dimensionality reduction, the RunPCA function was conducted on linear-transformation scaled data with 2000 variable features and we performed UMAP with the top 30 dimensions. Finally, we clustered cells by using the FindNeighbors and FindClusters (resolution = 0.5) functions, which got 16 clusters.

Cell-type annotation and cluster marker identification

Cell clusters were obtained, after dimensional reduction and projection of all cells into two-dimensional space by UMAP. The Seurat FindAllMarkers function was used to identify markers for each cell cluster with the default settings. Canonical markers of specific cell types were used for cluster annotation. Ovarian cell-type markers were selected from PanglaoDB database(https://panglaodb.se/), CellMarker database(http://biocc.hrbmu.edu.cn/CellMarker/) and literature review.

DEG identification and functional pathways enrichment

The FindMarkers function in Seurat was used to identify differential expressed genes (DEGs) between two groups of cells with default parameters (logfc.threshold = 0.25, test.use = ‘‘wilcox,’’ min.pct = 0.1), and the enrich GO function in the clusterProfiler package (v.3.15.2) [25], was used to perform functional analysis with differential gene sets annotated with GO database (http://geneontology.org/). Enrichment pathways were obtained with parameters (pvalueCutoff = 0.05, pAdjustMethod = “BH”, OrgDb = “org.Ss.eg.db”).

Cell trajectory construction by Monocle

The subpopulation of germ line or granulosa cell (GC) line was imported into Monocle (version 2.20.0) to dissect the cell differentiation fate, also termed “pseudotime analysis.” With the gene count matrix as input, the new dataset for Monocle object was created, and functions of “reduceDimension” and “orderCells” were carried out to generate the cell trajectory based on pseudotime. Particularly, the ordering genes were DEGs between clusters in each cell type calculated by “differentialGeneTest” function in Monocle. In addition, the root state (that is, a prebranch in the heatmap) was set and adjusted following consideration of the biological meanings of different cell branches.

The regulon activity of transcription factors with SCENIC

The SCENIC algorithm was utilized to analyze the activity of transcription factors (TFs) and identify regulons (TFs and their target genes) in individual cells. The gene expression matrix, with genes in rows and cells in columns, was input into SCENIC (version 0.9.1) [26]. Co-expressed genes for each TF were constructed using GENIE3 software, followed by Spearman’s correlation analysis between the TFs and their potential targets. The “runSCENIC” procedure was then employed to generate the Gene Regulatory Networks (GRNs), also known as regulons. Regulon activity was analyzed using the AUCell (Area Under the Curve) software, applying a default threshold to categorize specific regulons as “0” (off) or “1” (on). t-SNE parameters were set to visualize the data with 50 principal components and a perplexity of 50, determined through a consistency test across multiple perplexity values and number of PCs. Cell states were mapped using specific regulons, and the average binary regulon activity was calculated.

Cell–cell communication analysis

Cell–cell communication was analyzed using iTALK package. The input data were Seurat object, and the ligand–receptor pairs were detected from top 50% highly expressed genes. The communication types mainly included growth factors, cytokines, and checkpoint. The network plot was visualized using LRPlot function.

Statistical analysis

Analyses were performed in the SAS statistical software (SAS ver. 9.4, SAS/STAT, SAS Institute Inc., Cary, NC). Data were expressed as least square mean ± standard error of mean unless otherwise stated. Continuous variables and measured over time, such as body weight and size including length, height, heart girth, chest depth, chest width and abdominal circumference, and hormones were analyzed by using the MIXED procedure of SAS (SAS/STAT) with estimations carried out by the method of restricted estimated maximum likelihood. All the mixed models included the repeated statement as day or stage of the estrous cycle and random effect of the animal identification nested within treatment for the proper error term. Continuous variable but measured on single point, such as serum leptin levels, pubertal onset, inter-estrus interval, back fat thickness, serum circulating cholesterol esters profile, relative expression of mRNA and protein were analyzed by PROC TTEST in SAS. Statistical significance was defined as *P < 0.05, **P < 0.01.