Fresh and healthy uteri without endometritis from Holstein cows were acquired from a surrounding abattoir (Wuhan, China). The health and estimated estrous cycle of bovine were objectively confirmed by individual records of well-being and examination of veterinary surgeon from Wuhan Academy of Agricultural Sciences. This experimental processes were approved by the Institutional Animal Use Committee of Institute of Animal Science and Veterinary Medicine, Wuhan Academy of Agricultural Sciences.

Primary bovine endometrial epithelial cells were isolated and cultured in accordance with protocols previously described with some minor changes (24, 25). The uteri were eviscerated within 30 min after exsanguination, and immediately rinsed 3 times with ice-cold PBS to remove residual blood and dirt on the surface. After be immersed in sterile PBS, the uteri were kept on ice and transported to the laboratory within 1 h. The endometrium tissue was dissected into strips, which were placed in serum-free DMEM/F12 (Gibco, United States) supplemented with 50 U/L penicillin, 50 μg/mL streptomycin (Gibco, United States), 50 μg/mL gentamycin (Sangon Biotech, Shanghai, China) under sterile conditions. Then, the strips were cut into about 1mm³ pieces and washed 3 times with sterile PBS. Tissue pieces were digested with PBS containing 0.1% collagenase I (Sigma, United Kingdom) and incubated in a constant-temperature shaker at 80 r/min for 1 h at 37°C. The suspension was filtered with a 40 μm mesh, the filtrate was centrifuged to remove the supernatant at 1000 r/min for 10 min at room temperature, and resuspended in DMEM/F12 supplemented with 10% FBS (Gibco, United States). After washing five times, the cell pellets were resuspended with complete culture medium composed of 90% DMEM/F12, 10% FBS, 50 U/L penicillin and 50 μg/mL streptomycin, and dispersed into single cell suspension. Cells were seeded into a 10 cm culture dish which was kept in a water-jacked incubator at 37°C with 5% CO₂. The medium was changed every 48 h. Subculturing was conducted when the cells reached approximately 80–90% confluence. Cell morphology was recorded and photographed under an inverted microscope.

The endometrial epithelial cells were purified by differential attachment method following three to five passages. Then, the purity of epithelial cell culture was confirmed by immunofluorescence detection of cytokeratin 18. Briefly, the primary endometrial epithelial cells were plated on glass coverslips in six-well culture dishes. The cells were fixed in 4% paraformaldehyde for 15 min after reaching 50–60% confluence, and permeabilized with 0.5% Triton X-100 (prepared in PBS) for 20 min at room temperature. Non-specific binding was blocked by pre-treating the cells with 3% BSA (Sigma, United States) for 30 min at room temperature. Subsequently, the cells were incubated with primary antibody CK-18 (Servicebio, Wuhan, China) at 4°C overnight, followed by incubation with the secondary fluorescently-labeled Cy3 antibody (Servicebio, Wuhan, China) for 1 h at room temperature, after which the cells were washed three times in PBS. The nuclei were counterstained with DAPI (Thermo, United States). Fluorescent images were obtained using a laser scanning confocal microscope (NIKON Eclipse Ti, Tokyo, Japan).

As shown in Figure 1A, the BEND at logarithmic growth stage were inoculated in 6 cm dishes and cultured with DMEM/F12 containing 10% FBS and 1% penicillin–streptomycin in humidified atmosphere with 5% CO₂ at 37°C. Once reaching 50–60% confluence, the cells were divided into three groups of three replicates and rinsed with PBS. For the groups specified for E₂, P₄ and IFNT treatment, the medium were replaced by fresh complete culture medium supplemented with E₂ (10⁻⁹ M, Sigma) and P₄ (10⁻⁷ M, Sigma, United States), and continued culturing for 12 h. Subsequently, the cells were treated with IFNT (100 ng/mL, Cell sicences, United States) for an additional 12 h.

Total RNA extraction from BEND were performed with TRIzol reagent (Invitrogen, United States), and genomic DNA was removed using DNase I (Invitrogen, USA) according to the manufacturer’s instructions (HiPure Universal RNA Mini Kit, Magen). RNA degradation and contamination were monitored on 1% agarose gels. The quantity and purity of RNA was measured using Qubit 3.0 (Thermo Fisher Scientific, MA, United States) and Nanodrop One (Thermo Fisher Scientific, MA, USA). RNA integrity was accurately assessed using the Agilent 2,100 system (Agilent Technologies, Waldbron, Germany).

The mRNAseq library was generated using the ALFA-SEQ Directional RNaLib Prep Kit following the manufacturer’s recommendations. In brief, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was performed using RT Buffer interruption within the ALFA-SEQ Directional RNaLib Prep Kit. First-strand cDNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H), followed by second-strand cDNA synthesis using DNA polymerase I and RNase H. The remaining overhangs were converted into blunt ends by exonuclease/polymerase activity. After adenylation of the 3′ end of the DNA fragments, X-Fold diluted VAHTS Adapter-S with a hairpin loop structure were ligated to prepare for hybridization. To preferentially select cDNA fragments of 150–200 bp in length, these library fragments were selected by AMPure XP beads (GE Sera-Mag Magnetic Speed-beads carboxyl-modified). Then, PCR was performed with 2× PfuMax HiFi PCR ProMix (EnzyValley), primers with index VAHTS Multiplex Oligos Set 4 for Illumina. Finally, AMPure XP beads (GE Sera-Mag Magnetic Speed-beads carboxyl-modified) was used to purify the PCR products, and the library insert sizes were evaluated on the Qsep400 high-throughput nucleic acid protein analysis system (Houze Biotechnology Co., Ltd., Hangzhou, China). The clustering of the index-coded samples was performed on a cBot Cluster Generation System. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2,500 platform and 150 bp paired-end reads were generated. The sequencing data generated in this study were deposited in the NCBI GEO database (accession number GSE253481).

Use fastp (v0.21.0) (https://github.com/OpenGene/fastp) to process the raw data in fastq format to obtain clean data (clean reads). Clean reads were mapped to the NCBI Rfam database, and rRNA sequences were removed by Bowtie2 (v2.33) (https://github.com/BenLangmead/bowtie2). The reference genome (ARS-UCD1.2) and annotation files were downloaded directly from the genome website. The sequences after ribosome removal were aligned to the reference genome by Hisat2 (v2.2.1) software (https://github.com/ infphilo/ hisat2).

RSEM (v1.3.3) (https://github.com/deweylab/RSEM) was applied to obtain the read counts for each gene. Fragments per kilobase of exon per million mapped reads (FPKM) of genes in each sample were calculated to estimate gene expression levels. Differentially expressed genes between two groups (three biological replicates per group) were screened with DESeq2 package (v1.34.0) (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html). The resulting p-values were adjusted using the Benjamini-Hochberg method for controlling the false discover rate (FDR). |log2 (Fold Change)| > 1 and FDR < 0.05 were set as thresholds for significantly differential expressed genes (DEGs). These genes were used for subsequent bioinformatic analysis.

The function annotation of DEGs, including cellular component (CC), molecular function (MF) and biological process (BP), were analyzed using Gene Ontology (GO) database (http://www.geneontology.org). Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/) was engaged to systematically analyze the biochemical metabolic and signal transduction pathways in which the DEGs are involved. GO and KEGG enrichment analysis of DEGs were performed using the clusterProfiler (v4.2.2) (http://www.Bioconductor.org/packages/release/bioc/html/clusterProfiler.html) for function and signaling pathway annotation. GO terms and KEGG pathways with FDR < 0.05 were considered significantly enriched for DEGs.

A Protein–protein interaction (PPI) network was constructed for DEGs, based on the STRING database (http://string-db.org). The constructed networks were based on known interactions in the selected reference species (cow) with a combined interaction score above 0.4 as selection threshold. Cytoscape software (version 3.10.1) was applied to visualize the PPI networks.

The expression profiles of selected DEGs were validated to confirm the repeatability and reproducibility of gene expression data obtained by RNA sequencing, which were accomplished by SYBR Green-based qRT-PCR using specific primers designed by an online tool (https://primer3.ut.ee/) and a subset of samples from those utilized for RNA sequencing. 1 μg of RNA was reverse transcribed into cDNA, using a reverse transcription kit (with gDNA wiper, Vazyme, China) following the manufacturer’s instructions. The relative expression level of genes were determined by qRT-PCR using the AceQ qPCR SYBR Green Master Mix (Vazyme, China) with the following program: pre-degeneration at 95°C for 5 min, 40 cycles of degeneration at 95°C for 10s and annealing at 60°C for 30s. Data acquisition and analysis were performed on a CFX96 Connect real-time PCR detection system using CFX manager software (Bio-Rad, United States). The sequences of primers synthesized by Tsingke Biotechnology Co., Ltd. (China) are listed in Supplementary Table S1. Relative mRNA expression levels of target genes were normalized against the relative quantity of ACTB mRNA, which was based on threshold cycle (Ct) values using a comparative 2^-△△Ct method. qRT-PCR was carried out in three biological replicates.

RNA interference of DRAM1 was executed using small interfering RNA (siRNA). BEND was transiently transfected with control siRNA or siRNA-DRAM1 by Lipofectamine 3,000 (Invitrogen, United States) according to the manufacturer’s manual. The sequence of siRNA- DRAM1 was sub5’-UCCUGCAAUCCGUCAUCUCUUdTdT-3′. qPCR was employed to detect the knockdown efficiency of DRAM1 in BEND.

Experimental data are presented as mean ± standard deviation (SD). All data were representative for one of three independently repeated experiments. The statistical analyses were carried out by GraphPad Prism software 8.0 (GraphPad Prism Software Inc., CA, United States). Unpaired two-tailed t test (Student’s t-test) was devoted to make comparisons between two independent groups. For comparison of statistical differences among multiple groups (>2), one-way analysis of variance (one-way ANOVA) with multiple comparisons tested by Tukey’s post hoc test was performed. p < 0.05 was considered statistically significant.

The identification of primary BEND was performed by immunofluorescence, which showed that the epithelial-specific marker cytokeratin is positive in primary BEND (Supplementary Figure S1). The flow schema of sample preparation for sequencing is shown in Figure 1A. Triplicate mRNAseq libraries from three groups of BEND, which includes cells of control group (CON), cells treated with E₂ and P₄ (GA), cells treated with E₂, P₄ and IFNT (GB), were generated and sequenced to yield 65.15Gb of data. After filtering the sequencing adapter sequence and low-quality reads, a total of 62,289,909, 67,196,670, 63,332,621 clean reads with >93.67% of Q30 and exceeding 51.74% of the GC content were acquired from CON, GA and GB group, respectively, that were used for the succedent analyses. The mean value of clean reads ratio was 88.66, 88.25 and 89.55% for the CON, GA and GB group, respectively (Supplementary Table S2). Approximately 98.18–98.38% of clean reads were mapped to the reference genome after removal of reads aligned with rRNA (Supplementary Table S3). Unmapped reads (1.62–1.82%) were excluded from further analyses. Accordingly, 17,162 mRNA transcripts were detected in the 9 samples. Principal component analysis (PCA) which clusters similar samples together and accounts for the origin of variance in data was performed with two elements: PC1 (32.7%) and PC2 (17.3%), indicating high similarity between samples in the group (Figure 1B). All these results demonstrated the reliability of RNA sequencing data, that could be used in further analyses.

The differential gene expression profile between groups was determined with DESeq. A total of 608 DEGs were identified in GB group compared to the CON group, among which 409 genes were upregulated and 199 genes were downregulated; 529 DEGs composed of 396 upregulated genes and 133 downregulated genes were identified in GB group compared to the GA group; 62 DEGs containing 13 upregulated genes and 49 downregulated genes were detected in GA group compared to the CON group (Figure 1C). In addition, a venn diagram of DEGs in three comparison groups was plotted, which showed that 440 genes are co-differentially expressed in both comparison groups (GB vs. CON, GB vs. GA) (Figure 1D). The volcano plot was constructed to visualize the distribution and fold change of these DEGs identified in GB vs. CON and GB vs. GA (Figure 1E) groups.

The hierarchical cluster analysis of DEGs was applicated to identify the expression level of genes with high correlation among samples. Heatmaps depicting the DEGs of GB vs. CON, GB vs. GA groups are shown in Figures 2A,B, respectively. Meanwhile, we plotted a heatmap of 440 DEGs that were commonly present in both comparison groups, which revealed that most of the DEGs exhibit similar differential expression trends and patterns among samples (Figure 2C).

We performed gene ontology (GO) annotation of DEGs to characterize the functional roles of these genes in biological pathways. For GO enrichment analysis, GO terms with adjusted p-value <0.05 were designated as significantly enriched. A large number of DEGs were enriched in the GO biological process category. 17 GO terms from Process Ontology were significantly enriched in GB vs. CON group, such as immune system process, response to external biotic stimulus, defense response, response to cytokine, antigen processing and presentation (Supplementary Table S5). 43 GO terms from Process Ontology were significantly enriched in GB vs. GA group, including immune system process, multi-organism process, regulation of immune system process, regulation of response to stress, negative regulation of cellular protein metabolic process (Supplementary Table S6). The GO terms with most significant enrichment in GB vs. CON (Figure 3A), GB vs. GA (Figure 3B) are displayed as circle diagrams.

KEGG enrichment analysis of DEGs was conducted to calculate the abundant metabolic pathways and signal transduction pathways. KEGG pathways were considered significantly enriched with adjusted p-value <0.05. Overall, 348 DEGs were assigned to 24 significantly enriched KEGG pathways in GB vs. CON group, which contains herpes simplex infection, NOD-like receptor signaling pathway, influenza A, antigen processing and presentation, measles, parkinson disease, cytosolic DNA-sensing pathway, oxidative phosphorylation, RIG-I-like receptor signaling pathway, necroptosis (Supplementary Table S7). Moreover, 305 DEGs were allotted to 28 significantly enriched KEGG pathways in GB vs. GA group, such as herpes simplex infection, NOD-like receptor signaling pathway, human immunodeficiency virus 1 infection, human papillomavirus infection, influenza A, phagosome, cell adhesion molecules (CAMs), Toll-like receptor signaling pathway (Supplementary Table S8). Advanced bubble charts reveal the significantly enriched KEGG pathways of DEGs in GB vs. CON (Figure 3C), GB vs. GA (Figure 3D).

A PPI network composed of 260 nodes connected via 2018 edges was generated by Cytoscape software based on STRING database to illustrate the interactions between DEGs. The results showed that a few node genes interact closely with other genes, including STAT1, IFIH1, MX2, IFI44, USP18, DHX58, PSMB8, CASP8, DRAM1, CXCR4 (Figure 4). Some of these highly connected node genes are rarely reported to correlate with early pregnancy, that need to be concerned and further investigated.

Considering the results of GO and KEGG enrichment analysis, 8 DEGs were randomly selected for exploring the relative expression pattern to validate the accuracy and reliability of RNA-Seq results by performing qRT-PCR. The expression trend of all these selected genes in the results of qRT-PCR were proved to be identical with RNA-seq data (Figure 5A). Meanwhile, the results of Spearman’s correlation analysis of all selected genes showed a strong positive correlation between RNA-seq data and qRT-PCR experiments (r² > 0.65, p < 0.05), which confirmed the credibility of DEGs (Figure 5B).

To investigate the functional role of DRAM1 during early pregnancy in ruminants, we conducted RNA interference assay to downregulate the transcript expression of DRAM1 in BEND. The results indicated that the mRNA expression of DRAM1 is prominently reduced (Supplementary Figure S2). Additionally, the mRNA expression of HOXA11, PTGS1 and PTGS2 were dramatically downregulated by silencing DRAM1 upon IFNT stimulation (Figures 6B–D), but no significant change was observed on mRNA expression of HOXA10 by suppressing DRAM1 under IFNT treatment (Figure 6A). Furthermore, even in the absence of external intervention with IFNT, the mRNA expression of PTGS1 and PTGS2 were also significantly inhibited in siDRAM1 group compared with siNC group (Figures 6C,D). Nevertheless, interfering with DRAM1 exerted no distinct influence on the mRNA expression of HOXA10 and HOXA11 in BEND without IFNT treatment (Figures 6A,B).