Transcriptomic Analysis of STAT1/3 in the Goat Endometrium During Embryo Implantation

Interferon tau (IFNT), a pregnancy recognition signal in ruminants, promotes the establishment of embryo implantation by inducing the expression of interferon-stimulated genes (ISGs) via the Janus kinase/signal transducer and activator of transcription (JAK/STAT) signaling pathway. However, the precise regulatory mechanism of IFNT in goat embryo implantation remains largely unknown. In this study, we performed RNA sequencing of goat endometrial epithelial cells (gEECs) with or without 20 ng/mL IFNT treatment. Differential comparison showed that there were 442 upregulated differentially expressed genes (DEGs) and 510 downregulated DEGs. Bioinformatic analyses revealed that DEGs were significantly enriched in immune-related functions or pathways. The qRT-PCR validation results showed that the expression levels of STAT family members (STAT1, STAT2, and STAT3) were significantly upregulated in gEECs after IFNT treatment, which is in agreement with the RNA-seq data. Meanwhile, the protein levels of p-STAT1 and p-STAT3 increased significantly in gEECs after 6 and 24 h of IFNT treatment, respectively. Further in vivo experiments also confirmed that both mRNA and protein phosphorylation levels of STAT1 and STAT3 in the uterus on day 18 of pregnancy (P18) were significantly increased compared to those on day 5 (P5) and day 15 of pregnancy (P15). On P5, STAT1 and STAT3 proteins were primarily located in the uterine luminal epithelium (LE) and glandular epithelium (GE), and were also detected in the stromal cells. The intense immunostaining of STAT1 and STAT3 proteins were decreased on P15 and then increased on P18, especially in the superficial GE and subepithelial stromal cells. Moreover, p-STAT1 and p-STAT3 were highly expressed in the deep GE on P18. Collectively, these results highlight the role of IFNT in regulating endometrial receptivity in gEECs and uncover the temporal and spatial changes in the expression of STAT1/3 during embryo implantation in the goat endometrium.

Interferon tau (IFNT), a pregnancy recognition signal in ruminants, promotes the establishment of embryo implantation by inducing the expression of interferon-stimulated genes (ISGs) via the Janus kinase/signal transducer and activator of transcription (JAK/STAT) signaling pathway. However, the precise regulatory mechanism of IFNT in goat embryo implantation remains largely unknown. In this study, we performed RNA sequencing of goat endometrial epithelial cells (gEECs) with or without 20 ng/mL IFNT treatment. Differential comparison showed that there were 442 upregulated differentially expressed genes (DEGs) and 510 downregulated DEGs. Bioinformatic analyses revealed that DEGs were significantly enriched in immune-related functions or pathways. The qRT-PCR validation results showed that the expression levels of STAT family members (STAT1, STAT2, and STAT3) were significantly upregulated in gEECs after IFNT treatment, which is in agreement with the RNA-seq data. Meanwhile, the protein levels of p-STAT1 and p-STAT3 increased significantly in gEECs after 6 and 24 h of IFNT treatment, respectively. Further in vivo experiments also confirmed that both mRNA and protein phosphorylation levels of STAT1 and STAT3 in the uterus on day 18 of pregnancy (P18) were significantly increased compared to those on day 5 (P5) and day 15 of pregnancy (P15). On P5, STAT1 and STAT3 proteins were primarily located in the uterine luminal epithelium (LE) and glandular epithelium (GE), and were also detected in the stromal cells. The intense immunostaining of STAT1 and STAT3 proteins were decreased on P15 and then increased on P18, especially in the superficial GE and subepithelial stromal cells. Moreover, p-STAT1 and p-STAT3 were highly expressed in the deep GE on P18. Collectively, these results highlight the role of IFNT in regulating endometrial receptivity in gEECs and uncover the temporal and spatial changes in the expression of STAT1/3 during embryo implantation in the goat endometrium.
Keywords: goat, interferon tau, embryo implantation, RNA-Seq, STAT1, STAT3 INTRODUCTION Embryo implantation is crucial for successful pregnancy. During this period, the elongated conceptus and receptive endometrium recognize each other and have complex interactions, which change the expression of a number of genes, resulting in promotion of embryo implantation and placenta formation, eventually leading to the establishment of pregnancy (1,2). Implantation failure is the major cause of pregnancy loss in cattle, accounting for 30%-50% of all cases (3, 4). Although many studies have been conducted to identify the gene network in the endometrium or conceptus during embryo implantation (5)(6)(7)(8), the precise molecular mechanism has not been well-characterized.
Interferon tau (IFNT) is a type I IFN secreted from trophectoderm cells of the ruminant conceptus, which is instrumental in the pregnancy recognition of ruminants in the synchronization of the conceptus and maternal uterus (9). In addition to its anti-luteolytic action, IFNT regulates the function of the receptive endometrium and the elongation of the conceptus to facilitate the establishment of pregnancy (10). Upon binding to the receptor, IFNT performs these functions through the Janus kinase/signal transducer and activator of transcription (JAK/STAT) signaling pathway. In the JAK/STAT signal transduction process, JAK1 and Tyrosine kinase 2 (Tyk2) successively phosphorylate STAT1 and STAT2 (11,12). The phosphorylated-STAT1/STAT2 heterodimer forms the interferon-stimulated gene factor 3 (ISGF3) complex with interferon regulatory factor 9 (IRF9) and translocates to the nucleus, which binds to IFN-stimulated response elements (ISREs) in the promoter region of a group of interferon-stimulated genes (ISGs), leading to the transcription of ISGs (13). In addition, type I IFN induces STAT3 phosphorylation (14). Similar to other STATs, STAT3 forms complexes with other transcription factors, including the STAT3 homodimer (15) and STAT1/STAT3 heterodimer (16), and translocates to the nucleus for signal transduction.
Notably, the mRNA levels of critical signaling components of the JAK/STAT pathway (STAT1, STAT2, and IRF9) are low in the endometrial luminal epithelium (LE) during the peri-implantation period in sheep (17). In cattle, nuclear STAT1 expression is reduced in the LE during the peri-implantation period when compared with that during the pre-implantation period, indicating a decrease in phosphorylated STAT1 protein levels and subsequent inhibition of ISGF3 synthesis (18). However, a detailed analysis of the expression profile of STATs in the goat uterus during embryo implantation has not been reported. Additionally, although the JAK/STAT3 signaling pathway is one of the key pathways in decidualization (19,20), its role in ruminants remains unknown.
Therefore, to explore the potential function of JAK/STAT3 during early pregnancy in goats, based on the gene regulatory network of IFNT-treated goat endometrial epithelial cells (gEECs) by RNA sequencing, we determined the expression and localization of STATs in the goat endometrium during the peri-implantation period.

Tissue Collection
Mature Guanzhong dairy goats (n = 9, aged 2-3 years, average weight = 59.28±1.93 kg) were reared in the experimental animal center of Northwest A & F University, Yangling, China. The goats exhibiting at least two estrous cycles of normal duration were used in this study. At estrus, female goats were mated with fertile males to induce natural pregnancy, which was recorded as day 0 of pregnancy. Pregnancy was confirmed on day 5 (P5, in the pre-implantation period) by recovering blastocysts from the uterus. Pregnancy at day 15 (P15, in the pregnancy recognition period) and day 18 (P18, in the embryo adhesion period) was identified by observing the elongated tubular conceptus and linear conceptus in the uterus, respectively. The uteri of pregnant goats on P5 (n = 3), P15 (n = 3), and P18 (n = 3) were collected immediately after the goats were subjected to midventral laparotomy and hysterectomy. The uterine tissues were fixed in 4% (v/v) paraformaldehyde in phosphate-buffered saline (PBS) without Ca 2+ /Mg 2+ (PBS) or immediately frozen in liquid nitrogen.

Cell Culture and Treatment
Immortalized gEECs were obtained as previously described (21). The gEECs were cultured in Dulbecco's minimum essential medium nutrient mixture F-12 (DMEM/F-12, BasalMedia, China) supplemented with 10% fetal bovine serum (FBS, ZETA LIFE, USA) at 37 • C in a humidified atmosphere of 5% CO 2 . When gEECs reached 50% confluence, the cells treated with or without 20 ng/mL recombinant ovine IFNT (Sangon Biotech) for 6 h were used for RNA-seq analysis. Meanwhile, gEECs treated with or without 20 ng/mL IFNT for 6, 12, and 24 h were used for qRT-PCR confirmation.

RNA Isolation and Quality Assessment
Total RNA of gEECs or uterus tissue was isolated using RNAiso Plus (Takara, Japan), and genomic DNA was removed using the Evo M-MLV RT Kit (Accurate Biology, China) according to the manufacturer's protocols. The purity, concentration, and integrity of the total RNA were checked using a NanoPhotometer R spectrophotometer (IMPLEN, CA, USA), Qubit R RNA Assay Kit in Qubit R 2.0 Fluorometer (Life Technologies, CA, USA), and RNA Nano 6000 Assay Kit of the Bioanalyzer 2,100 system (Agilent Technologies, CA, USA).

RNA Sequencing
Library preparation, clustering, and sequencing for RNA sequencing were performed according to the manufacturer's recommendations. Briefly, ribosomal RNA was removed from 3 µg of total RNA using the Epicenter Ribo-zero TM rRNA Removal Kit (Epicenter, USA). Subsequently, RNA-seq libraries were generated using the Library Prep Kit for Illumina R (NEB, USA). The library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, USA). Next, 3 µL of USER Enzyme (NEB, USA) was used with adaptor-ligated cDNA before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. The products were then purified using the AMPure XP system and library quality was assessed using the Agilent Bioanalyzer 2,100 system. Clustering was performed on a cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina). After cluster generation, the RNA-seq libraries were sequenced on an Illumina HiSeq 4,000 platform. The sequencing data generated in this study were deposited in the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE184110).

RNA-seq Data Analysis
Raw reads were aligned to the goat genome using HISA T2 (v2.0.4) software. Reference genome and gene model annotation files were downloaded directly from the genome website. Transcriptome assembly was performed using StringTie (v1.3.1) software in a reference-based approach. Gene expression abundance was assessed by calculating the FPKMs of genes in each sample using Cufflinks (v2.1.1) software. Differentially expressed gene analysis was performed using the Ballgown suite. Transcripts with an adjusted p-value <0.05 were assigned as "differentially expressed." P-values were adjusted using the Benjamini-Hochberg procedure.

GO and KEGG Enrichment Analyses
Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) was performed using the GOseq R package. KEGG enrichment analysis of DEGs was performed using the KOBAS software. GO terms and KEGG with q-values <0.05 were considered significantly enriched by DEGs.

qRT-PCR
Reverse transcription was performed using the Evo M-MLV RT Kit (with gDNA Eraser, Accurate Biology, China) according to the manufacturer's protocol. PCR was performed using the ChamQ SYBR qPCR Master Mix (Vazyme, China). Data collection and analysis were performed on a CFX Connect machine (Bio-Rad, USA) using CFX Manager software. All quantitative PCRs were performed in triplicate. Primers were synthesized by Tsingke Biotechnology Co., Ltd. (China), and their sequences are listed in Supplementary Table 1. Relative mRNA expression levels of target genes were normalized against the relative quantity of GAPDH mRNA in the same sample, and calculated by the formula using the 2 − Ct method.

Western Blotting
The gEECs and uterus tissues were lysed using a protein extraction kit (KeyGEN BioTECH, China) to obtain total protein. After 10 min of lysis at 4 • C, the samples were centrifuged at 12,000 g at 4 • C for 10 min to collect the supernatant for concentration determination. The concentration of protein samples was determined using a BCA kit (KeyGEN BioTECH, China). Next, the protein samples were added to loading buffer (5×, CWBIO, China) and heated at 100 • C for 10 min using a DryBath (Thermo Fisher Scientific, USA). For western blotting, the protein samples were separated by SDS-polyacrylamide gel electrophoresis (SDS-PAGE) on gels and blocked with 5% nonfat milk (BBI, China) in Tris-buffered saline Tween (TBST)

Immunohistochemistry
The uterine tissues were fixed for 24 h, dehydrated through a graded ethanol series, and embedded in paraffin. Next, 5-µmthick sections were mounted onto glass slides pre-coated with poly L-lysine solution (Sigma, USA) and incubated overnight at 37 • C. After deparaffinization and rehydration, sections were placed in citrate buffer (pH 6.0). Antigen retrieval was performed by heating the sections for 20 min in a microwave oven at 95 • C, and the slides were washed in PBS. The sections were then stained with the streptavidin-peroxidase method using an UltraSensitiveTM SP (Rabbit or Mouse) IHC Kit (Maixin Biotechnologies, China). Briefly, the sections were pretreated with 0.3% (v/v) H 2 O 2 for 40 min at 37 • C to quench endogenous peroxidase activity. After washing with PBS, the sections were incubated with 10% pre-immune serum for 70 min at 37 • C. After blocking, the sections were incubated with anti-STAT1 (1: 800, diluted in PBS), anti-Phospho-STAT1 (Tyr701, 1: 100), anti-STAT3 (1: 300), and anti-Phospho-STAT3 (Tyr705, 1: 100) at 4 • C overnight, and rewarmed at 37 • C for 40 min, washed with PBS, and incubated with biotin-labeled goat anti-rabbit or anti-mouse IgG at 25 • C for 40 min. The sections were then washed with PBS and incubated with streptavidin-biotin peroxidase for 40 min at 25 • C. Thereafter, the sections were visualized with diaminobenzidine (DAB), lightly counterstained with hematoxylin for 25 s, dehydrated, and coverslipped. As a negative control, the primary antibody was substituted with pre-immune serum. Sections were imaged under a microscope (Nikon, Germany) after drying at 25 • C.

Statistical Analysis
Data are presented as the mean ± standard error of the mean (SEM). Statistical significance was set at p < 0.05. All data were representative of at least three different experiments and statistically analyzed using SPSS software (SPSS 20.0, IBM Corp., Armonk, NY, USA). Student's t-test was used for comparisons between the two groups. Statistical differences among different groups (>2) were evaluated using one-way ANOVA with multiple comparisons among groups tested by Tukey's post hoc test.

Differentially Expressed Genes Analysis
We generated triplicate libraries from IFNT-treated and untreated gEECs (control group) and sequenced them to a depth of 89.46 Gb clean data. Alignment of the clean data against the reference genome yielded 82.79-89.78% of the total mapped reads, >90.70% of the Q30, and exceeding 46.38% of the GC content. Multiple mapped reads (3.88-5.46%) were excluded from further analyses. Therefore, 43,779 mRNA transcripts were detected in the 6 samples, indicating that the RNA sequencing results were reliable. After the quantification of gene expression levels, a volcano plot was used to provide an overview of the DEGs (Figure 1A). A total of 952 genes were differentially expressed between the IFNT and control groups. Among these genes, 442 genes were upregulated and 510 were downregulated (Supplementary Table 2). A heatmap depicting the DEGs is shown in Figure 1B.

GO and KEGG Enrichment Analyses
GO and KEGG enrichment analyses were performed to identify the potential biological processes of the DEGs. In the GO analysis performed among IFNT and control groups, the number of upregulated DEGs annotated to GO terms was identified (downregulated DEGs annotated to a few GO terms), and the significantly enriched GO term was covered. The top 23 enriched functions, including three biological processes and 20 molecular functions (Figure 2). Protein ubiquitination (GO: 0016567), protein modification by small protein conjugation (GO: 0032446), and response to biotic stimulus (GO: 0009607) were the primary categories in biological processes. Hydrolase activity (GO: 0016787) and small molecule binding (GO: 0036094) were the main function categories.
The DEGs pathway analysis was based on the KEGG pathway database for calculation of the abundant pathways. The top 20 enriched pathways of downregulated and upregulated DEGs are shown in Figures 3A,B. Among the downregulated DEG enrichment pathways, pathways in cancer showed the highest degree of enrichment. Among the upregulated DEGenriched pathways, the main enriched pathways were those related to influenza A, herpes simplex infection, measles, and hepatitis C.

STATs Expression in gEECs by IFNT Treatment
In agreement with the RNA-seq data, the expression levels of STAT1, STAT2, and STAT3 mRNA were significantly upregulated by IFNT in gEECs measured using qRT-PCR (Figure 4A; p < 0.05). Moreover, STAT1 protein levels in gEECs gradually increased with increasing time of IFNT treatment ( Figure 3B). Meanwhile, p-STAT1 expression was markedly increased by IFNT treatment and peaked at 6 h ( Figure 3B; p < 0.05). Although STAT3 protein levels in gEECs were not significantly altered within 24 h of IFNT treatment, p-STAT3 levels were significantly increased after 24 h of IFNT treatment (Figure 4C;   To validate RNA-seq data, the mRNA level of STAT1, STAT2, and STAT3, which was detected by qRT-PCR after gEECs were treated with 20 ng/mL IFNT for 6 h. Gene expression was normalized to GAPDH. Values represent as the mean ± S.E.M. from three independent experiments, and compared with the control group. (B,C) The protein expression of STAT1, p-STAT1, STAT3, and p-STAT3 was analyzed using Western Blotting after 20 ng/mL IFNT treatment for 6, 12, and 24 h. The data are presented as the mean ± S.E.M. from three independent experiments, and bars with different letters (a/b/c) are significantly different (p < 0.05).
p < 0.05). Unfortunately, the STAT2 and p-STAT2 antibodies which are suitable for goat tissue, were not detected in this study.

Expression Analysis of STATs in the Goat Uterus During Embryo Implantation
As shown in Figure 5A, the mRNA levels of STAT1, STAT2, and STAT3 in the uterus were higher on P18 than on P5 and P15, respectively (p < 0.05). However, there was no significant difference in STAT1, STAT2, and STAT3 mRNA expression levels between P5 and P15. Furthermore, no obvious changes in the STAT1 and STAT3 proteins were detected among P5, P15, and P18, whereas p-STAT1 and p-STAT3 were markedly upregulated on P18 (Figures 5B,C; p < 0.05).

Localization of STAT1/3 in the Goat Uterus
During Embryo Implantation STAT1 protein was mainly localized in the LE and glandular epithelium (GE) on P5, and slight immunostaining of the STAT1 protein was detected in the stromal cells (Figure 6). On P15, weak STAT1 immunostaining was only weakly observed in superficial GE (sGE) and subepithelial stromal cells. STAT1 expression was significantly increased in endometrial cells and was particularly abundant in the GE and subepithelial stromal cells on P18. Meanwhile, no obvious immunostaining of p-STAT1 protein was detected in the endometrium on P5. The p-STAT1 protein was slightly expressed in the sGE on P15; the expression gradually increased and was strongly detected in the deeper GE on P18. The STAT3 protein was widely distributed in the endometrial cells on P5, P15, and P18 (Figure 7). Moreover, p-STAT3 was only expressed in the GE at P18 (Figure 7).

DISCUSSION
In this study, we first provide a relatively comprehensive DEG profile for gEECs. Based on the RNA-seq results, we found that the expression and localization of STAT1 and STAT3, which are downstream components of the JAK-STAT signaling pathway, showed spatiotemporal changes during embryo implantation process.
By comparing the previously reported transcriptome analysis with the DEGs in our study, we observed that few genes were uniformly regulated by IFNT between EECs and in vivo uteri. In the co-culture system of bovine EECs and stromal cells, the  authors identified 222 upregulated DEGs in bovine EECs by comparing IFNT-treated and untreated groups; of which 96 of 222 appeared in our DEG list and all showed upregulation (22). We also compared our data with the DEGs in the receptive uterus vs. unreceptive uterus reported by Zhang et al. (23), because there are few reports on the transcriptome of IFNTtreated cells. There are 5 of 73 upregulated and 53 of 737 downregulated DEGs in the goat receptive uterus vs. unreceptive uterus were matched to our DEG list, and only 2 of 5 and 42 of 53 of upregulated and downregulated DEGs, respectively, were consistent with the gene regulation trends in our data (23). Similarly, in an earlier report, Gray et al. used microarray analysis to identify 45 genes specifically regulated by IFNT in the ovine uterus, but it only yielded 10 of 45 genes, consistent with our data (24).
Using bioinformatics analysis, we found that the significant enrichment and higher number of GO terms induced by IFNT were involved mainly in the regulation of enzyme activity, protein modification, and defense responses. Notably, ISG15 is involved in the formation of nuclease activity and ubiquitin-like modification, which once executes its functional interaction with the DNA helicase and substrate protein and regulates the immune response by altering the genome stability and protein function (25,26). However, in vivo studies have shown that energy metabolism is the most important enriched function in the ruminant receptive uterus (23,27), and were also predicted to potentially improve pregnancy outcomes (28). KEGG analysis in this study also demonstrated that defense response pathways, such as RIG-I (DDX58) like receptor, NF-κB, and JAK/STAT pathway, were remarkably enriched pathway, which is similar to the results of GO analysis.
It has been previously reported that the JAK/STAT signaling pathway mediates the formation of endometrial receptivity through the induction of IFNT (29). Therefore, we selected the JAK/STAT signaling pathway among the DEGs enriched for validation. In line with these expectations, IFNT markedly induced the transcription of STATs and phosphorylation of STAT1. However, the expression of p-STAT3 was negatively regulated by IFNT from 6 to 24 h, which is consistent with the maintenance time of the phosphorylation state of STAT3 (30,31). Previous in vivo studies on STAT expression were limited to transcriptional and unphosphorylated protein levels (17,18). Therefore, we focused on the localization and expression of p-STATs in the goat endometrium. As described in the literature, the expression of classical ISGs in LE/sGE is inhibited by interferon regulatory factor 2 (IRF2) (31). Similarly, STAT1 and p-STAT1 were hardly expressed in LE/sGE cells on P15. However, in the GE, STAT1 and p-STAT1 were completely undetectable, which is contrary to previous studies reporting that classical ISGs are strongly expressed in the ovine endometrial GE during the peri-implantation period (32,33). The mechanism underlying this abnormality requires further study. Moreover, on P15, we observed that STAT3 was widely expressed in LE/GE/sGE cells, but we did not detect the expression of p-STAT3. This result may be related to the formation of unphosphorylated STAT3. Unphosphorylated STAT3 is also induced by IFN and can, in turn, induce the transcription of some classical ISGs through nuclear translocation (34).
In summary, these findings provide the first evidence of transcriptome changes in STATs via IFNT stimulation in vitro. We demonstrated that the distribution of STAT1 and STAT3 changes temporally and spatially in the goat endometrium during embryo implantation and that STAT1 and STAT3 may exert different effects on goat embryo implantation.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The experimental and surgical procedures were conducted according to protocols approved by the Committee for the Ethics on Animal Care and Experiments at Northwest A&F University.

AUTHOR CONTRIBUTIONS
HL and PL designed the study. HL, CW, and ZL performed the experiments. HL and ZL analyzed the data. HL wrote the manuscript and revised by PL, AW, and YJ. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to thank Editage (www.editage.cn) for English language editing.