Identification of photoperiod-induced specific miRNAs in the adrenal glands of Sunite sheep (Ovis aries)

In seasonal estrus, it is well known that melatonin-regulated biorhythm plays a key role. Some studies indicate that the adrenal gland plays an important role in reproduction in mammals, but the molecular mechanism is not clear. This study used an artificially controlled light photoperiod model, combined with RNA-seq technology and bioinformatics analysis, to analyze the messenger RNA (mRNA) and microRNA (miRNA) of ewe (Sunite) adrenal glands under different photoperiod treatments. After identification, the key candidate genes GRHL2, CENPF, FGF16 and SLC25A30 that photoperiod affects reproduction were confirmed. The miRNAs (oar-miR-544-3p, oar-miR-411b-5p, oar-miR-376e-3p, oar-miR-376d, oar-miR-376b-3p, oar-miR-376a-3p) were specifically expressed in the adrenal gland. The candidate mRNA-miRNA pairs (e.g., SLC25A30 coagulated by novel miRNA554, novel miRNA555 and novel miRNA559) may affect seasonal estrus. In summary, we constructed relation network of the mRNAs and miRNAs of sheep adrenal glands using RNA sequencing and bioinformatics analysis, thereby, providing a valuable genetic variation resource for sheep genome research, which will contribute to the study of complex traits in sheep.

Introduction

Low fecundity is the greatest limiting factor in the modern mutton sheep industry. Many element affect fecundity, such as litter size, oestrus frequency, and embryo survival rate. Because the most part area of China is located in the temperate zone, many animals are in oestrus in autumn and winter and give birth in spring and summer to ensure the survival of their offspring. Biologists are therefore divided into long-photoperiod breeds (LPs) and short-photoperiod breeds (SPs) with seasonal variations (1, 2). For example, the Chinese Sunite sheep is a typical SP breeder, which is specifically represented as estrus from August to March of the next year and anestrus from April to July (3, 4).

Concerning reproduction, the unavoidable theme is the upstream control center, such as the hypothalamus and pituitary gland. Studies have shown that light stimulates the paraventricular nucleus (PVN) and then affects the production of melatonin in the pineal gland, ultimately acting on the hypothalamus and affecting the reproduction of sheep through the hypothalamus-pituitary-gonad axis (HPG) (5–9). In addition, stimulating PVN also secretes corticotropin-releasing hormone (CRH), which then activates the pituitary gland to release corticotropin (ACTH). ACTH in turn stimulates the adrenal gland to secrete cortisol, which then provides negative feedback to the brain in the classic steady-state feedback loop to regulate hypothalamus-pituitary-adrenal gland axis (HPA) signals (10). However, the close relationship between HPA and HPG is more than that. When the fetus develops in the uterus, the hormone system that regulates the HPG and HPA axes plays an important role in the growth and development of its tissue. Fetal glucocorticoid concentrations increase in the third trimester of pregnancy, which is conducive to the modification of fetal key tissues or organs to promote fetal survival, including lung maturation and pulmonary surfactant production (11). In a study of rodents, exposure to glucocorticoids was found to closely affect fetal development during intrauterine development, such as gonadogenesis, the establishment of the HPG axis and the reproductive tract's morphogenesis (12, 13). In addition, the increase in the concentration of glucocorticoids interferes with the concentration of serum testosterone (T) because it inhibits T biosynthesis (14).

With the development of high-throughput sequencing, RNA sequencing (RNA-seq) is increasingly widely used in livestock. Including sheep (15) and cattle (16), this technique is used to obtain the expression profile information of mRNA and miRNA, which makes a great contribution to revealing some important traits and mining their candidate genes (17). Research on miRNAs has shown a trend of the great outbreak in recent years. Precursor miRNAs are well known to be transcribed mainly by RNA polymerase II, digested by dicer and then processed into mature miRNAs (18). Many studies have shown that miRNA plays a key role in regulating animal phenotype, such as affecting wool curvature (19), immunity to infectious diseases (20) and fat deposition (21). In this study, we mainly explored the key miRNAs affecting reproductive traits and adrenal-specific expression of miRNAs in sheep adrenal tissue under different photoperiod conditions. Finally, the relationship between differentially expressed mRNAs and miRNAs was predicted by bioinformatics software, and an interaction network is constructed, which is expected to mine effective information.

Materials and methods

Ethics statement

The experimental animals involved in this study were carried out after being examined by the Animal Experimental Welfare Ethics Committee of the Institute of Animal Sciences of Chinese Academy of Agricultural Sciences (IAS-CAAS, Beijing, China). In addition, the review acceptance number is No. IAS2018-3, and all the experimental procedures are executed in accordance with the relevant guidelines and regulations formulated by the Ministry of Agriculture and Rural Affairs of the People's Republic of China.

Preparation of animals

Nine non-pregnant adult Sunite ewes (aged 2–3 years old; weight 30–40 kg), which were randomly selected from a farm in Bayan Nur City (40°75′north latitude), Inner Mongolia Autonomous Region, China, were used for the study. The selected ewes were uniformly transported to a farm at the Tianjin Institute of Animal Sciences, Tianjin (39°13′ north latitude), China, and the follow-up experiments were carried out after a month of routine feeding to adapt to the local environment. The ovaries of these animals were removed by surgery, and an estrogen silicone tube (E₂, Sigma Chemical Co., St. Louis, MO, USA) was implanted subcutaneously in the neck of the sheep to maintain plasma estradiol levels of 3–5 pg/ml, as described previously (1, 22, 23). The ewe postoperative recovery lasted for 30 days before artificial light period control. The ewes were randomly divided into three groups: SP42 (short photoperiod for 42 days; n = 3), LP42 (long photoperiod for 42 days; n = 3) and SPLP42 (short photoperiod for 42 days followed by a long photoperiod for 42 days; n = 3). The conditions for the long photoperiod were 16 h of light per day and 8 h without light. The lighting duration setting for the short photoperiod was the opposite of the lighting duration setting for the long photoperiod. All sheep had ad libitum feeding and drinking in an enclosed climate control chamber with only artificial light sources.

Tissues acquisition, library construction and sequencing

Adrenal gland tissue from euthanized ewes was quickly preserved in liquid nitrogen. Then, the stored tissues were used for RNA extraction with TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions.

The mRNA library was constructed with 3 μg of high-quality RNA using the NEB Next Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, USA) according to the manufacturer's instructions. During this process, Ribo-Zero™ Gold Kits (TIANGEN, Beijing, China) were used to remove rRNA. In addition, we used the UNG enzyme to degrade the second strand of U-containing cDNA and performed polymerase chain reaction (PCR) amplification to obtain an RNA library. Then, a PE150 (paired-end 150 bp, PE150) sequencing approach for mRNAs was performed on a Hiseq X platform (Illumina, San Diego, CA, USA).

The miRNA library was built by a Small RNA Sample Pre kit (TIANGEN). We directly took total RNA as the starting sample, added adaptors at both ends of small RNA, and then reverse-transcribed the RNA to synthesize cDNA. After PCR amplification, the target DNA fragments were separated by polyacrylamide gel electrophoreses (PAGE), and the 140- to 160-bp ligation products were recovered to generate a cDNA library. In addition, an SE50 (single-end 50 bp, SE50) sequencing approach for miRNAs was performed on the Illumina HiSeq2500 platform (Illumina, San Diego, CA, USA). All sequencing data were outsourced to Annoroad Gene Technology Co., Ltd. (Beijing, China).

Data processing and transcriptome assembly

Bcl2fastq (v2.17.1.14) was used to process the offline data and convert the original image file into row sequencing reads on-base calling, which were row reads. Using an in-house Perl script made by Annoroad Genentech Co., Ltd. (Beijing, China) to remove low-quality reads, reads with adaptor contamination and reads with a rate of N > 5%, the clean mRNA reads were acquired from the raw reads. We used the Ovis aries reference genome (Oar_v4.0) and the genome annotation file from Ensembl. Cleaned reads were then mapped to the reference genome using HiSAT2 (v2.0.5) (24), and StringTie (v1.3.2d) was used to assemble the transcripts (25). HiSAT2 was run with “-rna - strandness RF” and “-dta -t -p 4,” String Tie with “-G ref.gtf -rf−1,” and the other parameters were set as the default. In addition to the above steps, the following steps were added to obtain the clean miRNA reads. These steps include removing reads without a 3'adaptor and insert fragment, removing the reads containing consecutive A/T/G/C bases, and removing the reads with abnormal final length. To ensure the accuracy of the subsequent analysis, the clean reads of sRNA sequencing were mapped to the reference genome (Oar_v4.0) by the comparison analysis software Bowtie v1.1.2.

Classification notes of sRNA and identification of miRNA

We obtained the situation, in which the sequence matched different regions in each sample by mapping the clean reads to the Ovis aries sequence in the miRBase database (Release 21) (26). At the same time, the known miRNA can be identified. The clean reads that were not annotated as a known miRNA were compared with the ncRNA sequence in Rfam (13.0) (27) to realize the annotation of rRNA, tRNA, snRNA, snoRNA and other ncRNA. The RepeatMasker program was used to comment on different types of repeats for clean reads that were not annotated as known miRNAs and ncRNAs (28). After identifying the above sRNA types and then using the matching (100% positional overlap) results with the location information of exons and introns of the gene, the sRNA from mRNA will be annotated (29). For sRNA reads that did not match the above-known annotation type, the software miRDeep2 (30) was used to predict the new miRNA, and the sequence, expression and structure information of each new miRNA were analyzed. Different mature body sequences, precursor sequences and positions will be considered different new miRNAs.

Differential expression and functional enrichment analysis of miRNA and mRNAs

The transcripts per million (TPM) and the fragments per kilobase per million mapped reads (FPKM) values were calculated to represent miRNA and mRNA expression levels based on the read number (31). The difference analysis was carried out by the software DESeq (32). The screening conditions of differentially expressed miRNAs were pval ≤ 0.05, padj ≤ 0.05 and log₂| (fold change) | ≥1. The differential expression criteria of mRNA were |log₂Ratio| ≥1 and q < 0.05. In addition, the target genes of known and novel differentially expressed miRNAs were predicted by miRanda and TargetScan software, respectively (33). The intersection of the prediction results of the two software programs was selected as the target gene of the miRNAs.

We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on the targeted genes of differentially expressed miRNAs (DEMs) and differentially expressed mRNAs (DEGs). The hypergeometric test method was applied to assess significantly enriched GO terms and KEGG pathways, and FDR <0.05 was considered to be significantly enriched.

Construction of integral miRNA–mRNA interaction networks

To further describe the interaction between miRNA and mRNA, we used miRanda and TargetScan software for prediction and took the intersection of the prediction results of the two software programs as the target genes of miRNA. According to the identified DEMs and differentially expressed target mRNAs of DEMs, a regulatory network diagram was drawn by Cytoscape software (34).

Data validation

Transcripts (n = 4) were randomly selected and the primers were designed by Primer 5.0 software. The miRNA primers were synthesized by Shanghai Sangon Biotech. The miRNA quantitative (q)PCR conditions were as follows: 95°C for 15 min, followed by 40 cycles of 94°C for 20 s and 60°C for 34 s (miRcute Plus miRNA qPCR Kit, TIANGEN). In addition, U6 small nuclear RNA (snRNA; for miRNA) were utilized as reference genes. The data obtained from the qPCR were evaluated using the 2^−ΔΔCt method.

Results

Overview of sequencing data in adrenal gland tissue

We obtained the global expression profile of mRNA and miRNA in the adrenal glands of Sunite ewes under different photoperiod treatments. The data of subsequent analysis, namely, clean reads, were based on the filtered original offline data, which had no adapter-polluted reads, and the base error rate was <0.1%. Overall, we obtained 1,029 million clean reads of mRNA and 244 million clean reads of miRNA. To obtain more accurate sequences and the accuracy of subsequent analysis, bioinformatics analysis software was used to map clean reads to the reference genome and compare the results statistically. The results are shown in Table 1. We found that the average comparison rate of miRNA was ~50%, while the average comparison rate of mRNA was approximately 95%. In addition, 123 million perfect match reads of miRNA and 978 million perfect match reads of mRNA were obtained separately (Table 1; Supplementary Table 1).

TABLE 1

Table 1. Summary of the mapping data from the adrenal gland tissues.

Identification of mRNAs and miRNAs in adrenal gland tissue

A total of 18,947 mRNAs were identified (Supplementary Table 2) after mapping to the sheep genome. According to the Ovis aries gene annotation files in the related database, we counted the number and proportion of the unique alignment sequences on the three functional elements of the gene (exon, intron and intergenic). In general, if the annotation information of the species is more comprehensive, most of the sequences should be aligned to the exon region, but alternative splicing and noise expression will cause some sequences to come from the intron region and the intergenic region. The following figure shows the distribution of unique alignment sequences in various regions of the genome and chromosome. Our results suggested that many mRNAs were situated in the intergenic region (more than 45%), followed by the exon (approximately 33%) and intron (nearly 20%) regions (Figure 1A; Supplementary Table 3). Furthermore, the chromosome distribution of mRNAs showed that chromosome 3 and chromosome 1 expressed the most genes, accounting for 9.62 and 9.65%, respectively, followed by chromosome 2 (7.47%) and chromosome 5 (5.67%) (Figure 1B; Supplementary Table 4).

FIGURE 1

Figure 1. Summary information of mRNA identification. (A) According to the sequence alignment information, the number of sequences aligned to exons, introns and intergenic regions was calculated, and the pie chart was made according to the proportion. (B) Chromosome distribution of identified genes from the adrenal gland.

Regarding miRNAs, RNA-seq generated approximately 244 million clean reads after filtering and mapping. We classified and annotated all clean reads, including known miRNAs, ncRNAs (rRNA\ tRNA\ snRNA\ snoRNA\ other rfam), repeats, and small RNAs annotated by perfect matching with mRNA exons and introns and novel miRNAs. The detailed classification and specific statistics are shown in Table 2. The proportion of clean reads matching known miRNA maturity in each sample was more than 70%, followed by the proportion of clean reads matching with other ncRNA in the rfam database being more than 15%. In addition, we also counted the number of clean read types that matched the miRNA maturity, namely, unique clean reads. The statistical results are shown in Figure 2A. We found that the following four types, namely, exon+, other, rRNA and intron+, accounted for more than 10% of unique clean reads in the three experimental groups. Compared with the statistical results of total reads, the proportion of unique clean reads with known miRNA typeS decreased significantly, from ~70 to 1% (Supplementary Table 6). The distribution of chromosomes is similar to the distribution of chromosomes of mRNA. The results of Figure 2B show that the X chromosome contains the most miRNA (12.25%), followed by chromosomes 3 (9.86%), 2 (8.73%) and 1 (7.89%) (Supplementary Table 4).

TABLE 2

Table 2. Summary of the classification statistics of sRNA total reads mapped to reference genome.

FIGURE 2

Figure 2. Summary information of sRNA identification. (A) The statistical results of unique clean reads from the adrenal gland in the three experimental groups. (B) Chromosome distribution of identified miRNA from the adrenal gland.

Furthermore, our main concern is the identification of known miRNAs and novel miRNAs. In total, we identified miRNAs (n = 861), including known miRNAs (n = 151) and novel miRNAs (n = 710). The numbers of total clean reads identified as known miRNAs and novel miRNAs in the three experimental group samples were 93,377,486 and 177,217, respectively, and the numbers of unique clean reads identified as known miRNAs and novel miRNAs in the samples were 18,691 and 1,054, respectively. Figure 3 shows an example of miRDeep2 output. The top of Figure 3 shows the scores assigned to each part of the miRNA and the total count. Different colors are used to represent different parts of the predicted hairpin secondary structure. All readings related to miRNA are also shown at the bottom of the figure. Therefore, about the detailed miRDeep2 output information of novel miRNAs, we described in detail the three biological repetitive samples of the three processing groups (LP42-1 n = 115, LP42-2 n = 120, LP42-3 n = 119, SP42-1 n = 94, SP42-2 n = 138, SP42-3 n = 110, SPLP42-1 n = 117, SPLP42-2 n = 120, SPLP42-3 n = 121), including the precursor secondary structure, count number and read sequence shown in the picture. Known miRNA is also described in detail (LP42-1 n = 106, LP42-2 n = 106, LP42-3 n = 106, SP42-1 n = 104, SP42-2 n = 106, SP42-3 n = 105, SPLP42-1 n = 106, SPLP42-2 n = 106, SPLP42-3 n = 106). The known miRNA specifically expressed in the adrenal gland was found by comparing with Zhang's Small Tail Han sheep hypothalamic miRNA transcriptome data (17), only six miRNAs (oar-miR-544-3p, oar-miR-411b-5p, oar-miR-376e-3p, oar-miR-376d, oar-miR-376b-3p, oar-miR-376a-3p) were specifically expressed in the adrenal gland, and one (oar-miR-1193-3p) was specifically expressed in the hypothalamus (Supplementary Tables 5, 6, miRDeep2 file).

FIGURE 3

Figure 3. An example of the miRDeep2 output. The figure illustrates the output for the miRNA of novel_555. The upper part of Figure 3 shows the read count and the total count of each part. The predicted secondary structure of the hairpin is also depicted, with mature (red), star (purple) and circular (yellow) sequences highlighted in different colors. The bottom of the map shows the results of miRNA alignment with predicted precursor sequences on the genome (OBS line) and experimental sequences reported in miRBase (EXP line). For each sequence, the frequency and mismatch with the genome sequence (mm sequence) are given. Mismatches are also highlighted in uppercase letters.

The analysis of differentially expressed miRNAs and mRNAs in adrenal gland tissue

The numbers of DEGs and DEMs identified from the LP42 vs. SPLP42 comparison group were 144 and 48, respectively. Among these DEGs and DEMs, 45 and 10 genes were upregulated, and 99 and 38 genes were downregulated. In the same way, we counted the DEGs and DEMs of the SP42 vs. LP42 (DEG n = 454, DEM n = 36) and SP42 vs. SPLP42 (DEG n = 506, DEM n = 55) comparison groups. After analyzing the differences between all known miRNAs and novel miRNAs, we found that in known miRNAs, only the oar-miR-148a of the LP42 vs. SPLP42 comparison group was differentially expressed, so the graph shows only the statistics of novel miRNAs. Overall details are shown in Figure 4; Supplementary Table 7.

FIGURE 4

Figure 4. Differentially expressed (DEM) microRNA and (DEG) mRNA analysis. Different color rectangles were used to represent the numbers of upregulated (blue), downregulated (yellow) and total (green) differential genes.

GO and KEGG enrichment analysis of DEGs and target gene of DEMs

To better understand the potential functions of the DEGs and target genes of DEMs, GO term and KEGG pathway analyses were performed. For DEGs, the three comparison groups we set up showed large differences. We found that the number of GO terms for biological processes significantly enriched in the three comparison groups was more than 50; in particular, the term number in the SP42 vs. SPLP42 comparison group was as high as 665. We selected only the top five terms in each group to draw the resulting chart. Second, we found a total of 10 terms in the three comparison groups regarding cellular components, so they were all drawn on the resulting graph. For the part of molecular function, except that we only found 4 significant terms in LP42 vs. SPLP42, the other two comparison groups chose the top five terms to draw the resulting map (Supplementary Table 9). Our main focus on BP terms included regulation of multicellular organismal process and regulation of the developmental process. The NAADP-sensitive calcium-release channel activity of the MF term also attracted our attention. The results of the KEGG pathway analysis are shown in Figure 5, and all the significantly enriched pathways found by the three comparison groups are plotted on the resulting map (Supplementary Table 10). The main pathways that attracted our attention were cytokine-cytokine receptor interactions, the TNF signaling pathway, the Jak-STAT signaling pathway and the MAPK signaling pathway.

FIGURE 5

Figure 5. Functional enrichment analysis of DEGs. The upper part of the picture shows the enrichment results of GO and the lower part shows the enrichment results of KEGG.

Regarding the target genes of miRNAs, LP42 vs. SPLP42 (n = 119), SP42 vs. LP42 (n = 245) and SP42 vs. SPLP42 (n = 224) differential target genes were found (Supplementary Table 8). The functional enrichment analysis of miRNA target genes was divided mainly into known miRNAs and novel miRNAs. The results are shown in Table 3. In the taxonomic list of known miRNAs, we found only that there was significant enrichment of GO terms in the LP42 vs. SPLP42 comparison group, except that the first five terms were selected in the BP, and all the significantly enriched terms about CC and MF are shown in the rest of the list. We are most interested in Wnt-activated receptor activity, Wnt-protein binding and cytokine binding. However, the results of GO enrichment analysis of novel miRNA target genes showed that only 13 terms in the MF part of the SP42 vs. SPLP42 comparison group were significantly enriched (Supplementary Table 11). Finally, we were also surprised to find that there was no significant enrichment of pathways in either known miRNAs classification or novel miRNAs classification in the three comparison groups in the KEGG enrichment analysis. However, there were some pathways with p ≤ 0.05, not Q ≤ 0.05 (Supplementary Table 12).

TABLE 3

Table 3. Functional enrichment analysis of target gene of DEMs.

Analysis of integrated miRNA–mRNA co-expression network

To fully understand the potential roles of miRNAs, we built interactome networks using DEMs and their targets (DEGs). In total, 16 DE miRNAs (novel miRNAs) in SP42 vs. LP42 were predicted to target 67 genes. An mRNA–miRNA co-expression network was then constructed, where 1 DEG was targeted by 5 novel miRNAs. Regarding SP42 vs. SPLP42, 35 novel DE miRNAs were predicted to target 61 genes. The other mRNA–miRNA co-expression network was then constructed, where 4 DEGs were targeted by 16 novel miRNAs (Figure 6; Supplementary Tables 8, 13).

FIGURE 6

Figure 6. Overview of mRNA–miRNA networks. Different colors represent DEM (green) and DEG (yellow), and different shapes represent different modes of regulation.

Data validation

To evaluate the accuracy of sequencing, qPCR was used to verify the RNA-seq data. The results showed that the expression pattern of miRNAs in sheep adrenal glands was similar to the expression pattern of sequencing (Figure 7), which proved the reliability of the data produced by RNA-seq.

FIGURE 7

Figure 7. Quantitative verification results of miRNA. The relative expressions of 4 miRNA randomly selected in different individuals and among comparison groups were verified. Rectangles with different patterns represent RNA-seq (A) and qPCR (B) results, respectively.

Discussion

The transcriptome is space-time specific, but it is a powerful means to mine key information in the process of animal genetic development. The purpose of this study was to establish an OVX model and assist the light-controlled environment to identify the key sheep miRNAs and differentially expressed genes affected by light changes. Our previous studies successfully obtained the miRNAome and mRNAome information of the hypothalamus and pituitary through a similar experiment (17, 22). Therefore, we selected an adrenal transcriptional group for detailed analysis based on the close relationship between HPG and the HPA axis and successfully identified new adrenal-specific miRNAs and differentially expressed genes.

In this study, a total of 861 miRNA and 17,712 coding protein genes were identified in the three comparison groups. Among the miRNAs we identified, novel miRNAs accounted for approximately 82.45% of the total number of miRNAs, and known miRNAs accounted for only 17.54%. Among the 151 mature known miRNAs, we obtained only 106 stem-loop sequences. Both the 5 p arm and 3 p arm of 47 stem-loop sequences produced mature miRNA. Both arms of a miRNA precursor may give rise to functional levels of mature miRNA (35, 36). The dominant products may change from species to species and have different tissue expression preferences, including normal vs. pathological tissue (37–43). At the same time, when identifying three biological repetitive samples in each treatment group, we found that except for SP42-1 and SP42-3, all the other individuals identified 106 stem-loop sequences. The stem-loop sequences that were not detected were oar-mir-654 (producing 3P mature miRNA) and oar-mir-1193 (producing 5p mature miRNA). mir-1193 has been found to inhibit the proliferation and invasion of cancer cells by directly acting on transmembrane 9 superfamily member 3 (TM9SF3) and insulin-like growth factor 2 mRNA binding protein 2 (IGF2BP2) (44, 45). In a mouse study, overexpression of mir-1193-5p was found to be able to increase the differentiation tendency of oligodendrocyte progenitor cells (46). Similarly, previous studies have found that mir-654 is significantly related to the development of a variety of malignant tumors, including lung cancer, rectal cancer, esophageal cancer and colon cancer (47, 48). Since we only detected the above two miRNAs in SP42-2 individuals, it is not possible to determine their relationship with the photoperiod and adrenal gland, so we will explore further.

After differential expression analysis, we initially identified 144, 454 and 506 DEGs and 48, 36 and 55 DE miRNAs (LP42 vs. SPLP42, SP42 vs. LP42 and SP42 vs. SPLP42). According to the grouping, the quantity distribution of DEGs was more consistent with the expectation of the experiment; that is, the number of differentially expressed genes in the SP42 vs. LP42 and SP42 vs. SPLP42 comparison groups was greater than the number of differentially expressed genes in the LP42 vs. SPLP42 comparison group. However, the results of miRNA analysis showed that there were the fewest DEMs in the SP42 vs. LP42 comparison group. When studying the data of the sheep pituitary transcriptome after different photoperiod treatments, with the increase in long photoperiod maintenance days, the number of differentially expressed genes and differentially expressed lncRNAs were found to increase significantly compared with the short photoperiod (22). What causes the number of differential miRNAs to appear in the small probability of the significant difference treatment group needs further analysis. Similarly, during a miRNA study of sheep hypothalamus, significant expression of the let-7 and oar-miRNA-200 families was found (17), and some studies also proved that those identified miRNAs were differentially expressed in seasonal and non-seasonal sheep breeds (49). However, it is unfortunate that although the let-7 and oar-miRNA-200 families were detected in this study, there was no differential expression in different comparison groups. In our known miRNA identification of sheep adrenal glands and Small Tail Han sheep hypothalamus (17), only six miRNAs (oar-miR-544-3p, oar-miR-411b-5p, oar-miR-376e-3p, oar-miR-376d, oar-miR-376b-3p, oar-miR-376a-3p) were specifically expressed in the adrenal gland, and one (oar-miR-1193-3p) was specifically expressed in the hypothalamus. The rest can be detected in both tissues. In addition, several miRNAs, such as miRNA-200 family members, were thought to be conserved in the hypothalamus of mice (16, 50), rats (51), zebrafish (52), and sheep (17). However, to sum up, it may be conserved in the species, not in sheep breeds and the tissue; perhaps it is just the wrong time to monitor.

All the DEGs were subjected to GO functional enrichment and KEGG enrichment analysis according to the set comparison group. Because there were too many enrichment items in the BP part and the MF part, we chose the top five terms that were significantly enriched. In the BP section, we find that there are coincident entries GO:0051239 and GO:0050793 in the SP42 vs. LP42 and SP42 vs. SPLP42 comparison groups. The overlapping genes GRHL2 (grainy head-like protein 2 homologs), CENPF (centromere protein F) and FGF16 (fibroblast growth factor 16) were found by further analysis of the above two terms. GRHL2 is one of three mouse homologs of Drosophila Grainyhead (53) and is expressed in diverse embryonic epithelial tissues during development (53, 54). GRHL2 and its paralogue Grhl3 play essential roles in neural tube closure in mice (55–58). Data indicate that a conserved GRHL2-coordinated gene network controls trophoblast branching morphogenesis, thereby facilitating the development of the site of feto-maternal exchange (59). Normal placental development is a key factor in reproductive success, so whether screening this gene in our different photoperiod comparison groups indicates that placental development is affected by photoperiod is unknown. Therefore, it is very important to explore its function in the process of seasonal reproduction of sheep. Another key gene, CENPF, related to early embryonic development was also screened out. CENPF is a member of the centromere protein family that regulates chromosome segregation during mitosis (60, 61). Zhou et al. used mice to study and indicated that farnesylation plays a key role during CENPF degradation and localization in early embryos. At the same time, a knockout test of the CENPF gene showed that the gene is very important for early embryos (62). The fibroblast growth factor 16 (FGF16) gene is correlated with oocyte quality (63). In summer, when oocyte quality is low, the expression of FGF16 is low. Conversely, in winter, when oocyte quality is high, the expression of FGF16 is high. We found that the expression of the FGF16 gene was upregulated in the comparative group of SP42 vs. LP42, indicating that in adrenal tissue, the expression of the FGF16 gene under a long photoperiod was nearly 10 times lower than that under short photoperiod.

Our experimental analysis results are not ideal, whether it is the distribution of differential miRNAs in different comparison groups or the functional enrichment analysis of GO and KEGG of DEGs or the target genes of DEMs. However, the key comparison group SP42 vs. LP42 we set up has only made exciting discoveries in the process of GO enrichment analysis of mRNA. All the significantly enriched KEGG items of DEGs or target genes of DEMs are shown in the results section, but the SP42 vs. LP42 comparison group unexpectedly did not have any significant entries. Therefore, we need to further construct a miRNA–mRNA coexpression network to study interaction regulation and mine gene expression regulation patterns. We identified the overlapping core gene SLC25A30 in the SP42 vs. LP42 and SP42 vs. SPLP42 comparison groups. SLC25A30 is a member of the mitochondrial transporter family (64). Mitochondrial transporters for inorganic anions/malate (SLC25A30), thiamine pyrophosphate (TPP) (SLC25A19) and iron (SLC25A28) are also considered conceptus-induced IFNT-dependent endometrial DEGs and are increased in endometrial epithelial cells by IFNT (65, 66). A study by Gorgoglione et al. (67) suggests that SLC25A30 may have the function of exporting sulfite and thiosulfate or transporting malate acid, which may contribute to the flow of malate–aspartate acid in mitochondria. Importantly, malic acid and malate dehydrogenase is present in the intrauterine environment of cattle, and malic acid has been shown to affect the early embryonic development of hamsters (68–70). SLC25A30 gene expression has been reported to be upregulated by cellular oxidative stress, and our study may also be a sign of increased oxidative stress (71). Novel miRNA554, novel miRNA555 and novel miRNA559 are the common miRNAs that regulate the core gene in the two comparison groups (SP42 vs. LP42 and SP42 vs. SPLP42). These three miRNAs are located on chromosome 7 and belong to the same miRNA cluster according to location analysis. However, the specific functions need to be further analyzed.

These results suggest that several key DEGs and DEMs in the adrenal gland are directly or indirectly involved in the process of photoperiod-changing reproductive activity, and further study of gene/miRNA knockout or overexpression is helpful for us to understand their real function in female reproductive traits.

Conclusions

We successfully obtained the mRNA and small RNA data of adrenal gland tissue samples of Sunite sheep under different photoperiod treatments by RNA sequencing and bioinformatics analysis. The key candidate genes of photoperiod affecting reproduction, such as GRHL2, CENPF, FGF16 and SLC25A30, were confirmed. The miRNA (oar-miR-544-3p, oar-miR-411b-5p, oar-miR-376e-3p, oar-miR-376d, oar-miR-376b-3p, oar-miR-376a-3p) was specifically expressed in adrenal tissue. The predicted mRNA-miRNA pairs (SLC25A30 regulated by Novel miRNA554, Novel miRNA555 and Novel miRNA559) showed significant differences in SP42 vs. the other two groups, indicating that the relationship may play an important role in the process of different photoperiod affecting adrenal function. Our results provide a new perspective for the study of sheep reproduction and help to deepen the understanding of ovine reproduction.

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 below: https://www.ncbi.nlm.nih.gov/sra/PRJNA756142; https://www.ncbi.nlm.nih.gov/sra/PRJNA811389.

Ethics statement

The animal study was reviewed and approved by the Animal Experimental Welfare Ethics Committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences.

Author contributions

QiuL and MC designed the research. XD wrote the paper and performed the study. RD, XD, and QinL collected the data. XD and XH analyzed data. MC revised the final manuscript. All authors reviewed the manuscript and approved the final version.

Funding

This work was financially supported by National Natural Science Foundation of China (32172704), China Agriculture Research System of MOF and MARA (CARS-38), and the Agricultural Science and Technology Innovation Program of China (CAAS-ZDRW202106 and ASTIP-IAS13).

Acknowledgments

We are grateful to Xiaosheng Zhang, Jinlong Zhang, and Xiaofei Guo in Tianjin Institute of Animal Sciences for the ovariectomy, estrogen implantation, and ewes feeding. Thanks to Xiaohan Cao, Xiaoyu Li, Xinlong Dong, and Qing Xia for their help in the sample collection.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2022.888207/full#supplementary-material

Abbreviations

mRNA, message RNA; miRNA, microRNA; HPA, The hypothalamic–pituitary–adrenal; LP, Long photoperiod; SP, Short photoperiod; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LP42, The LP lasts for 42 days; SP-LP42, The SP lasts for 42 days, and then convert it to LP for 42 days; SP42, The SP lasts for 42 days; GRHL2, grainy head-like protein 2 homologs; CENPF, centromere protein F; FGF16, Fibroblast growth factor 16; SLC25A30, solute carrier family 25, member 30.

References

1. La YF, He XY, Zhang LP, Di R, Wang XY, Gan SQ, et al. Comprehensive analysis of differentially expressed profiles of mRNA, lncRNA, and circRNA in the uterus of seasonal reproduction sheep. Genes-Basel. (2020) 11:301–19. doi: 10.3390/genes11030301

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ebling FJP FD. Photoperiod requirements for puberty differ from those for the onset of the adult breeding season in female sheep. J Reprod Fertil. (1988) 84:283–93. doi: 10.1530/jrf.0.0840283

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Dupre SM, Miedzinska K, Duval CV, Le Y, Goodman RL, Lincoln GA, et al. Identification of Eya3 and TAC1 as long-day signals in the sheep pituitary. Curr Biol. (2010) 20:829–35. doi: 10.1016/j.cub.2010.02.066

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Li XY, He XY, Liu QY, Wang XY, Guo XF, Xia Q, et al. Expression pattern analysis of TAC1 and PRLR genes in different reproductive states of sheep. Acta Vet Zootech Sin. (2018) 49:253–62. doi: 10.11843/j.issn.0366-6964.2018.02.004

CrossRef Full Text | Google Scholar

5. Batailler M, Chesneau D, Derouet L, Butruille L, Segura S, Cognie J, et al. Pineal-dependent increase of hypothalamic neurogenesis contributes to the timing of seasonal reproduction in sheep. Sci Rep. (2018) 8:6188. doi: 10.1038/s41598-018-24381-4

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Karsch FJ, Bittman EL, Foster DL, Goodman RL, Legan SJ, Robinson JE. Neuroendocrine basis of seasonal reproduction. Recent Prog Horm Res. (1984) 40:185–232. doi: 10.1016/B978-0-12-571140-1.50010-4

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Weems PW, Goodman RL, Lehman MN. Neural mechanisms controlling seasonal reproduction: principles derived from the sheep model and its comparison with hamsters. Front Neuroendocrinol. (2015) 37:43–51. doi: 10.1016/j.yfrne.2014.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Masumoto K, Ukai-Tadenuma M, Kasukawa T, Nagano M, Uno KD, Tsujino K, et al. Acute induction of Eya3 by late-night light stimulation triggers TSH beta expression in photoperiodism. Neurosci Res. (2011) 71:E172. doi: 10.1016/j.neures.2011.07.742

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hanon EA, Lincoln GA, Fustin JM, Dardente H, Masson-Pevet M, Morgan PJ, et al. Ancestral TSH mechanism signals summer in a photoperiodic mammal. Curr Biol. (2008) 18:1147–52. doi: 10.1016/j.cub.2008.06.076

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Acevedo-Rodriguez A, Kauffman AS, Cherrington BD, Borges CS, Roepke TA, Laconi M. Emerging insights into hypothalamic-pituitary-gonadal axis regulation and interaction with stress signalling. J Neuroendocrinol. (2018) 30:e12590. doi: 10.1111/jne.12590

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Seckl JR, Holmes MC. Mechanisms of Disease: glucocorticoids, their placental metabolism and fetal'programming' of adult pathophysiology. Nat Clin Pract Endoc. (2007) 3:479–88. doi: 10.1038/ncpendmet0515

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Pereira OCM, Arena AC, Yasuhara F, Kempinas WG. Effects of prenatal hydrocortisone acetate exposure on fertility and sexual behavior in male rats. Regul Toxicol Pharm. (2003) 38:36–42. doi: 10.1016/S0273-2300(03)00046-1

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Piffer RC, Garcia PC, Pereira OCM. Adult partner preference and sexual behavior of male rats exposed prenatally to betamethasone. Physiol Behav. (2009) 98:163–7. doi: 10.1016/j.physbeh.2009.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Hardy MP, Gao HB, Dong Q, Ge RS, Wang Q, Chai WR, et al. Stress hormone and male reproductive function. Cell Tissue Res. (2005) 322:147–53. doi: 10.1007/s00441-005-0006-2

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Jiang Y, Xie M, Chen W, Talbot R, Maddox JF, Faraut T, et al. The sheep genome illuminates biology of the rumen and lipid metabolism. Science. (2014) 344:1168–73. doi: 10.1126/science.1252806

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Correia CN, McLoughlin KE, Nalpas NC, Magee DA, Browne JA, Rue-Albrecht K, et al. RNA Sequencing (RNA-seq) reveals extremely low levels of reticulocyte-derived globin gene transcripts in peripheral blood from horses (Equus caballus) and cattle (Bos taurus). Front Genet. (2018) 9:278. doi: 10.3389/fgene.2018.00278

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Zhang ZB, Tang JS, Di R, Liu QY, Wang XY, Gan SQ, et al. Integrated hypothalamic transcriptome profiling reveals the reproductive roles of mRNAs and miRNAs in Sheep. Front Genet. (2019) 10:1296. doi: 10.3389/fgene.2019.01296

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Gebert LFR, MacRae IJ. Regulation of microRNA function in animals. Nat Rev Mol Cell Biol. (2019) 20:21–37. doi: 10.1038/s41580-018-0045-7

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Lv XY, Chen WH, Wang SH, Cao XK, Yuan ZH, Getachew TF, et al. Integrated hair follicle profiles of microRNAs and mRNAs to reveal the pattern formation of Hu sheep lambskin. Genes. (2022) 13:342. doi: 10.3390/genes13020342

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Chen S, Wang CQ, Chen QL, Zhao DT, Liu YB, Du L, et al. Downregulation of three novel miRNAs in the lymph nodes of sheep immunized with the brucella suis strain 2 vaccine. Front Vet Sci. (2022) 9:813170. doi: 10.3389/fvets.2022.813170

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Fei XJ, Jin ML, Wang YQ, Li TT, Lu ZK, Yuan ZH, et al. Transcriptome reveals key microRNAs involved in fat deposition between different tail sheep breeds. PLoS ONE. (2022) 17:e0264804. doi: 10.1371/journal.pone.0264804

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Xia Q, Chu MX, He XY, Zhang XS, Zhang JL, Guo XF, et al. Identification of photoperiod-induced lncRNAs and mRNAs in pituitary pars tuberalis of sheep. Front Vet Sci. (2021) 8:644474. doi: 10.3389/fvets.2021.753614

PubMed Abstract | CrossRef Full Text | Google Scholar

23. He XY, Tao L, Zhong YJ, Di R, Xia Q, Wang XY, et al. Chu MX: Photoperiod induced the pituitary differential regulation of lncRNAs and mRNAs related to reproduction in sheep. PeerJ. (2021) 9:e10953. doi: 10.7717/peerj.10953

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. (2016) 11:1650–67. doi: 10.1038/nprot.2016.095

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. (2015) 33:290. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Griffiths-Jones S. miRBase: the microRNA sequence database. Methods Mol Biol. (2006) 342:129–38. doi: 10.1385/1-59745-123-1:129

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Burge SW, Daub J, Eberhardt R, Tate J, Barquist L, Nawrocki EP, et al. Rfam 11.0: 10 years of RNA families. Nucleic Acids Res. (2013) 41:D226. doi: 10.1093/nar/gks1005

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Huda A, Jordan IK. Analysis of transposable element sequences using CENSOR and RepeatMasker. Methods Mol Biol. (2009) 537:323–36. doi: 10.1007/978-1-59745-251-9_16

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. (2010) 26:841–2. doi: 10.1093/bioinformatics/btq033

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Friedlander MR, Mackowiak SD Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. (2012) 40:37–52. doi: 10.1093/nar/gkr688

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. (2010) 28:511–U174. doi: 10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. (2010) 11:R106. doi: 10.1186/gb-2010-11-10-r106

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. (2003) 5:R1. doi: 10.1186/gb-2003-5-1-r1

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Du T, Zamore PD. MicroPrimer: the biogenesis and function of microRNA. Development. (2005) 132:4645–52. doi: 10.1242/dev.02070

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Rana TM. Illuminating the silence: understanding the structure and function of small RNAs. Nat Rev Mol Cell Biol. (2007) 8:23–36. doi: 10.1038/nrm2085

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Griffiths-Jones S, Hui JH, Marco A, Ronshaugen M. MicroRNA evolution by arm switching. EMBO Rep. (2011) 12:172–7. doi: 10.1038/embor.2010.191

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ro S, Park C, Young D, Sanders KM, Yan W. Tissue-dependent paired expression of miRNAs. Nucleic Acids Res. (2007) 35:5944–53. doi: 10.1093/nar/gkm641

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Ruby JG, Stark A, Johnston WK, Kellis M, Bartel DP, Lai EC. Evolution, biogenesis, expression, and target predictions of a substantially expanded set of Drosophila microRNAs. Genome Res. (2007) 17:1850–64. doi: 10.1101/gr.6597907

PubMed Abstract | CrossRef Full Text | Google Scholar

40. De Wit E, Linsen SE, Cuppen E, Berezikov E. Repertoire and evolution of miRNA genes in four divergent nematode species. Genome Res. (2009) 19:2064–74. doi: 10.1101/gr.093781.109

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Chiang HR, Schoenfeld LW, Ruby JG, Auyeung VC, Spies N, Baek D, et al. Mammalian microRNAs: experimental evaluation of novel and previously annotated genes. Genes Dev. (2010) 24:992–1009. doi: 10.1101/gad.1884710

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Marco A, Macpherson JI, Ronshaugen M, Griffiths-Jones S. MicroRNAs from the same precursor have different targeting properties. Silence. (2012) 3:8. doi: 10.1186/1758-907X-3-8

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Li SC, Tsai KW, Pan HW, Jeng YM, Ho MR, Li WH. MicroRNA 3' end nucleotide modification patterns and arm selection preference in liver tissues. BMC Syst Biol. (2012) 6:S14. doi: 10.1186/1752-0509-6-S2-S14

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Li XL, Li YH, Lu H. MiR-1193 suppresses proliferation and invasion of human breast cancer cells through directly targeting IGF2BP2. Oncol Res. (2017) 25:579–85. doi: 10.3727/97818823455816X14760504645779

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Shen LY, Du XJ, Ma HY, Mei SX. MiR-1193 suppresses the oroliferation and invasion of human T-cell leukemia cells through directly targeting the transmembrane 9 superfamily 3 (TM9SF3). Oncol Res. (2017) 25:1643–51. doi: 10.3727/096504017X14908284471361

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Ma Q, Matsunaga A, Ho B, Oksenberg JR, Didonna A. Oligodendrocyte-specific Argonaute profiling identifies microRNAs associated with experimental autoimmune encephalomyelitis. J Neuroinflamm. (2020) 17:297. doi: 10.1186/s12974-020-01964-5

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Pu JT, Hu Z, Zhang DG, Zhang T, He KM, Dai TY. MiR-654-3p suppresses non-small cell lung cancer tumourigenesis by inhibiting PLK4. Oncotargets Ther. (2020) 13:7997–8008. doi: 10.2147/OTT.S258616

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Li P, Cai JX, Han F, Wang J, Zhou JJ, Shen KW, et al. Expression and significance of miR-654-5p and miR-376b-3p in patients with colon cancer. World J Gastro Oncol. (2020) 12:492–502. doi: 10.4251/wjgo.v12.i4.492

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Zhai M, Xie Y, Liang H, Lei X, Zhao Z. Comparative profiling of differentially expressed microRNAs in estrous ovaries of Kazakh sheep in different seasons. Gene. (2018) 664:181–91. doi: 10.1016/j.gene.2018.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Choi PS, Zakhary L, Choi WY, Caron S, Alvarez-Saavedra E, Miska EA, et al. Members of the miRNA-200 family regulate olfactory neurogenesis. Neuron. (2008) 57:41–55. doi: 10.1016/j.neuron.2007.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Sangiao-Alvarellos S, Pena-Bello L, Manfredi-Lozano M, Tena-Sempere M, Cordido F. Perturbation of hypothalamic microRNA expression patterns in male rats after metabolic distress: impact of obesity and conditions of negative energy balance. Endocrinology. (2014) 155:1838–50. doi: 10.1210/en.2013-1770

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Garaffo G, Conte D, Provero P, Tomaiuolo D, Luo Z, Pinciroli P, et al. The Dlx5 and Foxg1 transcription factors, linked via miRNA-9 and−200, are required for the development of the olfactory and GnRH system. Mol Cell Neurosci. (2015) 68:103–19. doi: 10.1016/j.mcn.2015.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Wilanowski T, Tuckfield A, Cerruti L, O'Connell S, Saint R, Parekh V, et al. A highly conserved novel family of mammalian developmental transcription factors related to Drosophila grainyhead. Mech Dev. (2002) 114:37–50. doi: 10.1016/S0925-4773(02)00046-1

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Auden A, Caddy J, Wilanowski T, Ting SB, Cunningham JM, Jane SM. Spatial and temporal expression of the Grainyhead-like transcription factor family during murine development. Gene Expr Patterns. (2006) 6:964–70. doi: 10.1016/j.modgep.2006.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Rifat Y, Parekh V, Wilanowski T, Hislop NR, Auden A, Ting SB, et al. Regional neural tube closure defined by the Grainy head-like transcription factors. Dev Biol. (2010) 345:237–45. doi: 10.1016/j.ydbio.2010.07.017

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Pyrgaki C, Liu A, Niswander L. Grainyhead-like 2 regulates neural tube closure and adhesion molecule expression during neural fold fusion. Dev Biol. (2011) 353:38–49. doi: 10.1016/j.ydbio.2011.02.027

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Brouns MR, De Castro SC, Terwindt-Rouwenhorst EA, Massa V, Hekking JW, Hirst CS, et al. Over-expression of GRHL2 causes spina bifida in the Axial defects mutant mouse. Hum Mol Genet. (2011) 20:1536–46. doi: 10.1093/hmg/ddr031

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Werth M, Walentin K, Aue A, Schonheit J, Wuebken A, Pode-Shakked N, et al. The transcription factor grainyhead-like 2 regulates the molecular composition of the epithelial apical junctional complex. Development. (2010) 137:3835–45. doi: 10.1242/dev.055483

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Walentin K, Hinze C, Werth M, Haase N, Varma S, Morell R, et al. A GRHL2-dependent gene network controls trophoblast branching morphogenesis. Development. (2015) 142:1125–36. doi: 10.1242/dev.113829

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Wan X, O'Quinn RP, Pierce HL, Joglekar AP, Gall WE, DeLuca JG, et al. Protein architecture of the human kinetochore microtubule attachment site. Cell. (2009) 137:672–84. doi: 10.1016/j.cell.2009.03.035

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Karpen GH, Allshire RC. The case for epigenetic effects on centromere identity and function. Trends Genet. (1997) 13:489–96. doi: 10.1016/S0168-9525(97)01298-5

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Zhou CJ, Wang XY, Han Z, Wang DH, Ma YZ, Liang CG. Loss of CENPF leads to developmental failure in mouse embryos. Cell Cycle. (2019) 18:2784–99. doi: 10.1080/15384101.2019.1661173

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Ferreira RM, Chiaratti MR, Macabelli CH, Rodrigues CA, Ferraz ML, Watanabe YF, et al. The infertility of repeat-breeder cows during summer is associated with decreased mitochondrial DNA and increased expression of mitochondrial and apoptotic genes in oocytes. Biol Reprod. (2016) 94:1–10. doi: 10.1095/biolreprod.115.133017

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Hediger MA, Clemencon B, Burrier RE, Bruford EA. The ABCs of membrane transporters in health and disease (SLC series): Introduction. Mol Aspects Med. (2013) 34:95–107. doi: 10.1016/j.mam.2012.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Mathew DJ, Sanchez JM, Passaro C, Charpigny G, Behura SK, Spencer TE, et al. Interferon tau-dependent and independent effects of the bovine conceptus on the endometrial transcriptome. Biol Reprod. (2019) 100:365–80. doi: 10.1093/biolre/ioy199

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Ruprecht JJ, Kunji ERS. The SLC25 mitochondrial carrier family: structure and mechanism. Trends Biochem Sci. (2020) 45:244–58. doi: 10.1016/j.tibs.2019.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Gorgoglione R, Porcelli V, Santoro A, Daddabbo L, Vozza A, Monne M, et al. The human uncoupling proteins 5 and 6 (UCP5/SLC25A14 and UCP6/SLC25A30) transport sulfur oxyanions, phosphate and dicarboxylates. Biochim Biophys Acta Bioenerg. (2019) 1860:724–33. doi: 10.1016/j.bbabio.2019.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Forde N, Simintiras CA, Sturmey R, Mamo S, Kelly AK, Spencer TE, et al. Amino acids in the uterine luminal fluid reflects the temporal changes in transporter expression in the endometrium and conceptus during early pregnancy in cattle. PLoS ONE. (2014) 9:e100010. doi: 10.1371/journal.pone.0100010

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Tribulo P, Balzano-Nogueira L, Conesa A, Siqueira LG, Hansen PJ. Changes in the uterine metabolome of the cow during the first 7 days after estrus. Mol Reprod Dev. (2019) 86:75–87. doi: 10.1002/mrd.23082

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Chaney HL, Grose LF, Charpigny G, Behura SK, Sheldon IM, Cronin JG, et al. Conceptus-induced, interferon tau-dependent gene expression in bovine endometrial epithelial and stromal cells. Biol Reprod. (2021) 104:669–83. doi: 10.1093/biolre/ioaa226

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Haguenauer A, Raimbault S, Masscheleyn S, Gonzalez-Barroso MD, Criscuolo F, Plamondon J, et al. A new renal mitochondrial carrier, KMCP1, is up-regulated during tubular cell regeneration and induction of antioxidant enzymes. J Biol Chem. (2005) 280:22036–43. doi: 10.1074/jbc.M412136200

PubMed Abstract | CrossRef Full Text | Google Scholar