Bone Morphogenetic Protein 15 Knockdown Inhibits Porcine Ovarian Follicular Development and Ovulation

Bone morphogenetic protein 15 (BMP15) is strongly associated with animal reproduction and woman reproductive disease. As a multifunctional oocyte-specific secret factor, BMP15 controls female fertility and follicular development in both species-specific and dosage-sensitive manners. Previous studies found that BMP15 played a critical role in follicular development and ovulation rate in mono-ovulatory mammalian species, especially in sheep and human, but study on knockout mouse model implied that BMP15 possibly has minimal impact on female fertility of poly-ovulatory species. However, this needs to be validated in other poly-ovulatory species. To investigate the regulatory role of BMP15 on porcine female fertility, we generated a BMP15-knockdown pig model through somatic nuclear transfer technology. The BMP15-knockdown gilts showed markedly reduced fertility accompanied by phenotype of dysplastic ovaries containing significantly declined number of follicles, increased number of abnormal follicles, and abnormally enlarged antral follicles resulting in disordered ovulation, which is remarkably different from the unchanged fertility observed in BMP15 knockout mice. Molecular and transcriptome analysis revealed that the knockdown of BMP15 significantly affected both granulosa cells (GCs) and oocytes development, including suppression of cell proliferation, differentiation, and follicle stimulating hormone receptor (Fshr) expression, leading to premature luteinization and reduced estradiol (E2) production in GCs, and simultaneously decreased quality and meiotic maturation of oocyte. Our results provide in vivo evidence of the essential role of BMP15 in porcine ovarian and follicular development, and new insight into the complicated regulatory function of BMP15 in female fertility of poly-ovulatory species.


INTRODUCTION
In the past three decades, increasing studies have revealed the important role of the oocytespecific secreted factor, bone morphogenetic protein 15 (BMP15), in mammalian ovarian and follicular development through its multiple functions including promoting granulosa cells (GCs) proliferation and steroidogenesis (Otsuka et al., 2000(Otsuka et al., , 2001Moore et al., 2003;Moore and Shimasaki, 2005), preventing cell apoptosis and premature luteinization (Hussein et al., 2005;McNatty et al., 2005;Juengel et al., 2011;Chang et al., 2013;Zhai et al., 2013), regulating glycometabolism and lipid metabolism (Sugiura et al., 2007;Su et al., 2008), and controlling oocyte competence and ovulation (Fabre et al., 2006;Hussein et al., 2006). As a key signaling molecule mediating the dialogue between oocyte and its surrounding somatic cells (Gilchrist et al., 2008), BMP15 is expressed initially in the early follicle stage, and the expression gradually increases in subsequent follicle stages till the period of ovulation and/or luteinization (Paradis et al., 2009;Sun et al., 2010). This expression pattern differs with species, for example, the initial expression of BMP15 protein can be found in primary follicle (PF) stage of sheep, human, and pig, but not until the pre-ovulatory stage in mice (Paulini and Melo, 2011). BMP15 protein is secreted as BMP15/BMP15 homodimers and BMP15/growth differentiation factor 9 (GDF9) heterodimers, and it binds to the membrane bound type II serine/threonine kinase BMP receptor (BMPR2) and type I activin receptor-like kinase ALK6, resulting in the phosphorylation and activation of the SMAD pathways (Liao et al., 2003;Pulkki et al., 2012;Mottershead et al., 2013). In particular, BMP15 homodimers bind to ALK6 receptor to activate the Smad1/5/8 signaling pathway in some species, for example, in humans and sheep but not in rodents. However, BMP15/GDF9 heterodimers bind to BMPR2 receptor to activate Smad2/3 signaling pathway in all reported species (sheep, human, mouse, pig, etc.) (Reader et al., 2011(Reader et al., , 2016Peng et al., 2013;Lin et al., 2014). In most cases, BMP15/GDF9 heterodimers are more potent than BMP15 homodimers in the regulation of GCs, oocytes, and zygote development .
Bone morphogenetic protein 15 mutation or deficiency has been associated with altered female fertility in different species. As previously reported, natural mutations in BMP15 of sheep can lead to increased ovulation rate and litter size in heterozygotes but infertility in homozygotes due to bilateral ovarian hypoplasia (Braw-Tal et al., 1993;Smith et al., 1997;Galloway et al., 2000;Fabre et al., 2006). Altered fertility has also been reported in sheep immunized with BMP15 mature protein or different region of the peptides (Juengel et al., 2002(Juengel et al., , 2004(Juengel et al., , 2013McNatty et al., 2007). In humans, BMP15 mutations have been associated with primary ovarian insufficiency (POI) and infertility phenotype in women (Chand et al., 2006;Abir et al., 2014;Al-ajoury et al., 2015). However, in the polyovulatory mice, there was no significant difference between BMP15 ± females and wild-type in terms of the ovulation rate, and only a mild reduction of fertility was observed in BMP15 null female mice (Yan et al., 2001). Several studies have attempted to determine if there are species-specific differences in the BMP15 system that may play causal roles in the differences in fertility observed in poly-ovulatory mice and mono-ovulatory sheep and humans. One study has attributed the species-specific differences to the temporal variations in the production of the mature form of BMP15. The study found that mouse BMP15 mature protein was barely detectable until the pre-ovulatory stage, when it was markedly increased (Yoshino et al., 2006). Another study reported that defects in the production of mouse BMP15 mature protein correlate with species-specific differences (Hashimoto et al., 2005). Moreover, a phylogenetic analysis found a better conservation in areas involved in dimer formation and stability of BMP15 within mono-ovulatory species, but high variations in these areas within poly-ovulatory species, implying a correlation with altered equilibrium between homodimers and heterodimers, and modified biological activity that allows polyovulation to occur (Monestier et al., 2014). Hence, it seems that the role of BMP15 in the regulation of follicular development and ovulation rate was more critical in mono-ovulatory mammalian species than in poly-ovulatory animals. However, the role of BMP15 in ovarian and follicular development in polyovulatory mammalian species has remained unclear, as this has not yet been investigated in in vivo studies of non-rodent poly-ovulatory mammals.
In this study, we aim to investigate the function of BMP15 in female fertility and follicular development of non-rodent poly-ovulatory mammal by using a BMP15 knockdown transgenic (TG) pig model. The TG gilts had decreased female fertility with disordered estrous cycle, significant reduced ovarian size and follicle number, higher ratio of abnormal follicles, and none corpus lutein formed before 365 days old. We found that knocking down BMP15 can impair porcine follicle growth and cause dysovulation mainly by influencing oocyte quality and oocyte meiotic maturation, suppressing GCs proliferation and GCs functions, including inhibiting the expression of Fshr and E2 production, resulting in premature luteinization. These effects on follicular cell functions could finally lead to the absence of dominant follicle selection but appearance of abnormally enlarged antral follicles (AFs) with ovulation dysfunction in TG gilts. Our findings were evidently different from the unchanged fertility observed with BMP15 ± mice, strongly suggesting the important role of BMP15 in non-rodent poly-ovulatory mammals, thus providing the basis for further investigation on the different regulatory role of BMP15 in mono-and poly-ovulatory mammals. Table S1) targeting porcine Bmp15 mRNA was designed and selected by Invitrogen's web-based siRNA design software 1 . Human U6 promoter followed by each shRNA sequence was individually synthesized (Sangon Biotech, China) and cloned downstream of the EGFP expression cassette on pEGFP-N1 vector (Takara Bio, United States) to generate each pEGFP-Bmp15-shRNA expression vector ( Figure 1A). Meanwhile, a scramble shRNA expression vector was generated as negative control. To evaluate the RNA interference efficiency of shRNA, porcine Bmp15 CDS was synthesized (Sangon Biotech, China), and cloned into psiCheck II vector (Promega, FIGURE 1 | Generation of the BMP15 knockdown pig model. (A) Diagram of shRNA expression vector. Synthesized hU6-BMP15 shRNA fragment was inserted downstream of EGFP expression cassette on pEGFP-N1 vector. (B) RNA interference efficiency of five BMP15 shRNAs was examined by a dual-luciferase reporter system after 48 h transfection of h293T cells. NC, random shRNA plasmid. (C) PCR analysis of the muscle tissue proved that the integrated Bmp15 shRNA fragment had been transmitted to F1 gilts. +, TG gilt; -, sibling WT gilt; M, DNA Maker. (D) Southern blot analysis showed slightly less than 10 copies of constructed plasmids integrated in both F0 and F1 TG pigs, which was consistent with the result of about seven copies of qPCR analysis (data not shown). DNA with pEGFP-Bmp15 shRNA plasmid copies of 10, 20, and 40 were used as the positive control. (E) qPCR analysis of BMP15 mRNA level in 365 days old transgenic ovaries with two different phenotypes (TGF and TGS). TGF, transgenic ovary with many visible antral follicles on ovarian surface. TGS, transgenic ovary with streak phenotype. * p < 0.05. (F) Western blot analysis of BMP15 protein level in postnatal 30-day old TG ovaries. Three prominent, distinct bands were observed corresponding to apparent molecular weights of 34 kDa, 27 kDa, and 15 kDa. (G) Quantitative analysis of BMP15 protein levels based on the band intensity using Image J software. (H) F1 TG gilt showed a visible intense GFP fluorescence on the toes and muscle while subjected to sunlight.

Five shRNAs (Supplementary
United States) to generate psiCheckII-Bmp15 plasmid. Each pEGFP-Bmp15-shRNA plasmid was respectively co-transfected with psiCheckII-Bmp15 plasmid into HEK293 cells. After 48 h of culture, transfected cells were collected and subjected to RNA interference efficiency detection by using a dual-luciferase reporter system (Promega, United States). The shRNA with the most efficient RNA interference was selected for the generation of BMP15 knockdown pig model.

Generation of Bmp15 Knockdown Pig Model
Procedures for the generation of the Bmp15 knockdown gilts are illustrated in Supplementary Figure S1. Briefly, the selected pEGFP-Bmp15 shRNA plasmid was transfected into PEFs derived from a male Yorkshire pig. After G418 selection and fluorescence examination, EGFP-positive PEFs were used as donor cells for somatic cell nuclear transfer (SCNT). For SCNT, oocytes were collected from abattoir ovaries with a 20G needle connected to a syringe, and then cultured in HEPES-buffered tissue culture medium 199 and later maturation medium, until in vitro maturation. SCNT by handmade cloning and embryo transplantation were carried out by BGI Ark Biotechnology Company, China. After 114 days of pregnancy, we obtained two surviving F0 generation Bmp15-knockdown TG males. Next, we mated one TG boar with wild-type sows through artificial insemination (AI), and obtained F1 generation TG gilt for this study. Sibling gilts without pEGFP-Bmp15-shRNA integration were used as controls (WT) in this study. The protocol for the animal study was approved by the Institutional Animal Care and Use Committee (IACUC), Sun Yat-sen University (Approval Number: IACUC-DD-16-0901).

Tissue Collection
A total of 54 animals including 25 WT and 29 TG gilts were sacrificed at ages of 30-500 days. Among these TG gilts, six gilts contained eight TGS ovaries at age of 110, 160, 200, and 365 days. Tissues of ovary, pituitary, muscle, liver, kidney, heart, and uterus were collected. Tissues used for RNA extraction were directly soaked in TRIzol reagent (Promega) and frozen rapidly in liquid nitrogen. Muscle tissues for DNA extraction together with the 30-day ovarian tissues for protein detection were directly frozen in liquid nitrogen before being transported to the laboratory. All the other ovaries were washed in sterilized saline water and photographed. Ovaries at the ages of 60 and 90 days (n = 6) were used for mRNA detection, and 365-day old ovaries (n = 6) of WT and TGF gilts were used for mRNA detection in ovarian tissues and isolated follicle, respectively. Primers for qPCR analysis are shown in Supplementary Table S3. Ten ovaries (5 WT and 5 TGF) at the age of 60-170 days were frozen in OCT and stored at −80 • C before laser capture microdissection (LCM). Six ovaries (3 WT and 3 TGF) from different 365-day old gilts were used for dissection and follicular fluid collection. Furthermore, they were used for COCs collection and subsequent single-cell sequencing. The other ovaries (24 WT, 24 TGF, and 8 TGS) at age of 30-400 days were used for hematoxylin and eosin (HE) observation and immunohistochemistry (IHC) analysis; they were fixed in 10% (w/v) paraformaldehyde/0.02 M PBS (pH 7.2) on ice before transportation to the laboratory.

Identification and Characterization of Transgenic Gilts
Transgenic pigs were first screened by GFP fluorescence on the toes under sunlight, and the integration of pEGFP-Bmp15 shRNA plasmid in genome using genomic DNA extracted from muscle tissues was confirmed by PCR analysis (Supplementary Table S2). The copy number of integrated plasmid was determined by qPCR and Southern blot analysis. Bmp15 mRNA expression level in TG pigs was detected by qPCR in 365-day old ovaries, and BMP15 protein level was detected by western blot analysis in 30-day old ovarian tissues.

F1 Gilts Estrous Checking and Hormone Assays
About 50 F1 TG gilts at age of 170-400 days were checked daily for signs of estrous in the presence of an intact mature boar. TG and WT gilts (n = 2, each) at about 365 days old were chosen for daily vaginal smears analysis, and daily jugular venous blood collection at 9:00-11:00 AM for 24 days continuously. Vaginal cell smears analysis and estrous identification were performed as described in a previous report (Mayor et al., 2007). Daily blood samples were centrifuged at 1500 g for 15 min, and the serum samples were collected and stored at −80 • C. These serum samples were thawed on ice prior to use in the quantification of the concentration of estradiol (E2) and progesterone (P4) by chemiluminescence immunoassay (CLIA) (Siemens, Germany).

Histological Examination
Ovaries derived from gilts at age of 30-400 days were fixed in 10% (w/v) paraformaldehyde with 0.02 MPBS (pH 7.2) at 4 • C for about 2 h. Next, they were cut into vertical slices of about 0.5 cm thickness and fixed in fresh 10% (w/v) paraformaldehyde for a cumulative period of 24 h. These slices were mounted in paraffin, and serially cut into 5 µm-thick sections by Rotary Microtome (MICROM, Germany), and stained with HE. Ovarian HE sections were observed and photographed under a fluorescent microscope (Zeiss, Germany).

Ovary Dissection and Follicular Fluid Collection
Each of the three 365-day TGF and WT ovaries derived from different individuals were flushed with sterilized saline water and placed in the incubator at 38 • C during transportation to the laboratory. Later, visible AFs in these ovaries were dissected by scalpel blade and tweezers, and classified into three groups (1-3 mm, 3-5 mm, and >5 mm) according to their diameter, which was measured by a Vernier caliper. Total follicle number of each group was counted. Next, follicular fluid from AFs with diameter of 3-5 mm and diameter >5 mm were collected by a disposable 10 mL syringe. The concentrations of follicle stimulating hormone (FSH), LH, E2, and P4 in follicular fluid were quantified by the CLIA method (Siemens, Germany).

Laser Capture Microdissection (LCM)
A total of 10 ovaries from WT and TGF gilts of age 60-170 days were embedded in OCT and placed on a cryostat (MICROM, HM560, Germany). All ovaries were cut into 7 µmthick sections and mounted on RNAse free membrane slides (MMI, 50102). These membrane slides were then fixed in icecold 95% ethanol for 1 min, and later washed in 75% ethanol for 30 s. Subsequently, the sections were stained following a previously reported method (Golubeva et al., 2013). Briefly, the staining mixture was prepared with 1% cresyl violet in absolute ethyl alcohol, Eosin Y, RNAse free water, and 100% ethanol at the ratio of 3:1:4:4. Membrane slides were stained in this fresh staining mixture for 30 s and dehydrated thrice in 100% ethanol for 1 min, followed by 30 s incubation in xylene. Slides were finally dried for 5 min by a hair dryer blowing cold wind, and stored at −80 • C until use.
The follicles were distinguished from each other as follows: PF was defined by a clear monolayer of cuboidal GCs; secondary follicle (SF) was defined by more than two layers of GCs but without any antrum; small antrum follicle (SAF) was defined by obvious small antrum but not completely separated granulose and cumulus cells; AF was characterized by a big single central antrum with completely separated granulosa and cumulus cell layers. The entire PF, SF, and SAF, but only parietal granulosa and theca cells of AF (APC) were isolated by LCM. Each type of the follicles of each ovary was dissected under 20 × magnification microscopic visualization using MMI Cell Cut Plus system (MMI, Swiss). Later, the dissections were treated with 100 µL of TRK Lysis buffer of the MicroElute total RNA kit (Omega) and 2 µL 2-mercaptoethanol. Both TGF and WT lysates were, respectively mixed according to their follicle stages after lysing for 10 min at 16 • C, followed by storage on dry ice until RNA extraction. A total of eight LCM-derived RNA samples, PF WT , SF WT , SAF WT , APC WT , PF TGF , SF TGF , SAF TGF , and APC TGF , were used for transcriptomic analysis. RNA-seq was performed on an Illumina HiSeq 2000 using Illumina TruSeq SBS kit v2 (209 cycles including index) to obtain paired-end reads (2 × 100 bp).

Single-Cell RNA Sequencing on Cumulus-Oocyte Complex (COCs)
Cumulus-Oocyte Complex were collected from large AFs (diameter, 5-7 mm) of 365-day old TGF and WT ovaries (n = 3) derived from different gilts by using a 20-gauge needle fixed to a 10 mL disposable syringe. They were pooled and placed on a stereomicroscope (Nikon). COCs with several layers of cumulus cells and uniform cytoplasm were selected. The selected COCs from both TGF and WT ovaries (n = 10 from each) were used for RNA microextraction by MicroElute total RNA kit (Omega). Total RNA was pre-amplified by SMARTer R Ultra TM Low RNA Kit (Clontech), and sequenced on Illumina HiSeq 2000 sequencing system.

Data Analysis
All the mRNA detections were performed three times, and quantitative data are represented as the mean ± SEM. Measurements of hormones concentration were not repeated but were conducted on at least three individuals. Statistically significant differences were analyzed by t-test or one way analysis of variance, followed by Duncan's post hoc test. A p-value < 0.05 indicated statistical significance.
For RNA-seq data analysis, raw RNA-Seq clean reads were obtained by removing reads containing low quality reads and/or adaptor sequences from raw reads and mapped to the pig genome (Sus scrofa 10.2), allowing up to two base mismatches. Differential expression analysis was performed using the Benjamini approach; genes with an adjusted p-value < 0.05 and log2 fold change >1 were designated as differentially expressed genes (DEGs). DEGs lists were submitted to the databases of Novogene Company (China) for further enrichment analysis. GO analysis was performed with Webgestalt software. In all the tests, p values were calculated using the Benjaminicorrected modified Fisher's exact test, and p < 0.05 was taken as a threshold of significance. Venn diagrams were drawn using the web tool 2 . As determined by gene co-expression analysis, a correlation coefficient of 0.98 was set as the threshold value. Closely correlated genes were imported in cytoscape software to generate the co-expression network.

Generation and Identification of BMP15 Knockdown Pig Model
In our study, we generated a BMP15 knockdown model through the shRNA technique. We designed and constructed five pEGFP-BMP15-shRNA plasmids, in which BMP15 shRNA sequence was under the control of human U6 promoter (Czauderna, 2003) and inserted downstream of the EGFP expression cassette ( Figure 1A). Each shRNA expression plasmid was co-transfected into HEK293T cell with a psiCheckII-BMP15 plasmid to examine their RNA interference efficiency in vitro using a dual-luciferase reporter system (Boudreau et al., 2009). We found that shRNA1 was the most effective with a RNA inference efficiency reaching 76% (Figure 1B), and thus, this shRNA was selected for transfection into embryonic fibroblast cells (PEFs) derived from a male Yorkshire pig. Transfected PEFs were then subjected to G418 selection to select the cells with stable expression of EGFP as donor cells for SCNT (Baguisi et al., 1999;Polejaeva et al., 2000;Liu et al., 2019;Supplementary Figure S1A).
Both F0 and F1 TG pigs showed visible intense GFP fluorescence on the toes and muscle when subjected to sunlight ( Figure 1H), directly suggesting that the pEGFP-BMP15 shRNA plasmid was successfully integrated into the genome of F0 TG boar, and can be transmitted to the next generation through the germline. This was confirmed by PCR analysis of Bmp15 shRNA fragment in the muscle tissue of F1 TG gilts ( Figure 1C). The copy number of integrated plasmid was estimated to be approximately seven in F1 TG pigs as analyzed by the combination of qPCR using a transferrin receptor gene to normalize the genomic DNA (data not shown), and Southern blot analysis ( Figure 1D). More importantly, evident decreased levels of BMP15 mRNA ( Figure 1E) in 365-day old TG ovaries and BMP15 protein in 30-day old TG ovaries ( Figures 1F,G) strongly demonstrated the successful generation of a BMP15 knockdown model, and implied an in vivo BMP15 knockdown efficiency of about 50% in TGF ovaries. Our qPCR data also revealed that BMP15 mRNA was highly expressed in the ovary tissues, lowly expressed in the pituitary (Supplementary Figure S1B), but was undetectable in other porcine tissues (e.g., liver, muscle, and kidney).

Knockdown of BMP15 Was Associated With Disordered Reproductive Cycle of TG Gilts
All F1 TG gilts presented normal appearance and growth condition, and 50 of them at ages of 170-400 days were checked daily for signs of estrous in the presence of an intact mature boar. Surprisingly, we did not observe any obvious estrous behavior or vulvar appearance changes (e.g., increased redness, swelling, or mucus production) (Soede et al., 1994) in sexually mature TG gilts (Supplementary Figure S2A). About twenty gilts at ages of 240-400 days were bred by AI after treatment with PG600, but all failed to become pregnant. To determine if the disordered estrous cycle in TG gilts was related to reproductive hormonal changes, we measured the concentration of plasma estrogen (E2), progesterone (P4), and FSH in 365-day old gilts throughout the estrous cycle. The results showed that a typical peripheral E2 concentration peak before the onset of estrous can be observed in WT gilts, which is consistent with previous studies (Noguchi et al., 2010;Soede et al., 2011; Figure 2A). In contrast, irregular E2 concentration peaks were observed in TG gilts during a 24-day continuous measurement (Figure 2B). Cytological analysis of vaginal smears (Mayor et al., 2007) further proved that disordered estrous cycle occurred in TG gilts, as irregular cytologic changes was observed throughout a 16-day continuous examination (Supplementary Figure S2B). Furthermore, except for TG3 gilt, the average P4 level of the other two gilts was significantly lower than that of the WT gilts, with TG1 gilt containing bilateral asymmetric ovaries [one streak ovary (TGS ovary), and the other one TGF ovary] presenting the lowest serum P4 level. The other two TG gilts contain bilateral TGF ovaries (Figure 2C), and higher serum FSH concentration was found in TG2 and TG3 gilts ( Figure 2D). In addition, we found over 2-fold up-regulated expression of Fsh mRNA level in the pituitary of both 150-day and 260-day old TG gilts ( Figure 2E), but no significant difference in the expression level of luteinizing hormone (Lh) (Figure 2F). These results indicate a disordered reproductive cycle and potential ovarian dysfunction in TG gilts.

Knockdown of BMP15 Led to Inhibition of Follicular Development and Ovulation in TG Ovaries
Because the estrous cycle is determined by ovarian and follicular development (Foxcroft and Hunter, 1985;Noguchi et al., 2010), the disordered estrous cycle of TG gilts is potentially caused by impaired ovarian follicular development. In this regard, ovaries from gilts of different ages were collected and processed for morphological examination. Surprisingly, we found remarkably decreased size in TG ovarian and number of AFs on the surface of TG ovaries of 140-365-day old gilts. In addition, apparent size difference was observed between bilateral TG ovaries ( Figure 3A). Corpus lutein was not observed in TG ovaries from 140-365-day old gilts, but can be found in 400and 500-day TGF ovaries ( Figure 4E). Besides, the weight of TG ovaries before sexual maturity was markedly lower than that of WT ovaries (Supplementary Table S4). Among the TG ovaries, we discovered eight streak ovaries in six gilts, denoted as TGS ovaries with an incidence of about 14%, whereas no streak ovary was found in WT ovaries (Figures 3A,B). These TGS ovaries contained none or less than three visible AFs on the ovarian surface. In the cortex of TGS ovaries from 110-and 200-day old gilts, most follicles were arrested in the primary stage (Figures 3B,E). In the cortex of TGS ovaries from 365day old gilts, most follicles were arrested in the secondary stage, and the degradation of follicles became apparent ( Figure 3B). Unlike TGS ovaries, the rest of TG ovaries contained many visible large AFs on the surface, and we denoted them as TGF ovaries (Figures 3A,C). Different stage of follicles can be found in these TGF ovaries; however, the follicle number decreased drastically during follicular development (Figures 3C-E). Notably, during the early follicle stage, the significant decline in primordial and PF number led to a much thinner ovarian cortex in the TGF ovaries of pre-puberty gilts. In addition, structural abnormality of SFs was evident, particularly in the TGF ovaries of gilts undergoing puberty (Figures 3C,G,H and Supplementary Figure S3A). We observed abnormally enlarged ( Figure 3H) or degenerated oocytes (Figure 3Gii), multioocytic follicles (Figure 3Gi), highly irregular GC layers (Figures 3Gii,iv), degraded GCs (Figure 3Giv), and abnormally thickened theca layers (Figure 3Giii). Furthermore, ovarian sections from five TG gilts at the age of 160-400 days showed a markedly reduced proportion of normal SFs in the TGF ovaries ( Figure 3F and Supplementary Figure S3A).
Histological observation also revealed some striking features of TGF AFs. Most notably, AF number declined remarkably, but its antrum was enlarged substantially ( Figure 3C) and surrounded with loosely organized smaller GCs ( Figure 3I). We further isolated AFs from three 365-day TGF ovaries derived from different gilts for statistical analysis. The results showed that TGF ovaries contained less total number of AFs and less number of small AFs (diameter <5 mm); however, it contained FIGURE 2 | Transgenic gilts presented disordered estrous cycle and reproductive hormones. Plasma E2 and P4 concentration of 365-day old WT and TG gilts were measured at a 24-h interval for 24 days. (A) During the estrous cycle, representative WT gilts showed typical serum E2 concentration peak before ovulation, accompanied with marked decreased P4 concentration (n = 3), data are presented as mean ± standard deviation. (B) Irregular plasma E2 concentration peaks was observed in three representative TG gilts in a continuous 24-day measurement. (C) The average P4 concentration of two of the three TG gilts in a continuous 24-day measurement was significantly lower than that of WT gilts (p < 0.05). Each point stands for a P4 concentration value. (D) Two of the three 365-day TG gilts showed higher serum FSH concentration. Serum FSH concentration was measured at a 24-h interval for 24 days. (E) Expression level of Fsh mRNA in the pituitary of both 150-and 260-day old TG gilts were two-fold higher than that in WT gilts (p < 0.05). (F) The average level of Lh mRNA level in pituitary was not significantly different in TG and WT gilts. * Stands for statistical significance (P < 0.05).
substantially more large AFs (diameter >5 mm) (Figures 4A,B). We found that the diameter of the largest AF in TGF ovary can be up to 9 mm, whereas it is about 7 mm in WT ovaries ( Figure 4C). Normally, porcine AFs stop growth at a diameter of about 5 mm, and only selected follicles continue to grow through the accumulation of follicular fluid and ovulate at a diameter of about 7 mm (Soede et al., 2011). Thus, the increased numbers of abnormally enlarged AFs in the TGF ovaries may be related to dysovulation and the disordered serum reproductive hormones found in TG gilts. Subsequent measurement of the concentration of reproductive hormones in follicular fluid of TGF large AFs (diameter >5 mm) showed that the E2 concentration was remarkably lower (Figure 4D), but the concentrations of the other three hormones including P4 (Figure 4F), FSH (Supplementary Figure S3B), and LH (Supplementary Figure S3C) was not significantly different from those in WT large AFs. Reduced E2 production in the TGF large follicles may imply an absence of dominant follicle selection (Clement and Monniaux, 2013). Taken together, these results provide convincing evidence that the knockdown of BMP15 FIGURE 3 | BMP15 knockdown-induced changes in ovarian and follicular development. (A) Representative photograms of ovaries collected from gilts of different ages showing reduced ovarian size and less visible follicles on the surface of TG ovaries as compared to ovaries from WT sibling. Bilateral TG ovaries were significantly different in size at ages of 200 and 365 days. Two ovarian phenotypes were found in TG ovaries, TGF ovaries had many visible large antral follicles on the ovarian surface. TGS ovaries contained none or less than three visible antral follicles (red arrows). (B) Histological observation of TGS ovaries showed that 110-and 200-day old TGS ovary presented major primary-like follicles sparsely scattered on the cortex (green arrows), while 365-day old TGS ovarian section was predominantly occupied by degraded secondary follicles (black arrow). (C) On the 200-day TGF ovarian section, decreased number of follicles and enlarged antral follicles (black arrow) was observed. In addition, the degradation of GCs in abnormally organized GC layer structure of secondary follicles was observed (black arrows). (D) In 30-and 80-day old TGF ovaries, drastically decreased number of early stage follicles led to thinner ovarian cortex (blue line and blue arrows). (E) Comparison of three ovarian phenotypes at the age of 110 days showed less number of early stage follicles in TGF ovarian cortex, and the minimum number of follicles in TGS ovaries. (F) Markedly reduced proportion of normal secondary follicles in the TGF ovaries. Secondary follicles in three sections of each ovary were counted. (G) Representative images of abnormal TGF secondary follicles (black arrows), including multiovular follicle with highly irregularly organized theca cell layers (i); follicle with oocyte-free structure, and abnormally thickened zona pellucida surrounded by highly degraded GCs (ii); follicle with abnormally thickened theca layers (iii); follicle with enlarged oocyte surrounded by highly irregularly organized GC layers with holes formed by degradation of GCs (iv). (H) TGF follicle showed larger oocyte in the early secondary follicle stage (black arrow head). (I) Smaller GCs were loosely organized in TGF antral follicles (black arrow head).
could severely inhibit both follicular development and ovulation of poly-ovulatory pig.

Knockdown of BMP15 Caused Premature Luteinization and Impaired Oocyte Quality in TGF Follicles
We examined the expression and activation of factors associated with follicular development. The results confirmed that BMP15 protein is abundantly located in both WT oocytes and GCs of primary to pre-ovulatory follicles. Both normal and abnormal TGF follicles showed slightly decreased BMP15 protein accumulation in the less degraded GCs than in WT. However, markedly reduced BMP15 expression level was noted in deteriorated oocytes of TGF abnormal (TGFA) follicles. TGS ovaries exhibited the minimum BMP15 protein level in the PFs of 110-day TGS ovaries and highly degraded SFs of 365-day TGS ovaries (Figure 5A and Supplementary Table S5). Thus, we speculated that the in vivo BMP15 interference efficiency is different in TG individuals, which possibly explains the two TG ovarian phenotypes (TGF and TGS). Besides, TGS ovaries displayed a phenotype of high degradation and serious inhibition of follicular development and cellular activity in the arrested SFs, similar to the phenotypes of BMP15 homozygotes mutations in sheep (Braw-Tal et al., 1993) and in women with POI (Luisi et al., 2015). Hence, we focused on the effects of TGF follicles.
In TGF follicles, we found that the expression patterns of both GDF9 and FSHR, the cooperator and down regulator of BMP15, respectively, corresponded with that of BMP15 ( Figure 5A, Supplementary Figure S4, and Supplementary Table S5) in TGF follicles. However, there were no changes in the expression of BMP15 receptors (ALK6 and BMPR2) (Supplementary Figure S4). Unlike BMP15, the expression of LHR was higher in TGF follicles than in WT ( Figure 5A). The excess expression of LHR suggested premature luteinization in TGF follicles, which was also demonstrated by the significantly raised expression of 3β-hydroxysteroid dehydrogenase (Grasa et al., 2016) in TGF SFs. In consideration of the striking features of reduced follicle number and degraded GCs in the TGF follicles, we determined the expression levels of caspase 3 and Ki67 to assess cell apoptosis and proliferation activity. Surprisingly, there was no change in the expression of both caspase 3 and Ki67 in the degraded TGFA follicles (Figure 5B). Subsequent investigation of BMP15-mediated signaling pathways suggested an underlying mechanism. We observed notably weakened Smad1/5/8 activity in TGFA follicles when compared to TGF normal SFs (Figure 6A and Supplementary Figure S5), but slighter attenuated Smad2/3 phosphorylation was shown in the abnormal follicles ( Figure 6B and Supplementary Figure S5). It is likely that Smad1/5/8 mainly contributed to the inhibition of follicular development in the TGFA follicles, and Smad2/3 activation by the BMP15-independent pathway played a role in the growth of the less degraded follicles. Except for GCs, we also discovered impaired oocyte quality in TGFA follicles. As shown in Figure 5C, an undetectable level of autophagy-related protein LC3B (microtubule-associated protein 1 light chain 3) (Jiang et al., 2017) was observed in the oocytes of TGFA SFs, whereas the oocytes of WT and TGF normal follicles displayed a normal level of LC3B. This result demonstrates that autophagy activity of the oocytes of TGFA follicles, which is fundamental to the cellular processes of the oocytes, was largely weakened (Su et al., 2017).

Knockdown of BMP15 Resulted in Dynamic Transcriptomic Alteration During TGF Follicular Growth
To further investigate the regulatory role of BMP15 in porcine follicular development, RNA-seq was carried out on follicles or GCs captured by LCM method from frozen sections of both WT and TGF ovaries (Fend et al., 1999;Bonnet et al., 2015). LCMcaptured follicles were categorized into three stages of follicular development: PF, SF, and SAF stages. For large antrum follicle, only parietal granulosa and theca cells (APC) were captured by LCM for RNA-seq. The follicles or APCs were captured from the frozen sections of each of the five TGF and WT ovaries of gilts at age ranging from 60-170 days. The gene profiles of 34,640 genes generated by RNA-seq were used for the identification of DEGs, as well as GO and pathway enrichment analysis basing on intra-(between each two continuous follicle stages in either WT or TGF sample) and inter-(between each follicle stage of WT and TGF sample) effect comparisons ( Table 1). In intra-effect comparisons, the largest number of DEGs (3503 DEGs) was found in SAF WT /SF WT comparison, and the least number of DEGs (350 DEGs) was found in APC WT /SAF WT comparison. However, in contrast to WT, during TGF follicular development, the lowest number of DEGs was Statistical data showed reduced total number of antral follicles, and substantially increased number of follicles with a diameter >5 mm in 365-day TGF ovaries. Antral follicles were isolated from three ovaries of different gilts, and then classified into three groups according to their diameter (d 1-3 mm, d 3-5 mm, and d >5 mm). (C) Comparison of the three largest follicles isolated from WT and TGF ovaries. (D) E2 concentration in follicular fluid of TGF large antral follicles was significantly lower than that in WT pre-ovulatory follicles. * p < 0.05. (E) Many corpus lutea were observed on the surface of both 400 and 500-day old TGF ovaries, but the ovaries were apparently smaller than 400-day old WT ovaries. (F) P4 concentrations in the follicular fluids of TGF and WT were not significantly different. The hormones in follicular fluid were measured in follicles of three 365-day old ovaries of different gilts.
found in SAF TGF /SF TGF comparison, and the largest number of DEGs was found in SF TGF /PF TGF comparison. A striking difference in the dynamics of transcriptions between WT and TGF follicular development was also found in GO ( Figure 7A) and pathway (Supplementary Figure S6A) enrichment. In WT ovary, more DEGs were presented in the enriched GO during the dynamical transition from early primary to SF stage, and from secondary to SAF stage, but less DEGs were presented in enriched GO during SAF to large antrum follicle stage transition. However, in TGF ovary, more DEGs were presented in enriched GO during the dynamical transition from SAF to large antrum follicle stage, but less DEGs were presented during secondary to SAF stage transition. These results in GO enrichment were consistent with that found in pathway enrichment. Based on the intraeffect analysis, GCs differentiation during the late secondary stage and early antrum formation may have been delayed during FIGURE 5 | Abnormal TGF follicles showed premature luteinization and impaired oocyte quality. (A) Immunostaining of ovarian sections indicated that the expression of BMP15 decreased remarkably in TGS abnormal follicles, but was only slightly reduced in TGF follicles, compared with WT follicles. FSHR shared a similar expression pattern with BMP15, and the expression pattern of LHR was different from that of BMP15. It expressed higher in TGF GCs of both preantral and antral follicles. Scale bar = 100 µm. (B) Follicular cell apoptosis, proliferation, and premature luteinization were evaluated by immunostaining with Caspase3, Ki67, and 3βHSD, respectively. Notably, higher expression level of 3βHSD was discovered in abnormal TGF follicles. However, the expression of Caspase3 and Ki67 in abnormal TGF follicles was not significantly different from that in WT follicles. Scale bar = 100 µm. (C) Immunofluorescence images demonstrates intensive expression of autophagy-related protein, LC3B, in oocytes of normal follicles of TGF and WT ovaries, but it was barely expressed in oocytes of abnormal follicles of TG (TGF and TGS) ovaries. Purple arrow, oocytes in normal follicles; Orange arrow, oocytes in abnormal follicles. Scale bar = 100 µm.
TGF follicular dynamical development (Hennet and Combelles, 2012). BMP15 probably played a more important role during the dynamic development of secondary and subsequent follicle stages than in the early stages.
Considering that the expression and function of BMP15 during follicular development are stage-specific (Paradis et al., 2009;Sun et al., 2010), we conducted the RNA-seq data analysis based on inter-effect comparisons. Both DEGs number and their clustering results revealed distinct effect on the knock down of BMP15 in each follicle stage, where the largest number of DEGs was found in the SF stage ( Figure 7B). Based on these DEGs, we enriched a total of 15 up-regulated and 26 down-regulated pathways (Supplementary Figure S6B) in all the four follicle stages. Ten pathways were enriched during the three dynamical transitions of each of the two continuous follicle stages, in which DEGs were identified based on the combination of the intraand inter-effect comparison (Supplementary Figure S7B). These pathways are illustrated in Figure 7D according to their relevant function. Interestingly, the knock down of BMP15 seemed to lead to an increased cell growth in PF, due to the significant up-regulation of estrogen and MAPK pathways and the downregulation of ubiquitin protein degradation pathway. However, knocking down BMP15 likely inhibited GCs proliferation and growth from the SF stage onward, because of the significant down-regulation of pathways such as cell cycle, Hippo, P53, and PI3K, which also implies the potential involvement of BMP15 in these pathways. Though GCs proliferation and growth were inhibited by knocking down BMP15, the significantly upregulated MAPK pathways in primary and SF stages, and Notch pathway in secondary and SAF stages may compensate for the FIGURE 6 | Smad1/5/8 signaling pathway was inhibited in abnormal follicles of TGF ovaries. (A) Immunofluorescence images showed Smad1/5/8 pathway was evidently less activated in abnormal follicles of TGF ovaries, and severely inhibited in highly degraded 365-day old TGS follicles, compared with normal follicle of TGF and WT ovaries. (B) Immunofluorescence images demonstrated a mild decrease in Smad2/3 signaling in both normal and abnormal follicles of TG ovaries. Smad2/3 signaling was remarkably inhibited in highly degraded 365-day TGS follicles. TGFN, normal follicles in TGF ovary; TGFA, abnormal follicles in TGF ovary. Scale bar = 100 µm.
knock down by promoting cell proliferation and growth in preantral follicles to support the continuous development of certain number of follicles to maturity. Furthermore, knocking down BMP15 did not result in significant pathway alteration during the transition from secondary to SAF stage ( Figure 7D). In addition to the fact that the minimum DEGs was presented in SAF TGF /SF TGF comparison (Table 1), our findings also suggest that knocking down BMP15 may cause an inhibition of GCs differentiation and abnormal development in preantral follicles. In the large antrum stage, DEGs involved in ovarian steroidogenesis (Lhr,Cyp17,3βHsd,etc.) were significantly up-regulated in theca cells of TGF follicles (Supplementary Table S6), possibly associated with premature luteinization of TGF follicles (Figure 5B). Except for GCs, we also enriched FIGURE 7 | Altered follicle dynamic transcriptomes during TGF follicular development. (A) Different numbers of DEGs were presented in enriched 12 biological processes (GO) of WT and TGF dynamic transcriptomes during follicular development. In TGF ovary, the least DEGs was presented in enriched GO during the transition from secondary to small antrum follicle stage; in WT ovary, the least number of DEGs was presented in enriched GO during the transition from small antrum to antral follicle. The highest number of DEGs was presented in enriched GO during the transition from primary to secondary follicle stage in TGF ovary, whereas in WT ovary, this was presented during the transition from secondary to small antrum follicle stage. Gene ontology-based analysis was conducted with Webgestalt software. (B) Both DEGs number and their clustering pattern revealed a significant difference between the gene expression of WT and TGF follicle during each follicle stage. (C) A subset of DEGs was validated by qPCR. 60 d stands for 60-day ovarian tissue, which was mainly composed of early stage follicles (primordial, primary, and early secondary follicles). 90 d represents 90-day ovarian tissue, which was mainly composed of secondary follicles without antral follicles. AF <5 mm, antral follicle with a diameter <5 mm. AF 5-6 mm, antral follicle with a diameter of 5-6 mm. AF >7 mm, antral follicle with diameter >7 mm. (D) Summary of important pathways (based on inter-effect analysis) involved in follicular development. Up-regulated pathways are shown in red, and down-regulated pathways are shown in green. * Stands for statistical significance (P < 0.05).
significantly down-regulated DEGs involved in oocyte meiosis and maturation in TGF follicles beyond the PF stage, which was possibly related to the impaired oocyte quality ( Figure 5C).
Moreover, we found that in vivo knock down of BMP15 resulted in significant decrease in Bmp15 expression level from primary to SAF stage (Supplementary Table S6), which was confirmed by qPCR analysis (Figure 7C). Through a correlation analysis of the DEGs, we predicted 13 downstream regulated genes of Bmp15, in which six genes (Atrx, Amd1, Dtd2, etc.) were positively correlated, and seven genes (Fgf9, Igfbp7, Cmpk2, etc.) were negatively correlated (Supplementary Figure S7A and Supplementary Table S7). Furthermore, an unexpected significantly decreased expression of Fshr in TGF follicles was detected by transcriptomic analysis (Supplementary Table S6), and confirmed by qPCR (Figure 7C), which seemed to be inconsistent with the previous perspective that BMP15 played a role in the suppression of Fshr expression (Otsuka et al., 2001;McMahon et al., 2008;Abir and Fisch, 2011;Shimizu et al., 2019). This discrepancy was probably due to species difference or the abnormal development of TGF GCs, including degradation and premature luteinization. In addition, mRNA level of mitotic spindle assembly checkpoint protein (Mad2l1) decreased significantly in TGF ovarian tissues and antrum follicles, implying the inhibition of TGF cell mitosis. Increased expression of Mkp1 (Dual specificity protein phosphatase 1) in 60day old TGF ovarian tissues but decreased expression of this gene in 60-day TGF ovarian tissues and antrum follicles ( Figure 7C) seemed to be related to the up-regulated MAPK pathway in early TGF preantral follicles ( Figure 7D).

Knockdown of BMP15 Caused Reduced Capacity of TGF Follicles to Ovulate
Evidence of disordered estrous cycle, abnormally enlarged AFs, and no corpus lutein observed in sexually mature TG gilts until 365-day old demonstrate that knockdown of Bmp15 could cause dysovulation. To investigate the underlying factors causing dysovulation, COCs were isolated from antrum follicles with a diameter of 5-7 mm for single-cell RNA sequencing. As expected, sequencing results showed a drastic reduction in BMP15 in TGF COCs (Table 2), which was confirmed by qPCR analysis in antrum follicles ( Figure 7C). Interestingly, GDF9, the closely related homologous protein of BMP15, was down-regulated with the knock down of BMP15 ( Figure 7C and Table 2). However, the expression of other BMP proteins (BMP4 and BMP6) was not affected by the knock down of BMP15. Surprisingly, BMP15 receptors (Bmpr2 and Alk6) as well as its signaling protein, Smad8, were significantly up-regulated, which was probably related to the increased expression of Fst and Inhα or another activin. About four-folds decreased expression of both Fshr and Hsd17β (Luu-The et al., 2006; Figure 8C and  Table  S6), which might have contributed to premature luteinization (Jakimiuk et al., 2001;Stocco, 2001;Magoffin, 2005;Long et al., 2009). It has been reported that the down-regulated expression of Amhr2 (Anti-Mullerian hormone receptor type 2) and Cx43 (Gap junction protein alpha 1) induced by LH in preovulatory follicles was important to ovulation (Norris et al., 2008;Pierre et al., 2013). Thus, the increased expression of these two genes in TGF COCs potentially resulted in a decreased capacity of oocyte meiosis resumption and ovulation. In addition, markedly decreased expression of oocyte quality related genes (Bmp15, Gdf9, Zp2, Zp3, Zar1, and Irf6) strongly implies a reduced oocyte competence in TGF COCs (Picton et al., 1998;Wu et al., 2003Wu et al., , 2007Hussein et al., 2006;Sabel et al., 2009). A total of 2,820 DEGs (885 up-regulated, 1,935 downregulated) was generated for pathway enrichment. The significantly up-regulated AMPK and ovarian steroidogenesis (Figures 8A-C) pathways were likely to contribute to the greater number of large antrum follicles in TGF ovaries, according to the findings of previous studies in sheep (Foroughinia et al., 2017) and sow (Knox, 2005). However, pathways (e.g., Cell cycle, P53) involved in cell proliferation and growth were significantly down-regulated (Figures 8A,B), which was consistent with the results of the dynamic transcriptomic analysis of TGF follicles (Supplementary Figure S6B). Furthermore, four pathways including oocyte meiosis, oocyte maturation, cell cycle, and ovarian steroidogenesis, which were closely involved in the regulation of oocyte maturation and ovulation, presented DEGs enrichment of more than 20% (Figures 8A,B). In total, these results revealed both impaired function of cumulus cells and oocyte maturation in TGF COCs, suggesting a reduced capacity to ovulate. However, surprisingly, the expression of inosine monophosphate dehydrogenase 2 (Impdh) and natriuretic peptide receptor 2 (Npr2) was not affected. These two genes have been reported in mice to be up-regulated by BMP15 and GDF9 during the activation of maturation promoting factor (MPF) [Cyclin B and Cyclin dependent kinase 1 (CDK1)] and stimulation of oocyte meiotic resumption in vitro . Instead, we discovered significantly decreased expression of Cyclin B and Cdk1 in TGF follicles from secondary stage onward (Figure 8D and Supplementary Table S8), which implies the involvement of BMP15 in the modulation of porcine oocyte meiosis possibly by regulating the expression of MPF.

DISCUSSION
The effect of BMP15 mutations on ovarian follicular development and ovulation rate was first discovered in Inverdale (FecX) sheep (Davis et al., 1992;Braw-Tal et al., 1993;Smith et al., 1997). In these sheep, ewes with single allele of inactive BMP15 gene showed increased ovulation rate and a higher incidence of twin or triplet births, whereas ewes with bi-alleles of inactive BMP15 gene were sterile with primary ovarian failure phenotype (Galloway et al., 2000). Studies on animals immunized with different regions of BMP15 peptide revealed that an increased ovulation rate can be found in females where some of the BMP15 have been neutralized, but an inhibition of follicular growth and ovulation was found in females where most of the active BMP15 has been neutralized (McNatty et al., 2007). Therefore, different extent of reduction of biologically active BMP15 protein level seems to have different influence on fertility. In this study, we found two different ovarian phenotypes (TGS and TGF) in our BMP15 knockdown gilts. The different ovarian phenotypes might be caused by different in vivo expression level of BMP15. Indeed, in TGS ovaries, the marked reduced level of both mRNA ( Figure 1E) and protein ( Figure 5A) of BMP15 revealed that FIGURE 9 | Summary of biological functions and possible regulatory mechanism of BMP15 in TGF follicular development. Compared to WT gilt, knockdown of BMP15 caused remarkable reduction of follicle number during the TGF follicular developmental process, accompanied by impaired oocyte quality, degradation, and premature luteinization in TGF preantral follicles. This may lead to increased expression of LHR and dramatically decreased E2 production in TGF antral follicles, resulting in lack of dominant follicle selection but abnormally enlarged antral follicles, presenting a dysovulation phenotype. Decreased E2 production in TGF antral follicles may have a negative feedback to pituitary to increase the expression of Fsh. However, the decreased expression of Fsh receptor, Fshr, in GCs beyond the PF stage, by knocking down BMP15, attenuates the stimulation of antral follicles growth. In addition, BMP15 knockdown leads to the inhibition of cell proliferation and differentiation, declined oocyte quality, and meiotic maturation during TGF follicular development, contributing to the appearance of enlarged follicle in TGF ovary and dysovulation. Colored ellipses represent results based on morphological observation and molecular examination. Dotted ellipses represent results based on transcriptomic analysis. the majority of BMP15 was knocked down by integrated shRNA. In contrast, BMP15 protein accumulated abundantly in TGF ovaries, and was only slightly lower than that in WT ovaries, despite the fact that the mRNA level of BMP15 had decreased to half of that of the wild-type as detected in 365-day ( Figure 1E) and 30-day (Figures 1F,G) old ovarian tissues. Therefore, the low interference of BMP15 expression in TGF ovaries may confer on them a less severely impaired ovarian phenotype, as TGF ovaries contained each stage of follicle but presented remarkably reduced follicle number, increased ratio of abnormal follicle, disordered reproductive hormones, and ovulation dysfunction. We found that TGF and TGS ovaries could be concurrently present in a single TG gilt (Figure 3A), and this may be caused by epigenetic mechanisms including DNA methylation and histone modifications (Spencer et al., 2015), which could cause dynamic and heterogeneous expression of shRNA from integrated TG construct in different individuals and even in different ovaries within a single individual, explaining the varied in vivo levels of active BMP15 in each ovary of each individual. We suspect that in some critical stages of early follicle development, there may be a threshold level of active BMP15 that supports the development of follicles into subsequent stages, indicating that the lack of BMP15 in a specific ovary may hamper the transition of the follicles into the subsequent developmental stages. Previous reports in sheep revealed that heterozygous mutations in BMP15 inhibit GCs growth but increase GCs sensitivity to FSH, leading to increased ovulation in smaller matured follicles with reduced amounts of E2 and inhibin (Shackell et al., 1993;Otsuka et al., 2000Otsuka et al., , 2001Fabre et al., 2006). However, this was inconsistent with our results. The in vivo mRNA level of BMP15 in TGF ovaries was knocked down to half of that in wild-type, and thus, TG gilts with TGF ovaries could be considered as pigs with heterozygous mutations in BMP15. In contrast to ewes heterozygous for mutations in BMP15, our TG gilts presented impaired ovulation, as corpus lutein could not be found in TGF ovaries from TG gilts younger than 365 days, though they can be observed in TGF ovaries from 400-and 500-day old TGF gilts (Figure 4E), indicating a delayed ovulation in TGF gilts. In addition, significantly decreased amount of total and smaller AFs in TGF ovaries (Figures 4A,B) seemed incapable of supporting an increased ovulation rate. Moreover, the appearance of abnormally enlarged AFs with lower FSHR expression level and E2 production ( Figure 4D) may indicate that TGF follicles could not ovulate at normal size probably due to the attenuated FSH sensitivity and a lack of dominant follicle selection in TG gilts. Besides, premature luteinization of GCs possibly caused by insufficient FSH stimulation may also contribute to the dysovulation of TGF gilts. Thus, the results in TGF puberty gilts were different from that of the poly-ovulatory Bmp15 −/− mice, which showed normal follicular development and could ovulate at puberty (Yan et al., 2001). Hence, we suggest that the role of BMP15 in ovulation is species-specific and different not only between mono-and poly-ovulatory mammals, but also between poly-ovulatory species.
Bone morphogenetic protein 15 suppresses Fshr expression in GCs, which affects GCs proliferation and steroidogenesis in the AFs of rodents (Otsuka et al., 2001;McMahon et al., 2008) and humans (Abir and Fisch, 2011;Shimizu et al., 2019). Previous studies on sheep indicated that heterozygous mutations in BMP15 could increase the sensitivity of GCs to FSH stimulation in the AF, leading to increased ovulation rate (Fabre et al., 2006). However, a recent study showed that treatment with BMP15 caused increased expression of Fshr in bovine preantral follicles after 12 days of culture (Passos et al., 2013). These conflicting reports are probably caused by species-specific differences and the different response to BMP15 stimulation in each follicular development stages. In this study, we found that, in contrast to the results found in sheep, the knockdown of BMP15 did not increase but significantly inhibited Fshr expression in both preantral and AFs ( Figure 7C, Table 2 and Supplementary  Table S6). Furthermore, the inhibition in GCs proliferation and differentiation, increased expression of genes involved in steroidogenesis (StAR, Cyp11a, and 3βHSD) (Figures 5B, 8C, Table 2 and Supplementary Table S6), drastic reduction of E2 production (Figure 4E), and subsequent absence of dominant follicle selection in the TGF follicles, were likely the consequences of the declined sensitivity of GCs to FSH. Suppression of FSHR expression in preantral follicles of TG gilts implies that BMP15 could stimulate porcine follicle growth and development at an earlier follicle stage rather than at gonadotropin-dependent period (Mori, 2016). Given the degradation of GCs and abnormal structure of GCs layers observed in TGF ovaries, another possible reason for the decreased expression of FSHR might be related to the impaired development of GCs caused by BMP15 deficiency.
As a paracrine and autocrine factor of oocyte, BMP15 can promote not only the development of follicular somatic cells, but also the development of oocyte itself. Studies on the in vitro maturation of oocyte have demonstrated that BMP15 are able to stimulate cumulus cell expansion (Braw-Tal et al., 1993;Sugiura et al., 2010;Peng et al., 2013;Lin et al., 2014;Sudiman et al., 2014), promote the signaling of LH-induced maturation of the cumulus-oocyte complex , improve oocyte quality (Hussein et al., 2006;Caixeta et al., 2013), and increase blastocyst rate and embryonic development of fertilized oocytes (Wu et al., 2007;Gode et al., 2011). Recently, the expression level of BMP15 has been suggested as a diagnostic marker of oocyte quality (Wu et al., 2007). In this study, we confirmed the important role of BMP15 in porcine oocyte development in vivo. We found that the knockdown of BMP15 could cause oocyte degeneration or abnormal enlargement (Figures 3G,H), and lack of normal autophagy activity in oocytes of abnormal preantral follicles ( Figure 5C). Further transcriptomic analysis of follicle and COCs also implies that genes and pathways involved in oocyte meiosis and maturation in TGF follicles were affected from secondary stage onward (Figures 7D, 8A). Previous studies have indicated the possible mechanisms of action of BMP15 in the regulation of oocyte meiosis. One study reported that BMP15 and GDF9 can promote oocyte meiotic resumption in mice through the up-regulation of Npr2 and Impdh . Another study showed that inhibiting BMP15 signaling pathway by Smad2/3 phosphorylation inhibitor resulted in significantly decreased expression of Cdc2 and Cyclinb1 during porcine oocyte in vitro maturation (Lin et al., 2014). However, our transcriptomic results showed that the expression level of Npr2 and Impdh were not affected in both TGF follicles and COCs, instead, the expression levels of MPF (Cdc2 and Cyclinb1) decreased significantly (Figure 8D and Supplementary Table  S8). Therefore, our findings demonstrate that the regulation of MPF expression is associated with the role of BMP15 in porcine oocyte meiosis and maturation; however, this requires further elucidation.
In summary, the knockdown of BMP15 caused markedly reduced fertility of TG gilts mainly through the inhibition of both GCs and oocyte development (Figure 9). The suppression of GCs proliferation and differentiation led to decline in the number of early follicles, GCs degradation, and reduced sensitivity of GCs to FSH stimulation, consequently leading to premature luteinization, higher LHR expression, and lower E2 production in large AFs. The effect on oocyte development directly led to impaired oocyte quality and oocyte meiotic maturation. Consequently, large AFs become abnormally enlarged, resulting in dysovulation and disordered reproductive hormones. Our results revealed a remarkable physiological suppression of porcine ovarian follicular development and ovulation in BMP15 knockdown gilts, demonstrating an essential role of BMP15 on porcine female reproduction and providing new insights into the regulatory role of BMP15 in poly-ovulatory mammals. Our findings provide the basis for further investigation on the complicated regulatory function of BMP15 in female fertility of poly-ovulatory species and the development of possible strategies for improving porcine female fertility through the modulation of BMP15 expression.

DATA AVAILABILITY STATEMENT
The authors confirm that all data underlying the findings are fully available without restriction. Data are available from the Sequence Read Archive (SRA) under Accession no. PRJNA561679.

ETHICS STATEMENT
All procedures were performed in strict accordance with the recommendations of the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The protocol was approved by the Institutional Animal Care and Use Committee (IACUC), Sun Yat-sen University (approval Number: IACUC-DD-16-0901).

AUTHOR CONTRIBUTIONS
ZH, YC, and XHL conceived and supervised the project. YQ, ZL, XY, and ZH contributed to the transgenic gilts generation. YQ, TT, WL, XS, GS, XFL, MW, and XYL contributed to the tissues and blood sample collection. YQ, TT, and WL designed and performed most of the molecular experiments. YQ contributed to the analysis of RNA-seq data. YC, XHL, PC, and DM provided support and technical assistance. YQ and ZH wrote the manuscript. YC provided critical revision of the manuscript. All authors read and approved the final version of the manuscript.

FUNDING
This work was jointly supported by the National Transgenic Major Program of China (2016ZX08006003-006) for generation of pig models, National Key R&D Programmes of China (2018YFD0501200), and Natural Science Foundation of Guangdong Province (2016A030313310) for investigation of the roles of BMP15 in porcine ovarian follicular development.