Differentially Expressed MiRNAs of Goat Submandibular Glands Among Three Developmental Stages Are Involved in Immune Functions

Submandibular glands (SMGs) are one of the primary components of salivary glands in goats. The proteins and biologically active substances secreted by the SMGs change with growth and development. Our previous studies showed that most of the differentially expressed genes in the SMGs of goats at different developmental stages are involved in immune-related signaling pathways, but the miRNA expression patterns in the same tissues are unknown. The aim of this study was to reveal the expression profile of miRNAs at three different developmental stages, detect differentially expressed miRNAs (DE miRNAs) and predict disease-related DE miRNAs. SMG tissue samples were collected from groups of 1-month-old kids, 12-month-old maiden goats and 24-month-old adult goats (three samples from each group), and high-throughout transcriptome sequencing was conducted. A total of 178, 241 and 7 DE miRNAs were discovered between 1-month-old kids and 12-month-old maiden goats, between 1-month-old kids and 24-month-old adult goats, and between 12-month-old maiden goats and 24-month-old adult goats, respectively. Among these DE miRNAs, 88 DE miRNAs with medium or high expression levels (TPM ≥50) were classified into five expression pattern clusters. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses indicated that some of the predicted target genes of the DE miRNAs in the five clusters were enriched in disease-related GO terms and pathways. MiRNA target genes in significant pathways were significantly enriched in Hepatitis B (FDR = 9.03E-10) and Pathways in cancer (FDR = 4.2E-10). Further analysis was performed with a PPI network, and 10 miRNAs were predicted to play an important role in the occurrence and prevention of diseases during the growth and development of goats.

Submandibular glands (SMGs) are one of the primary components of salivary glands in goats. The proteins and biologically active substances secreted by the SMGs change with growth and development. Our previous studies showed that most of the differentially expressed genes in the SMGs of goats at different developmental stages are involved in immune-related signaling pathways, but the miRNA expression patterns in the same tissues are unknown. The aim of this study was to reveal the expression profile of miRNAs at three different developmental stages, detect differentially expressed miRNAs (DE miRNAs) and predict disease-related DE miRNAs. SMG tissue samples were collected from groups of 1-month-old kids, 12-month-old maiden goats and 24-month-old adult goats (three samples from each group), and high-throughout transcriptome sequencing was conducted. A total of 178, 241 and 7 DE miRNAs were discovered between 1-month-old kids and 12-month-old maiden goats, between 1-month-old kids and 24-month-old adult goats, and between 12-month-old maiden goats and 24-month-old adult goats, respectively. Among these DE miRNAs, 88 DE miRNAs with medium or high expression levels (TPM ≥50) were classified into five expression pattern clusters. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses indicated that some of the predicted target genes of the DE miRNAs in the five clusters were enriched in disease-related GO terms and pathways. MiRNA target genes in significant pathways were significantly enriched in Hepatitis B (FDR = 9.03E-10) and Pathways in cancer (FDR = 4.2E-10). Further analysis was performed with a PPI network, and 10 miRNAs were predicted to play an important role in the occurrence and prevention of diseases during the growth and development of goats.

INTRODUCTION
The submandibular glands (SMGs) are one of the primary components of salivary glands in mammals. SMGs lie under the tongue and develop in set positions within the mouth (Tucker, 2007). In goats, the SMGs are mixed gland containing both serous and mucous acinar cells; in addition, the SMGs are the first major salivary gland to develop in the embryo. Mucous acinar cells mainly secrete mucins, which increase the viscosity of saliva, and serous acinar cells secrete a large number of proteins, which participate in physiological and biochemical reactions during growth and development (Kameyama et al., 2020).
The SMGs change with growth and development, and these changes may influence the secretion process and result in secretory granules with irregular shapes (D'Avola et al., 2006). Studies have shown that the proteins and biologically active substances secreted by SMGs change with age. The levels of hormones and epidermal growth factor (EGF) in rat SMGs gradually increase with age from 4 to 8 weeks (Hiramatsu et al., 1994). The level of carbonic anhydrase VI (CA VI) expression in ovine SMGs between days 7 and 59 is dramatically increased (∼600-fold), as determined by radioimmunoassay (Penschow et al., 1997). In the SMGs of postnatal rats, substance P increases in surges during the first 8 weeks of life (Ekström et al., 1994). In addition to the types of substances, their proportions in the SMGs change with growth and development. Kameyama et al. found that both types of mucins secreted in the mouse SMGs and the proportions of the mucin glycan species expressed throughout life change during the aging process (Kameyama et al., 2020).
MiRNAs, a class of non-coding RNA molecules with lengths of ∼22 nucleotides (nt), exert a negative regulatory effect on gene expression through the inhibition and degradation of target mRNAs by binding to target mRNAs through sequencespecific base pairing after gene transcription (Hudder and Novak, 2008). MiRNAs play an important role in cellular metabolism, growth, proliferation and the immune response. In cells, protein expression levels are regulated and further influenced by gene expression; some physiological and biochemical indicators of biological functions are regulated in this way. Previous studies have shown that miRNA expression exhibits certain tissue specificity, and miRNA expression is different in different developmental stages. Conserved and novel miRNAs were identified in dairy goat mammary glands (Ji et al., 2012(Ji et al., , 2017, goat endometrium during embryo implantation (Song et al., 2015), testes (Wu et al., 2014), skin, and hair follicles by Solexa sequencing (Yuan et al., 2013). Specific changes in urinary miRNA expression in male goats during the breeding season may be related to their high testosterone levels during breeding (Longpre et al., 2014). Serum miRNA expression profiles are different in Chongming white goats pre-and post-weaning (Liao et al., 2019).
MiRNAs 103 and 107 are potent regulators of glucocorticoid receptor (GR) function in hypoxia (Yang et al., 2020). Goat skeletal muscle satellite cell proliferation can be regulated by miR-101a (Li et al., 2015). Many miRNAs were found in exosomes and milk fat globules of human and cow milk; these miRNAs include miRNA-148a, which attenuates the expression of DNA methyltransferase 1, and miRNA-125b, which targets p53, which is considered the guardian of the genome (Melnik and Schmitz, 2017). In the goat mammary gland, fat metabolism is suppressed by miR-181b, which regulates multiple genes in the Hippo pathway and targets insulin receptor substrate 2 (IRS2) (Chen Z. et al., 2016). In goat mammary epithelial cells, fat metabolism is suppressed by the miR130b-mediated regulation of PPARγ coactivator-1α . Fatty acid metabolism is regulated by miR-148a, which cooperates with miR-17-5p, repressing PPARGC1A and PPARA in goat mammary epithelial cells (GMECs) (Chen et al., 2017). MiR-34c promotes apoptosis in goat male germline stem-cells (mGSCs) in a p53-dependent manner and reduces their proliferation (Li et al., 2013). In dairy goats, spermatogonial stem cell proliferation is regulated by Sirt1, which is regulated by miR-204 (Niu et al., 2016).
MiRNAs play a significant role in immune responses. MiRNAs are critical for lymphocyte development and differentiation, and a deficiency in miRNAs also affects T and B cell development (Xiao and Rajewsky, 2009). miRNA-deficient mice display severe immune deficiency, particularly impaired B cell and T cell responses (Turner and Vigorito, 2008;Roy and Sen, 2011). The miR-146 family was predicted to regulate the gene expression of TNF receptor-associated factor 6 and IL-1 receptor-associated kinase 1 by binding to sequences in the 3 ′ UTRs of these genes. UTRs inhibit the expression of associated reporter genes and Toll-like and cytokine receptors (Taganov et al., 2006). MiRNAs induced by proinflammatory stimuli, such as miR-146, miR-155, and miR-21, have been of great interest in investigations of inflammatory and immune responses. In the immune system, miR-155, which is predicted to regulate a broad range of inflammatory mediators, has emerged as a central regulator (Roy and Sen, 2011). Increased expression of miR-221 and miR-222, which regulate brahma-related gene 1 (Brg1), correlates with increased immunoparalysis and organ damage, and specific miRNAs may serve as biomarkers of immunoparalysis (Seeley et al., 2018).  were all detected in goat milk. miR-223 was found at higher levels in goat colostrum than in human colostrum, while the expression levels of miR-146 and miR-155 were absolutely opposite, and these miRNAs are significantly more highly expressed in colostrum than in mature milk (Na et al., 2015). In dairy cows with Staphylococcus aureus-induced mastitis, immune cytokines are regulated by miR-145, which targets and regulates FSCN1, and downregulating miR-145 and upregulating FSCN1 expression promotes mammary epithelial cell proliferation and accelerates the recovery of damaged tissues (Chen Z. et al., 2019). However, studies on the effects of miRNAs of SMGs were mainly focused on humans, mice and other animal models, the miRNA of goat SMGs is rarely reported.
In our previous study, the mRNA expression in the SMGs of goats were distinctly different at different developmental stages, and immune-related genes/pathways were identified . Numerous mRNAs, including mRNA encoding 8 antimicrobial peptides, were highly expressed in the SMGs of 1month-old kids and hardly expressed in 12-month-old maiden goats and 24-month-old adult goats, while the immunoglobulin expression pattern was the opposite, with low expression in the SMGs of 1-month-old kids and high expression in 12-monthold maiden goats and 24-month-old adult goats. Since miRNA data of the SMGs of goats at different developmental stages are lacking, it is difficult for us to study the relationship between miRNA and gene expression in this tissue. Therefore, in this study, the miRNA profiles of SMG samples from 1-month-old kid goats, 12-month-old maiden goats, and 24-month-old adult goats were described, and DE miRNAs were investigated. These studies are helpful for exploring the potential roles of miRNAs during goat development. The results expanded the repertoire of known goat miRNAs and may help to better understand the function and underlying mechanisms of goat SMGs.

Ethics Statement
All the animal experiments, which were conducted in accordance with the "Guidelines for Experimental Animals" of the Ministry of Science and Technology (Beijing, China), were approved by the Committee on the Management and Use of Laboratory Animals of Shandong Agricultural University (Permit Number: SDAUA-2018-052). Every effort was made to minimize pain.

Experimental Sample Collection
The experimental animals were obtained from Shandong Fengxiang Animal Husbandry and Seed Industry Technology Co., Ltd. (Laiwu District, Jinan, Shandong Province, China). In this experiment, goats were divided into three different groups based on their developmental stages: the 1-month-old kids (group A), the 12-month-old maiden goats (group B) and the 24-month-old adult goats (group C). Three healthy female goats (Lubo goats) were included in each group (nine goats in total). All the goats were raised under the same conditions with natural lighting and free access to food. Before slaughter, all the goats were fasted from food for 24 h and fasted from liquids for 2 h, and then, xylazine hydrochloride (Huamu Animal Health Products Co., Ltd., China) was injected intramuscularly to induce general anesthesia. After slaughter and exsanguination, the SMGs were immediately collected surgically, placed in liquid nitrogen and then placed in a −80 • C freezer for long-term storage. To distinguish the samples, the submandibular salivary gland samples from 1-month-old kids were labeled A1-s, A2-s and A3-s; the samples from 12-month-old maiden goats were labeled B3-s, B4-s and B5-s; and the samples from 24-month-old adult goats were labeled C2-s, C3-s and C5-s.

RNA Extraction, Library Construction, and Sequencing
Total RNA was extracted from SMGs tissues using TRIzol (Invitrogen, Carlsbad, CA, USA). The RNA molecules in a size range of 18-30 nt were then enriched by polyacrylamide gel electrophoresis (PAGE). The 3 ′ adapters were added, and the 36-44 nt RNAs were enriched. The 5 ′ adapters were then also ligated to the RNAs. The ligation products were reverse transcribed by PCR amplification, and the 140-160 bp PCR products were enriched to generate a cDNA library and sequenced using Illumina HiSeqTM 2500 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Sequencing Data Quality Control
To obtain clean reads, raw reads were further filtered with FastQC (Cock et al., 2010) according to the following rules: remove low-quality reads containing more than one low-quality (Q-value ≤20) base or containing unknown nucleotides (N); remove reads without 3 ′ adapters; remove reads containing 5 ′ adapters; remove reads containing 3 ′ and 5 ′ adapters but no small RNA fragment between them; remove reads containing polyA in small RNA fragments (More than 70% of the bases in a single read are A); and remove reads shorter than 18 nt (not including adapters).

Sequencing Data Alignment and miRNA Identification
With Bowtie (version 1.1.2) (Langmead and Salzberg, 2012), we aligned our clean reads with small RNAs in the GenBank database (Release 209.0) and Rfam database (11.0) to identify and remove rRNAs, scRNAs, snoRNAs, snRNAs, and tRNAs. We also mapped our clean reads with the goat reference genome (Capra hircus ARS1), and the reads that mapped to exons, introns and repeat sequences were removed. To identify known miRNAs, we mapped the clean reads with the miRBase database (Release 22). All of the unmapped reads were used for novel miRNA identification according to their genome positions and hairpin structures predicted by Mireap_v0.2 software (Chen L. et al., 2019). The parameters of Mireap_v0.2 were established as follows: minimal miRNA sequence length is 18 nt, maximal miRNA sequence length is 26 nt, minimal miRNA reference sequence length is 20 nt, maximal miRNA reference sequence length is 24 nt, minimal depth of Drosha/Dicer cutting site is 3, maximal copy number of miRNAs on reference is 20, maximal free energy allowed for a miRNA precursor is 18 kcal/mol, maximal space between miRNA and miRNA * is 35 nt, minimal space between miRNA and miRNA * is 14 nt, maximal bulge between miRNA and miRNA * is 4 nt, maximal asymmetry of miRNA/miRNA * duplex is 5 nt, and flank sequence length of miRNA precursor is 10 nt. To distinguish known miRNAs in goat and other species, the 5 p and 3 p in the names of miRNA were denoted by "x" and "y, " respectively.

MiRNA Expression Profile and PCA Analysis
For both known miRNAs and novel miRNAs, the miRNA expression level was calculated and normalized by transcripts per million (TPM). The formula was as follows: TPM = Actual miRNA counts/Total counts of clean tags * 106. The principal component analysis was performed using TPM of all miRNAs by R package (http://www.r-project.org/).

Differentially Expressed miRNA Analysis
DE miRNAs across different developmental stages were identified with DESeq2 (version 1.10.1) (Love et al., 2014). A fold change ≥2 and a false discovery rate value (FDR) <0.05 were set as cutoffs to determine significant differential expression. Then, we performed target gene prediction for all the DE miRNAs with RNAhybrid (v2.1.2) (Krüger and Rehmsmeier, 2006), Miranda (v3.3a) (Enright et al., 2003), and TargetScan (v7.0) (Sethupathy et al., 2006). The intersection results from three software programs were chosen as predicted miRNA target genes. For RNAhybrid, threshold were as follows: forces structures to have a helix from position 2-8 with respect to the query; the number of hits per target is 1; maximal bulge loop size is 3 nt; maximal internal loop size (per side) is 3 nt; the cutoff P-value is 0.05; the cutoff energy is −10 kcal/mol; and maximal query length is 24 nt. For Miranda, the set score threshold is 140; the set energy threshold is −10 kcal/mol; demand strict 5 ′ seed pairing; gap-open penalty is−4.0; and the gap-extend penalty is−9.0. For TargetScan, 2-8-nt sequences starting from 5 ′ small RNA were chosen as seed sequences to predict the 3 ′ UTR of the transcripts. After removing the miRNAs with low expression levels (TPM < 1), the DE miRNA expression scale of nine groups of samples was standardized by Z-Score.
Family analysis was conducted based on sequence similarity and the conservation of miRNAs among different species. According to the average TPM values of three samples in each developmental stage, the miRNAs were divided into four groups: extremely high-expression group (TPM ≥ 5,000), highexpression group (5,000 >TPM ≥ 500), medium-expression group (500 > TPM ≥ 50), and low-expression group (50 > TPM ≥1).

Expression Pattern Analysis of DE miRNAs
To further understand the expression pattern of the DE miRNAs with TPM ≥ 50 in the SMGs of goats at three different developmental stages, expression pattern analysis of the DE miRNAs was performed with TCseq software (1.14.0) (http://b ioconductor.org/packages/release/bioc/html/TCseq.html)  using the c-means method with the parameter dist = "Euclidean, " algo = "cm, " k = 5. Line graphs were drawn with the expression levels of the DE miRNAs of the goats in three different developmental stages.

Target Gene Prediction and Functional Enrichment Analysis
Schober et al. reported that−1.0 <Pearson correlation coefficient <-0.90 was strong or high correlation, and the correlation coefficient of 0.65 could either be interpreted as a "good" or "moderate" correlation (Schober et al., 2018). To generate a comprehensive description of the biological characteristics of the DE miRNAs, we mapped all of the candidate target genes to the Gene Ontology (GO) database (http://www. geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.genome.jp/kegg/) using the R package clusterProfiler (Alexa and Rahnenführer, 2009;Yu et al., 2012) with the screening standard FDR, calibrated after the pvalue. The terms and pathways with FDR <0.05 were considered significantly enriched pathways.

Network Analysis
PPI network analysis of the target genes in significantly enriched pathways was performed by Metascape (http://metascape.org/) and Cytoscape 3.7.2. Ten target genes were chosen as key genes according to interactions from large to small. According to a Pearson correlation coefficient <-0.7 and a p-value <0.05, the top 10 core gene targets of the miRNAs above were predicted. The top 10 miRNAs were selected as core miRNAs according to the order of miRNA expression and fold change in differential expression from large to small. Then, regulatory networks of the miRNAs and potential target genes that were significantly enriched in the disease-related KEGG pathways were drawn by Cytoscape 3.7.2.

Quantitative Real-Time PCR (qRT-PCR) Validation
To further verify the differential expression of the miRNAs in the nine goat SMG samples as determined by highthroughput sequencing, we randomly selected 20 differentially expressed miRNAs for RT-qPCR detection. The experiment was carried out using a LightCycler 480 real-time system (Roche, USA). The qRT-PCR primers used for validation are listed in Supplementary Table 1 and were designed according to the sequences of the selected DE miRNAs using Primer Premier 5.0 (Premier, Canada). qRT-PCR was performed according to the instructions of the Mir-X TM miRNA qRT-PCR TB Green TM Kit (Clontech). A reaction volume of 25 µl with 12.5 µl TB Green Advantage Premix (2×), 0.5 µl ROX Dye (50X), 0.5 µl miRNA-specific primer (10 µM), 0.5 µl mRQ 3 ′ Primer (10 µM), 2.0 µl cDNA, and 9 µl ddH 2 O was used. The reaction mixtures were incubated in a 96-well plate at 95 • C for 10 s, followed by 40 cycles of 95 • C for 5 s and 60 • C for 20 s, then 95 • C for 60 s, 55 • C for 30 s and 95 • C for 30 s. U6 was used as the reference gene. All the reactions were performed in triplicate, and the results were calculated using the 2 − Ct method (Livak and Schmittgen, 2001). Tukey's honestly significant difference (HSD) test was used to test the significance of mRNA expression, and miRNA expression histograms of different groups were generated by GraphPad Prism software Inc., San Diego, CA.

MiRNA Identification and Expression Profile
All the clean reads were aligned to the small RNAs in the miRBase 22 database; the alignment rates for the 9 libraries with goat known miRNAs were all over 60%, and the alignment rates for other known miRNAs were all over 10%. All of the unannotated reads were aligned with the reference genome. According to their genome positions and hairpin structures predicted by Mireap_v0.2 software, novel miRNA candidates were identified. The alignment rates for novel miRNAs in the nine libraries were in the range of 0.04-0.1% (Supplementary Table 2). All data have been uploaded to the Gene Expression Omnibus (GEO) database under the accession number PRJNA702107. According to the alignment results, miRNA annotation and identification were performed. In this study, 418 known miRNAs in goats, 1,292 known miRNAs in other species, and 479 new miRNAs were identified ( Table 1). As shown in Table 2, among the known miRNAs in goats, the expression level of 26, 55, 90, and 247 miRNAs reached extremely high expression (TPM ≥ 5,000), high expression (5,000 > TPM ≥ 500), medium expression (500 > TPM ≥ 500), and low expression (50 > TPM ≥ 1), respectively. Among the known miRNAs in other species, the number of miRNAs that reached extremely high expression, high expression, medium expression, and low expression was 9, 31, 60, and 1,192, respectively. For the novel miRNAs, only two miRNAs reached medium expression, and the other 477 miRNAs had low expression. Among all the miRNAs, the 10 miRNA families with the highest expression levels were miR-26, miR-148, let-7, miR-30, miR-99, miR-200, miR-125, miR-143, miR-27, and miR-181. The proportions of the expression of these miRNAs in each group are shown in Supplementary Figure 1a, and their sequences are shown in Table 3. The five miRNAs with the highest average expression levels were all goat known miRNAs, as shown in Table 4, and they were chi-miR-26a-5p, chi-miR-148a-3p, chi-miR-143-3p, chi-miR-99a-5p, and chi-miR-27b-3p. Together, their expression accounted for nearly one-third of the expression of all miRNAs (33.21%).
The results of miRNA length statistical analysis showed that the miRNAs in each of the three developmental stages ranged in length from 20 to 24 nt, of which the 22 nt miRNAs accounted for more than 40, and 47.08% of these miRNAs were expressed in the SMGs of 24-month-old adult goats (Supplementary Figure 1c). The PCA results showed that the first principal component (44.3%) could clearly distinguish the difference between group A and the other two groups (Supplementary Figure 1b), the results may indicate the similar expression of miRNAs between group B and group C.

Identification of Differentially Expressed miRNAs
The DE miRNA analysis revealed 178 and 241 DE miRNAs between group A and group B and between group A and group C, respectively, but only 7 DE miRNAs between group B and group C. These results further indicate the similar expression of miRNAs between group B and group C. It is worth noting that more than 80% of the DE miRNAs were significantly highly expressed in group A and expressed at low levels in groups B and C, and the results are shown in Figure 1. Conservation analysis and species source analysis were performed on DE miRNAs. As shown in Supplementary Figure 2a, the top 20 miRNA families of all DE miRNAs with high expression were let-7, miR-30, miR-9, miR-181, miR-92, miR-29, miR-124, miR-125, miR-10, miR-7, miR-2, miR-133, miR-34, miR-19, miR-1, miR-26, miR-2285, miR-199, miR-135, and miR-219. Conservation analysisof the DE miRNAs showed that most of the DE miRNAs matched the species of Homo sapiens, Bos taurus, Mus musculus, with the identified miRNA number of 417, 391, and 358, respectively (Supplementary Figure 2d). To goat known DE miRNAs, the top 20 miRNA families with high expression were different from the other species known DE miRNA  families. There were no miR-2 and miR-2285, except for miR-23 and miR-27 (Supplementary Figures 2b,c), and no Gallus gallus in the top 20 species sources, including Pan paniscus (Supplementary Figures 2e,f).

Expression Pattern Analysis of DE miRNAs
The results showed that all the DE miRNAs were significantly clustered among the different groups (Figures 2A,B). Among all the DE miRNAs, we focused on the miRNAs whose expression levels reached a medium and high level. 59 and 86 DE miRNAs were identified between group A and B and group A and C, respectively. Only chi-miR-206 met the requirement between group B and C, and this result may be related to the small changes that occur in goat SMGs after puberty. We identified a total of 88 DE miRNAs after consolidation, including 54 known miRNAs of goats, 26 known DE miRNAs of other species and eight novel miRNAs. Expression pattern analysis was performed to confirm the changes in the expression patterns of these miRNAs in the three different age groups ( Figure 2C). The results showed that the 88 DE miRNAs were grouped into 5 categories: cluster 1, cluster 2, cluster 3, cluster 4, and cluster 5. 7 miRNAs were included in cluster 1, and these miRNAs were expressed at low levels in group A and expressed at high levels in groups B and C. Cluster 2 also contained 7 miRNAs, and their expression gradually increased with growth and development. Ten miRNAs were included in cluster 3, and the expression trend was exactly opposite that of cluster 2; that is, the expression of these miRNAs gradually decreased with growth and development. Cluster 4 contained the largest number of miRNAs, 63, and these miRNAs were highly expressed in group A but expressed at low levels in groups B and C. Only one miRNA was included in cluster 5, that is, chi-miR-206, which was highly expressed in group B but expressed at low levels in groups A and C. The expression pattern results indicate that the expression of the DE miRNAs was time specific, and miRNAs with similar expression patterns may have similar functions.

Target Gene Prediction and Functional Enrichment Analysis
To analyze the functions of miRNAs with medium and high expression levels and with different expression patterns, the target genes of the DE miRNAs in goat SMGs were predicted with the standard of −1.0 < Pearson correlation coefficient <-0.68 and p-value < 0.05. Then, GO and KEGG pathway enrichment analyses of the miRNA target genes with different expression patterns were performed. GO enrichment analysis of the miRNA target genes with different expression patterns is shown in Figure 3. Among the significantly enriched terms of the five expression pattern clusters, more than half were enriched in biological processes. As shown in Figures 3B, 4 and Supplementary Table 3, most of the cluster 1-4 miRNA target genes that significantly enriched in the first three terms were different. These results suggested that the expression of differentially expressed genes in the SMGs of goats at 1, 12, and 24 months of age was regulated by miRNAs. With the growth and development of goats, the gene functions of SMG mRNAs were also altered by miRNA regulation.
The results of the KEGG pathway enrichment analysis of target genes are shown in Figure 5. For Cluster 1 and 2, most of the significantly enriched pathways were related to disease, and the three most significantly enriched pathways in Cluster 1 were Hepatitis B (FDR = 0.00054), Herpes simplex infection (FDR = 0.00054), HTLV-I infection (FDR = 0.00054); the most significantly enriched pathways in cluster 2 were Epstein-Barr virus infection (FDR = 0.0024), Hepatitis B (FDR = 0.014), and MiRNAs in cancer (FDR = 0.024). The pathways that were not only significantly enriched in cluster 1 but also significantly enriched in cluster 2 were Epstein-Barr virus infection, Hepatitis B, Herpes simplex infection, HTLV-I infection, miRNAs in cancer, and Pancreatic cancer. For cluster 4, the three most significantly enriched pathways were Protein processing in the endoplasmic reticulum (FDR = 0.0024), Valine, leucine, and isoleucine degradation (FDR = 0.0133) and Lysosomes (FDR = 0.0133). This may indicate that the DE miRNAs expressed at low levels in the SMGs of 1-month-old goats were related to the occurrence of infectious diseases, while the DE miRNAs expressed at high levels were related to protein formation processes and amino acid metabolism or other pathways.

Network Analysis
Protein-protein interaction analysis (PPI) was performed on 181 target genes that were significantly enriched in the KEGG pathways. The results are shown in Figure 6A. There were 612 interactions for 168 genes. All the genes in the network were ranked according to the degree of interaction, and the top 10 potential key target genes, namely, CDK6, CSNK2A1, EGLN3, MAPK9, PCNA, PRKCB, PTEN, RB1, CREBBP, and YWHAB, were predicted to be core genes. KEGG pathway enrichment analysis was further performed on the potential key target genes, and it was found that the significantly enriched pathways were Hepatitis B (FDR = 9.03E-10) and Pathways in cancer (FDR = 4.2E-10).
According to the screening criteria of Pearson correlation coefficient <-0.68 and p-value < 0.05, a total of 21 miRNAs were shown to target the 10 core target genes identified in the RNAseq data of goat SMG miRNAs and mRNAs. Network analysis was performed with Cytoscape 3.7.2 with the 10 core target genes and 21 miRNAs, and the results are shown in Figure 6B. The miRNAs with the most target genes were miR-141-y and chi-miR-141, both of which targeted four genes, namely, CSNK2A1, PCNA, PRKCB, and YWHAB, followed by chi-miR-20a-5p and FIGURE 4 | The venn diagram of target genes that significantly enriched in the first three GO terms in cluster 1-4. miR-29-y, which targeted three and two genes, respectively. The top 10 miRNAs were chosen as core miRNAs according to their expression levels, and their fold change from large to small with FDR < 0.05. The results are shown in Tables 5, 6. The 10 miRNAs were miR-29-y, miR-141-y, miR-142-x, chi-miR-93-5p, chi-miR-20b, chi-miR-20a-5p, chi-miR-17-5p, chi-miR-16b-5p, chi-miR-15b-5p, and chi-miR-130b-3p. Among the 10 miRNAs, miR-29-y and miR-141-y were gradually overexpressed in three developmental stages, while the other eight miRNAs were all highly expressed in 1-month-old goat SMGs and expressed at low levels in 12-and 24-month-old goat SMGs. This result may suggest that the miRNAs highly expressed in 1-month-old goat SMGs are related to diseases, but the specific regulatory mechanism requires further verification.

Experimental Validation
To verify the accuracy of the sequencing and analysis results, 20 DE miRNAs were selected for RT-qPCR validation (Figure 7).
The qRT-PCR validation results for most of the selected miRNAs were consistent with the RNA-seq results.

DISCUSSION
In this study, we performed miRNA sequencing on the SMGs of goats at three different developmental stages. Currently, functional studies on salivary gland miRNAs are mainly focused on humans and mice. A number of studies have already investigated the role of multiple salivary gland miRNAs in the etiopathogenesis of Sjögren's syndrome (Argyropoulou et al., 2018;Reale et al., 2018). Mature miRNAs in the fetal salivary gland have also been found to affect the development of other tissues and organs by being loaded into exosomes. The above studies indicate that miRNAs play important roles in salivary gland development and functional maintenance. However, currently, reports on miRNA function in the salivary glands of ruminants are rare. To the best of our knowledge, this   is the first miRNA sequencing and expression profiling study on goat and ruminant salivary glands. The overexpression or underexpression of miRNAs affects gene expression, thus affecting the related pathways that may cause lesions. miRNAs are considered biomarkers of disease due to their ability to interfere with physiological or pathological processes. Some diseases can be prevented and diagnosed by the detection of miRNAs. As vertebrate-specific miRNAs, miR-26 family members have been detected in many different kinds of tumors and play important roles in cardiovascular disease, Parkinson's disease, neurodegenerative diseases, and many other non-tumor diseases (Gao and Liu, 2011;Icli et al., 2014;Goh et al., 2019;Su et al., 2019). The validated target genes of these miRNAs are involved in cell metabolism, proliferation, differentiation, apoptosis, autophagy, invasion, and metastasis (Lagos-Quintana et al., 2002;Li et al., 2020). In some tumors, the differential expression pattern of miR-26 might be related to specific tumorigenic processes. miR-125a-5p significantly inhibits the proliferation of head and neck cancer cells, whereas inhibition of miR-125a-5p enhances their proliferation (Xu et al., 2020). The expression of miR-181 family members is significantly related to colorectal cancer and correlated with the lymph node metastasis of oral squamous cell carcinoma, and miR-181 may play a role in the pathogenesis of Sjögren syndrome. The let-7 family was one of the first tumor suppressor miRNAs to be identified (Boyerinas et al., 2010). The upregulation of let-7 expression is associated with the progression and poor prognosis of hepatocellular carcinoma . let-7 as a regulator of myofibroblast differentiation through the regulation of IL-6 in mice (Di et al., 2021). Jiang et al. found that let-7 is expressed in a variety of immune cells, including B cells and macrophages, which indicates that let-7 may play physiological and pathological roles in vivo (Jiang, 2019). In human and murine glioblastoma/glioma, differentially expressed let-7 miRNAs are signaling molecules that induce microglial inflammatory cytokine release, modulate antigen presentation, and attenuate cell migration in a TLR7-dependent manner (Buonfiglioli et al., 2019;Song et al., 2020). During cardiac hypertrophy, the overexpression of miR-99 facilitates the transition from physiological hypertrophy to pathological hypertrophy, and silencing of miR-99 has the opposite effect, through the regulation of the Akt-1 pathway (Ramasamy et al., 2018). Renal clear cell carcinoma-derived miR-27aloaded exosomes inhibit secreted frizzled-related protein 1 (SFRP1) expression and accelerate tumor angiogenesis (Hou et al., 2020). In Cr(VI)-induced carcinogenesis, miR-143 inhibits tumor growth and angiogenesis by regulating IL-6/HIF-1α in vivo . The miRNA-30 and miRNA-199 profiles in serum and tissue saliva samples distinguish benign from malignant salivary gland tumors (Cinpolat et al., 2017).
In the SMGs of goats, most of these top 10 miRNA families are related to diseases, especially cancer. Among them, let-7, miR-30, miR-181, and miR-125 are differentially expressed RNAs, which may indicate their potential functions in SMGs with growth. This requires further study. Among the top 10 miRNA families, two are involved in some metabolic processes. miR-27a significantly affects the expression of mRNAs related to milk fat metabolism (Lin et al., 2013). miR-148a-3p inhibits alpaca melanocyte pigmentation by targeting MITF .
In the present research, chi-miR-26a-5p, chi-miR-148a-3p, chi-miR-143-3p, chi-miR-99a-5p, and chi-miR-27b-3p were found to be the top five most highly expressed miRNAs. All five most highly expressed miRNAs were not differentially expressed among the three developmental stages, which may suggest that these miRNAs play important and stable roles in the SMGs of goat with growth. MiR-26a-5p was reported to be significantly downregulated in the B lymphocytes purified from primary Sjögren's syndrome patients (Wang-Renault et al., 2018). MiR-26a promoted endometrial epithelium cells proliferation and induced stromal cells apoptosis via the PTEN-PI3K/AKT pathway in dairy goats (Zhang et al., 2018). Dey et al. found that miR-26a is required for skeletal muscle differentiation and regeneration in vivo of mice. However, the function of miR-26a in the SMGs has not been reported (Dey et al., 2012). MiR-148a can negatively regulate cell proliferation, migration, and invasion by downregulating MTA2 in salivary adenoid cystic carcinoma cells (Li et al., 2018). In goat rumen, the expression of miR-148a-3p was significantly different between stages, and involved in the development of the rumen epithelial cells by targeting QKI, meanwhile proliferation of GES-1 cells were induced by miR148a-3p inhibitor (Zhong et al., 2020). So we hypothesize the stable expression of chi-miR-148a-3p in goat SMGs is also related to the growth and development of SMGs cells, which need further study. MiR-143-3p can affect epithelial cell homeostasis by targeting SERCA2b, RyR2, and AC9 in Sjogren's syndrome (Cortes et al., 2018). Chi-miR-143-3p influences the development of mammary gland epithelial cells in dairy goats via target Ndfip1 in vitro (Ji et al., 2019). MiR-99a-5p was found to play a role in the metastasis of salivary adenoid cystic carcinoma (Feng et al., 2017). Therefore, all five most highly expressed miRNAs except miR-27b-3p have already been reported to be involved in cell proliferation and apoptosis or homeostasis to maintain the salivary gland function. These miRNAs may play were the samples for goat of 1-month-old kids (group A), B were the samples for goat of 12-month-old maiden goats goats (group B), C were the samples for goat of 24-month-old adult goats (group C). Significant differences are indicated by * and ** (* indicates p < 0.01, ** indicates p < 0.05), NS indicates there is no significant differences.
important roles in the salivary glands of goats and are not related to the differential regulation of transcriptional levels among different developmental stages. Their underexpression may cause pathological or physiological changes and lead to diseases.
DE miRNA analysis showed that there were more than 150 DE miRNAs between the SMGs of group A and group B and the SMGs of group A and group C; in addition, more than 80% of the DE miRNAs were significantly downregulated, while there were only seven DE miRNAs between the SMGs of group B and group C. The results indicated that miRNA expression in the SMGs significantly changed during the growth and development of goats from kids to adolescent goats, while miRNAs were steadily expressed in the SMGs from adolescent goats to adult goats. This suggested that miRNAs play a stable role in maintaining salivary gland function in goats after puberty.
To investigate the functions of the DE miRNAs with different expression profiles, expression pattern analysis was performed on the medium and highly expressed DE miRNAs, and five clusters were identified. The expression and profiles of the DE miRNAs were correlated with the developmental stages of the goats. GO and KEGG pathway enrichment analyses were performed on target genes in the five clusters. The results showed that the target genes of the miRNAs in the different clusters were significantly enriched in different GO terms and KEGG pathways. The target genes of the DE miRNAs, expressed at low levels in the SMGs of kids and highly expressed in the SMGs of adolescent and adult goats (clusters 1 and 2) were significantly enriched in infectious diseases, including Hepatitis B, Herpes simplex infection, and HTLV-I infection. The growth and development of the SMG is carried out under epithelial-mesenchymal interactions (Wang and Laurie, 2004), which are regulated by a complex of signaling molecules, including growth or differentiation factors, cell receptors and extracellular matrix (ECM) proteins. Activation of hepatic stellate cells, exaggerated expression of profibrogenic genes, and/or suppression of antifibrogenic genes can cause excessive and persistent accumulation of the ECM, which further leads to fibrosis (Wynn, 2008;Loureiro et al., 2020). This indicates that the target genes of the DE miRNAs in clusters 1 and 2 may be related to the growth and development of SMG ECM proteins and that their abnormal expression can cause infectious diseases. Epithelial-to-mesenchymal transition (EMT) is implicated in a wide array of malignant behaviors in cancer, including proliferation, invasion, and metastasis (Pan et al., 2021). Many studies have reported that miRNAs play a crucial role in the malignant behaviors of cancer cells (Volinia et al., 2010;Gianpiero and Carlo, 2013), including EMT-related cancer metastasis (Aiello and Kang, 2019), and some miRNAs have been predicted to be potential tumor biomarkers. We noticed that cancer-related pathways were significantly enriched by the target genes of the cluster 1 and 2 DE miRNAs of goat SMGs. Therefore, the cancer-related DE miRNAs in the goat SMGs may have potential diagnostic, prognostic, and therapeutic value in cancers, and this requires a better understanding. The target genes of the DE miRNAs that are highly expressed in the SMGs of kids and expressed at low levels in the SMGs of unbred adolescent and adult goats (cluster 4) were significantly enriched in metabolismrelated and secretion-related pathways. Therefore, we speculated that the function of the goat SMGs changes after puberty compared with that at kids, and the SMGs mainly function in metabolism and secretion at later stages, which might be related to the gradual improvement of the goat immune system with growth and development (Jaskoll et al., 2001;Ghezzi and Ship, 2003). The results were basically consistent with our previous transcriptome study of goat SMGs, in which the genes that were highly expressed in the kid SMGs were significantly enriched in immune-related pathways . Therefore, we predicted that the biological function of goat SMGs changes with growth and development. The mechanism by which RNAs function in the goat SMGs is complex, and genes play their roles through the combined action of multiple types of RNAs.
Ten core miRNAs were screened by PPI analysis. The 10 core target genes were significantly enriched in Pathways of Hepatitis B (FDR = 9.03E-10) and Pathways in cancer (FDR = 4.2E-10), suggesting that miRNAs in the SMGs may be associated with the occurrence and progression of cancers and infectious diseases. Target gene prediction was performed on the 10 predicted core miRNAs, and the results showed that the expression of the target genes CSNK2A1, PCNA, PrKCB, YWHAB, and CDK6 could be regulated by miRNAs expressed at low levels in the 1-month-old kids, namely, chi-miR-141, miR-141-x, miR-141-y, and miR-29-y. As a known target for tumor, casin kinase 2α1 (encoded by CSNK2A1) play key roles in cell cycle regulation, cellular differentiation, proliferation and apoptosis regulation (Tsuyuguchi et al., 2019). It's activation is required for transforming growth factor β-induced epithelialmesenchymal transition (Kim et al., 2018). Proliferating cell nuclear antigen (PCNA) plays an essential role in controlling a variety of reactions as a component of the replication and repair machinery (Kelman, 1997;Maga and Hübscher, 2003). As a part of a transcription complex, CDK6 induces the expression of tumor suppressors (Kollmann et al., 2013;Sherr et al., 2016). Protein kinase Cβ (PRKCB) plays a role in autophagy regulation (Patergnani et al., 2013), CREBBP inactivation promotes the development of HDAC3-dependent lymphomas (Jiang et al., 2017), and abnormal expression of Casein kinase 2α1 (CSNK2A1) is linked to neurological diseases and cancer (Trinh et al., 2017). We predicted that chi-miR-141, miR-141-x, miR-141-y and miR-29-y mainly play a role in regulating the proliferation and differentiation of SMG cells. Our results also indicate that the growth and development of goat SMGs are basically completed before puberty. The expression of the target genes MAPK9, CREBBP and EGLN3 could be regulated by the miRNAs chi-miR-15b-5p, chi-miR-16b-5p, chi-miR-17-5p, chi-miR-20a-5p, chi-miR-93-5p, and miR-142-x, which were highly expressed in the SMGs of 1-month-old kids and expressed at low levels in the SMGs of 12-month-old unbred adolesent goats and 24month-old adult goats. JNK2 (encoded by MAPK9) is implicated in apoptosis (Sabapathy and Wagner, 2004;Dhanasekaran and Reddy, 2008;Wang et al., 2010). EGLN3 and CREBBP also play important roles in the occurrence and development of disease, and EGLN3 can promote cell apoptosis (Lee et al., 2005). The loss of CREBBP promotes the development of HDAC3-dependent lymphoma (Jiang et al., 2017). We predicted that chi-miR-15b-5p, chi-miR-16b-5p, chi-miR-17-5p, chi-miR-20a-5p, chi-miR-93-5p, and miR-142-x play a role in cell apoptosis of goat SMGs. Our results indicate that the fuction of goat SMGs was basically stable after puberty. Apoptosis has been considered as one of the main factors that may be related to loss of salivary gland secretory function (Hayashi, 2011). Apoptosis and inflammation are inseparably linked (Cullen et al., 2013).
In conclusion, a set of DE miRNAs in the SMGs of goats at different developmental stages were identified, and functional enrichment analysis showed that the target gene of them are mainly focused on cell proliferation, growth and apoptosis, and significantly related to disease. The DE miRNAs identified in our research might contribute to disease prevention and be involved in the gene regulation of submandibular glands at different developmental stages. Our results also indicated that the growth and development of goat SMGs are basically completed before puberty.

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: NCBI [Accession: PRJNA702107].

ETHICS STATEMENT
The animal study was reviewed and approved by Committee of the Management and Use of Laboratory Animals of Shandong Agricultural University (Permit Number: 2,007,005).

AUTHOR CONTRIBUTIONS
JW, TC, and AW contributed to conception and design of the study. XZ and LH organized the database. RX performed the statistical analysis. AW wrote the first draft of the manuscript. QL and YC wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

ACKNOWLEDGMENTS
We thank Shandong Fengxiang Animal Husbandry and Seed Industry Technology Co., Ltd., for providing the experimental goats. We thank the editors and reviewers for their contributions to this work.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.