*2.2. Ovary Histology*

The ovaries were fixed in 4% paraformaldehyde for 24 h before washing with 1 × PBS and dehydrated using a graded ethanol series (80% ethanol for 1 h, followed by 95% ethanol for 1 h, followed by 100% ethanol for 1 h). Transparency was improved using xylene (pure ethanol: xylene (1:1) for 1 h, then xylene for 1 h). Samples were infiltrated with paraffin (xylene: paraffin (1:1) at 62 ◦C for 1 h, then paraffin at 62 ◦C for 2 h) and processed for paraffin embedding. Sections were cut to 6 μm before staining with hematoxylin and eosin. Samples were scanned using a microscope slide scanner (Pannoramic MIDI, 3DHISTECH Ltd., Budapest, Hungary).

### *2.3. RNA Isolation, Library Construction and Illumine Sequencing*

TRIzol® reagent (Invitrogen, San Diego, CA, USA) extracted the total RNA by following the manufacturer's instructions. The total RNA was treated with DNase I to implement DNA digestion and obtain pure RNA products. Finally, the RNA purity was examined using NanoPhotometer® spectrophotometer NanoDrop 2000 (IMPLEN, Westlake Village, CA, USA). The RNA integrity and concentration was determined using Agilent 2100/4200 system software (Agilent Technologies, Santa Clara, CA, USA). Next, equal amounts of RNA from the same group of different individuals were pooled in order to construct the library. Next–generation sequencing library preparations were constructed according to the manufacturer's instructions (NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Ipswich, MA, USA)). After the mRNA library successfully passed the quality inspection, PE150 sequencing was performed using the Illumina NovaSeq 6000 platform (Thermo, Waltham, MA, USA).

#### *2.4. Basic Analysis of Sequencing Data*

We removed any technical sequences or their fragments (including adapters, polymerase chain reaction (PCR) primers), and any low quality sequences (with a base quality < 20). We achieved this by filtering the data (in the fastq format) using Trimmomatic (v0.30) to provide only high–quality, clean data in our analysis. Firstly, the whole genome sequence of *E. carinicauda* assembled by our research group was taken as the reference genome (the data have not been published yet). Next, Hisat2 (v2.0.1) was indexed to the reference genome sequence. Finally, we mapped the clean reads to the Silva database to remove the rRNA. All subsequent date analyses were based on clean data without rRNA.

#### *2.5. Differential Expression Genes (DEGs) Analysis and Enrichment Analysis*

The DESeq2 and edgeR [26] methods were used to perform the differential expression analysis. The influence of gene length and sequencing volume on the calculated gene expression level was eliminated by the fragments per kilobase per million reads (FPKM) method. The calculated gene expression level can be directly used to compare the expression variations between different genes. We used FPKM to normalize the data. The Benjamini and Hochberg's approach controlled the false discovery rate. |log2FC| ≥ 1 and *p*-values < 0.05 were set to detect any significant DEGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) enrichment analyses of differentially expressed gene sets were implemented by the topG and KOBAS packages 3.0 [27], respectively.

#### *2.6. RNA-Seq Data Validation by Real–Time Quantitative PCR*

Ten DEGs were randomly selected for quantitative reverse transcription–PCR (qRT– PCR) analysis to validate the differential expression of mRNAs. The qRT–PCR assay was performed using SYBR Green PCR Master Mix (Life Technologies, MASS, Waltham, MA, USA) in the 7500 fast Real–Time PCR system (Applied Biosystems, Foster, CA, USA) according to the manufacturer's agreement. The 18S rRNA of *E. carinicauda* was used as the internal reference [28]. All primer sequences and 18S rRNA sequences are listed in Table S1. The relative expression of target genes was calculated with 2−ΔΔCT methods. One–way ANOVA method and Duncan's test in the statistical software SPSS 22.0 (SPSS, Chicago, IL, USA) were used for statistical analysis. The results are presented as the mean ± standard error, and differences in gene expression were considered statistically significant at *p* < 0.05.
