Dynamic mRNA and miRNA expression of the head during early development in bighead carp (Hypophthalmichthys nobilis)

Head of fish species, an exquisitely complex anatomical system, is important not only for studying fish evolution and development, but also for economic values. Currently, although some studies have been made on fish growth and body shapes, very limited information is available on the molecular mechanism of head development. In this study, RNA sequencing (RNA–Seq) and small RNA sequencing (sRNA–Seq) technologies were used to conduct integrated analysis for the head of bighead carp at different development stages, including 1, 3, 5, 15 and 30 Dph (days post hatch). By RNA-Seq data, 26 pathways related to growth and bone formation were identified as the main physiological processes during early development. Coupling this to sRNA–Seq data, we picked out six key pathways that may be responsible for head development, namely ECM receptor interaction, TNF signaling pathway, osteoclast differentiation, PI3K–Akt signaling pathway, Neuroactive ligand–receptor interaction and Jak–STAT signaling pathway. Totally, 114 important candidate genes from the six pathways were obtained. Then we found the top 20 key genes according to the degree value by cytohubba, which regulated cell growth, skeletal formation and blood homeostasis, such as pik3ca, pik3r1, egfr, vegfa, igf1 and itga2b. Finally, we also acquired 19 key miRNAs playing multiple roles in the perfection of various tissues in the head (such as brain, eye and mouth) and mineralization of head bone system, such as let–7e, miR–142a–5p, miR–144–3p, miR–23a–3p and miR–223. Results of this study will be informative for genetic mechanisms of head development and also provide potential candidate targets for the interaction regulation during early growth in bighead carp.


Introduction
The head of vertebrate is an extremely complex organ, which is one of the classic research topics of evolutionary morphology and developmental biology [1]. Species evolved into different head shapes and sizes (skull morphologies) to adapt to their living environment [2].
The evolutionary morphology and developmental biology of head have been reported in a variety of animals, such as dogs [3,4], bird [5,6], and ruminant [7,8]. Head shapes and sizes of fish are very important not only for the studies on fish evolution and development, but also for the economic value of aquaculture industry, which directly affect the fillet yield of economic fish [2]. For most aquaculture fish, smaller head and uniform body shape provide a higher percentage of meat yield. Therefore, selecting smaller head and streamlined body shape is of great significance in aquaculture [9]. Because skull morphology of fish is one of the commercial traits, it is very important to study the genetic regulation of head development. The ultimate goal of researchers is to produce fish with very large or very small head for ornamental or economic value [10]. Genome-wide association study (GWAS) and quantitative trait loci (QTL) mapping of head shapes and sizes were performed in fish species, including common carp [10,11], German mirror carp [12] and catfish [2]. Although some studies have been made on the genetic architecture of the shape and size of fish head, little is known about the genetic control basis of the head development in aquaculture fish.
Bighead carp (Hypophthalmichthys nobilis) is one of the most important cultured fishes in Asia, which provides a large amount of high-quality protein for human beings [13]. Principally produced in China, the global production of bighead carp reached 3.14 million tons in 2018 (FAO). For most aquaculture fish, smaller head and uniform body shape provide a higher percentage of meat yield [9]. Unlike other fish, the head of bighead carp is a raw material of many traditional Chinese dishes (Steamed Fish Head with Diced Hot Red Peppers, and so on), and the head is more popular than its fillet with even higher economic value. In addition, the nutrients of brain pulp, phospholipid and the soft skin with collagen in the head of bighead carp are rich, and therefore the head has higher nutritional value [14]. To meet the market demand for bighead carp, selection for fast-growing and big-headed varieties is prevalent among current breeding programs. At present, although some researchers have made some progresses on the genetic basis of fish growth and body shape, there are few reports on the molecular mechanisms of head development.
MicroRNAs (miRNAs) are small, endogenous and noncoding RNAs, about 22 nucleotides long; MiRNAs regulate gene silencing through partial messenger RNA (mRNA) excision and translation inhibition by binding to target mRNA [15]. Generally, a single miRNA regulates multiple target genes, and a single gene is also regulated by multiple miRNAs, so the interaction between miRNAs and target genes is complicated [16]. MiRNAs play important roles in regulation and control of growth, development, immunity and reproduction [17][18][19][20]. Benefiting from the advancement of sequencing technology, next-generation sequencing technology (NGS) provides powerful tools for studies on fish life activities and environmental adaptation, such as RNA sequencing (RNA-Seq) and microRNA sequencing (sRNA-Seq) [21]. Currently, the integrated analysis of mRNA and miRNA expression profiles has been performed on several fish species, including zebrafish [19,22], blunt snout bream [18], Japanese flounders [21,23], largemouth bass [24] and yellowfin seabream [20].
In the present study, we used NGST to conduct the RNA-Seq and sRNA-Seq from the head of bighead carp at different development stages, including 1Dph (days post hatch) group, 3Dph group, 5Dph group, 15Dph group and 30Dph group. This is the first report on integrated analysis of RNA-Seq and sRNA-Seq in bighead carp and as such offers deeper insight into the genetic mechanisms of head development. Furthermore, we identified a set of differentially expressed unigenes (DEGs) and differentially expressed miRNA (DEmiRs) from early growth, bone formation and blood homeostasis signaling pathways, which potentially services as valuable resources for future molecular breeding studies on head-related traits in fish species.

Sample characterization
To identify cues controlling molecular pathways associated with head development and skull morphology in bighead carp, we characterized differentially expressed mRNAs and miRNAs of the head in bighead carp samples from different development stages. To select suitable head tissue samples, the total length (TL), body length (BL), and head length (HL) of 15 bighead carp individuals in each development stage, including 1Dph (days post hatch) group, 3Dph group, 5Dph group, 7Dph group, 10Dph group, 15Dph group, 20Dph group, 30Dph group and 50Dph group, were measured in this study ( Table 1). The ratios of HL and BL (HL/BL) from these different development stages were calculated simultaneously. The HL/BL value showed an initial increase from 1Dph to 3Dph, then slowly increased from 3Dph to 15Dph, reached the peaking at Dph30 (Fig. 1). During the development stages from 3Dph to 5Dph, the HL/BL values were similar while the morphological characterizations changed greatly, such as the transformation of food habit from endogenous to exogenous, the progressive  Table S2. Based on the annotation, hundreds of genes and KEGG pathways related to head development and skull morphology were identified in the head transcriptome (Table 2), such as TGF-β signaling pathway (tgfb1), Insulin-like growth factors (Igf ) signaling pathway (igf1) and BMP pathway (bmp 3). A correlation analysis and principal component analysis (PCA) were performed on these 15 samples and the repeatability of samples from five time points was acceptable (Fig. 2).

Identification and function analysis of DEGs of head transcriptome during different developmental stages
Based on the criteria of |fold-change| ≥2 and FDR ≤ 0.01, the differentially expressed unigenes (DEGs) were identified of head transcriptome during different developmental stages. Comparison of adjacent development stages, Dph3vsDph1, Dph5vsDph3, Dph15vsDph5 and Dph30vsDph15 detected 7324, 833, 4238 and 16,061 DEGs, with 95 overlapping DEGs detected in all the four comparison groups (Fig. 3 (Table 3). Through Venn diagram analysis, 7 and 110 overlapped differentially expressed miRNAs were identified in the successive development stages and the following four development stages using Dph1 as reference stage, respectively (Fig. 4).

Analysis of miRNA-mRNA integration at five developmental stages
From the first day to the 30th day after hatching, significant morphological changes have taken place in the head of bighead carp, including the continuous perfection of various tissues in the head (such as brain, eye, mouth and other tissues) and the mineralization of the head skeletal system (brain skeleton and pharyngeal bone). Head development of bighead carp involves many vital movements, such as cell growth and death, organ development (such as immune system), signal transduction and nutrition metabolism (such as carbon metabolism), which indicates the complexity of head development of bighead carp.
In the significant enrichment pathways of different comparison groups, 10 pathways and 12 pathways overlapped with the above 26 pathways (Table 2) from transcriptome and small RNA data were shown in Fig. 5a and Fig. 5b, respectively. The common six pathways enriched by differential genes and target genes of differential miRNA (Fig. 5) were ECM-receptor interaction, TNF signaling pathway, osteoclast differentiation,  PI3K-Akt signaling pathway, neuroactive ligand receptor interaction and JAK-STAT signaling pathway. These six pathways were recognized as the key pathways of head development. A total of common 114 gene were obtained from differential genes and target genes of differential miRNA in the key pathways. Using string online  database, the protein interaction map of these genes was drawn (Fig. 6). Based on the protein interaction map, the top 20 key genes (hub genes of 114 candidate genes) were found according to the degree value by cytohubba, which regulated cell growth, skeletal formation and blood homeostasis (Fig. 7a). The heatmap of key genes was shown in Fig. S2. Finally, we constructed the miRNA-mRNA interaction network based on the 19 differential miRNAs whose target genes were the top 20 key genes using Cytoscape software (Fig. 7b).

Quantitative analysis of miRNA and mRNA expression
To validate the RNA-Seq and sRNA-Seq results, qRT-PCR was performed to detect the expression of 6 randomly selected differentially expressed mRNAs and 6 miRNAs. As shown in Fig. S3, the results of qRT-PCR were consistent with those of RNA-Seq and sRNA-Seq at 1Dph, 3Dph, 5Dph, 15Dph and 30Dph, respectively (Fig.  S3).

Discussion
Unlike other intensively-cultured fish, the head of bighead carp is a raw material of many traditional Chinese dishes, such as fish head with minced peppers, fish head hot pot, casserole fish head soup and fish head tofu soup. Therefore, the head is more valued than its fish meat [14]. To explore the molecular mechanisms of head development, we examined the RNA-Seq and sRNA-Seq data of bighead carp at different development stages, including 1Dph, 3Dph, 5Dph, 15Dph and 30Dph. Firstly, 26 pathways were identified as important pathways during early stages in the RNA-Seq data based on previous studies, which were closely related to growth development and bone formation [10,18,25,26]. Then the common six pathways enriched by differential genes and target genes of differential miRNA were obtained. The top 20 key genes were found from the 114 candidate genes in the six key pathways. Finally, we also acquired 19 key miRNAs playing multiple roles in the perfection of various tissues in the head (such as brain, eye and mouth) and mineralization of head bone system.

Key pathways regulating head development
The common six pathways enriched by differential genes and target genes of differential miRNA (Fig. 5) were ECM-receptor interaction, TNF signaling pathway, osteoclast differentiation, PI3K-Akt signaling pathway, neuroactive ligand receptor interaction and JAK-STAT signaling pathway. These six pathways were identified as the key pathways of head development in bighead carp. Wu et al. (2020) studied the transcriptome of breast muscle of Jinghai yellow chicken with different growth rate at early growth stage, and found that ECM -receptor interaction pathway was one of the enrichment pathways of differentially expressed genes, indicating that this pathway may play an important role in the growth and development of chicken [27]. As an important cytokine, tumor necrosis factor (TNF) can induce apoptosis, cell survival, inflammation and immunity. Wang and Guo (2021) studied the effects of TNF -α on the autophagy and NF -kappaB signaling pathway of fibroblast like synoviocytes in rheumatoid arthritis. The results showed that inhibition of autophagy can improve the imbalance of proliferation / apoptosis of fibroblast like synoviocytes aggravated by TNF -α to a certain extent, thus delaying the degree of rheumatoid arthritis [28]. In addition, Hurem et al.
(2017) also found TNF -α took part in regulation of γ radiation on the early development of zebrafish [29].
Osteoclasts are the multi-nuclear red cells from the mononuclear macrophage line, responsible for bone absorption [30], and participate in the development and regeneration of organs. PI3K -Akt signaling pathway is activated by a variety of cellular stimuli or toxic damage, and regulates basic cellular functions such as transcription, translation, proliferation, growth and survival. TNF signaling pathway, osteoclast differentiation and PI3K -Akt signaling pathway are involved in the formation of intermuscular bone in blunt snout bream, and play important roles in bone formation and development [25]. Similar to this study, neuroactive ligand receptor interaction was also found in the development of chicken crowns [31] and biomineralization and shell formation of pearl oyster Pinctada [32], indicating that this pathway plays an important role in animal development and biomineralization. JAK-STAT signaling pathway is one of the few pleiotropic cascade pathways, which is used for development and steady-state signal transduction in a variety of animals from human to fly [33]. In mammals, JAK-STAT signaling pathway is the main signaling mechanism of a series of cytokines and growth factors. Therefore, these six pathways play important roles in the growth, development and bone formation of animals. They may be involved in the continuous perfection of various tissues in the head (such as brain, eye, mouth and other tissues) and the mineralization of head bone system (cranial skeleton and pharyngeal bone).
Phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha (PIK3CA) and phosphoinositide-3-kinase regulatory subunit 1(PIK3R1) are members of PI 3-Kinases [34,35], which are involved in many cellular functions, including growth factor receptor signaling, cytoskeletal organization, apoptosis and so on [36]. Epidermal growth factor receptor (EGFR) is a kind of cell surface protein binding epidermal growth factor, which can promote cell proliferation by inducing receptor dimerization and tyrosine autophosphorylation [37]. Insulin like growth factor 1 (igf1) and hepatocyte growth factor (hgf) are closely associated with growth and development, which play important roles in cell growth, cell motility and morphogenesis [38,39]. Studies on igf1 gene have shown that its mutation can affect the growth and development and reduce body size and weight of animal, such as human, dog and mouse [3,40,41]. Mast/stem cell growth factor receptor kit (KIT) acts as a cell surface receptor. It regulates many essential biological processes, including cell survival and proliferation, hematopoiesis, stem cell maintenance and melanogenesis [42]. These kinases, growth factors and growth factor receptors may play important roles in the head development of bighead carp by regulating cell proliferation, cell growth, cell motility and morphogenesis.
Fibronectin 1 (FN1) is capable of binding cell surface and various compounds, such as fibrin, collagen, actin, DNA and heparin. Fn1 regulates cell adhesion, cell movement, wound healing and so on [43]. Moreover, fn1 plays an important role in the process of osteoblast compaction, which is necessary for osteoblasts mineralization [44]. Integrin alpha-IIb/beta-3 (ITGA2B:ITGB3) acts as a receptor for multiple proteins, such as fibronectin, prothrombin, and thrombospondin. It is involved in chondrocyte growth and proliferation [45]. Signal transducer and activator of transcription 3 (STAT3) protein changes the expression level of various genes to adapt to cell stimuli from many cytokines and growth factors, including HGF, EGF, BMP2 and IL5 [46]. It has been reported stat3 gene plays a key part in cell growth and apoptosis [47]. Stat1 gene is also identified in the head development of bighead carp, which is an important paralog of stat3. Studies on thrombospondin 1 (THBS1) protein have shown that THBS1 plays an important role in mediating cell-to-cell and cell-to-matrix interactions [48]. Nie et al. (2019) found stat1, thbs1 and thrombospondin at the early stage of intermuscular bone development, believing that they play an active role in osteoblast differentiation and intermuscular bone formation [25]. The identification of genes regulating bone formation suggested that development and mineralization of bone system coordinated with head development of bighead carp.
Vascular endothelial growth factor A (vegfa) encodes a heparin-binding protein, which exists as a disulfidelinked homodimer [49]. Kinase insert domain receptor (KDR), known as vascular endothelial growth factor receptor 2, acts as a cell-surface receptor for VEGFA, VEGFC and VEGFD [50]. Vegfa and kdr play essential roles in the regulation of angiogenesis, vascular development, vascular permeability, and embryonic hematopoiesis [51]. In this study, we also identified coagulation Factor II, thrombin (f2) [52], coagulation factor II thrombin receptor (f2r) [53], von willebrand factor (vwf) [54], janus kinase 1(jak1) [55] and glucagon receptor (gcgr) [56], which function in blood homeostasis, inflammation and wound healing. Platelet derived growth factor receptor beta (pdgfrb) is essential for normal development of the cardiovascular system and rearrangement of the actin cytoskeleton, and has roles in the regulation of many biological processes including embryonic development, angiogenesis, cell proliferation and differentiation [57]. Angiogenesis, vascular development, blood homeostasis and inflammation are involved in the tissue/organ development of animals. These genes associated with angiogenesis may regulate head development of bighead carp, providing blood homeostasis for the perfection of various tissues (brain, muscles, lipids, etc.) and the development of head bone system.

Key miRNAs regulating head development
Finally, we also acquired 19 key miRNAs (Fig. 7b) playing multiple roles in the regulation of cellular proliferation and differentiation, metabolism, immunity, regeneration and many others. let-7e and let-7i belonging to the let-7 family demonstrated to be associated with the development of cancer [58,59]. Hu et al. reported that let-7i was highly correlated with physiological and growth characteristics of velvet antlers [60]. MiR-223 is a hematopoietic cell-derived miRNA that plays multiple roles in regulation of monocytemacrophage differentiation, neutrophil recruitment, innate immunity, cancer and optic nerve regeneration [22,61]. Xie et al. have reviewed miR-223 exerts various effects on bone metabolism, especially in the processes of osteoclast and osteoblasts differentiation and highlighted the potential clinical applications of miR-233 in bone metabolism disorders [62].
Some studies have confirmed that the teleost-specific miR-462/731 cluster was involved in some biological processes, such as the immune responses to viral pathogens in Atlantic Salmon [71] and Atlantic cod [72], hypoxia [73] in zebrafish and blunt snout bream [74]. Huang et al. (2020) found significant reduction of digestive organs, especially liver and exocrine pancreas after the miR-462/miR-731 knockdown during zebrafish liver development, and those phenotypes could be partially rescued by corresponding miRNA duplex [75]. Besides these nine key miRNAs mentioned above and two new miRNAs which are not reported in the literature, the rest eight miRNAs also played roles in development, cell proliferation, cell cycle and tumorigenesis [76][77][78].
Key miRNAs identified in this study regulating a broad spectrum of genes involved in a variety of cellular processes, including development, cell proliferation, bone metabolism, immunity, tumorigenesis, and tissue regeneration, were thus considered to be related to the continuous development and perfection of various tissues in the head (such as brain, eye, mouth and other tissues) and the formation and mineralization of head bone system (cranial skeleton and pharyngeal bone). Further functional characterization of key genes and miRNAs determined by using over expression, knockout and knockdown strategies may help elaborate the molecular mechanisms for head development of bighead carp.

Ethics statement
All experimental protocols involved in fishes in this study were conducted in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the Institute of Hydrobiology and ARRIVE guidelines (https:// arriv eguid elines. org). This study was also approved by the Ethics Committee of Experimental Animals of the Hubei Provincial Committee for Animal Welfare. All efforts were made to minimize suffering of the fish.

Morphological observation and Fish sample collection
Fish samples used in this study originated from one mixed family of bighead carp, which were cultured at the Zhangdu Lake Fish Farm (Wuhan, China). In May 2019, the offspring were collected at the following developmental stages, including 1Dph (days post hatch) group, 3Dph group, 5Dph group, 7Dph group, 10Dph group, 15Dph group, 20Dph group, 30Dph group and 50Dph group. The total length (TL), body length (BL), and head length (HL) of 15 bighead carp individuals in each development stage were measured in this study. The ratios of HL and BL (HL/BL) from these different development stages were calculated simultaneously. Based on morphological characterizations and the HL/ BL values, the whole head in samples of five different developmental stages (1Dph, 3Dph, 5Dph, 15Dph and 30Dph) were selected for the integrated analysis of RNA-seq and sRNA-seq in bighead carp. After anesthetizing, the head of bighead carp from 1Dph, 3Dph and 5Dph were dissected by needles under a dissecting microscope (Zeiss Stemi 2000-C), while the head from 15Dph and 30Dph were dissected by scissors and tweezers. All head tissues (1Dph, 3Dph, 5Dph, 15 Dph and 30Dph) were sampled, quickly frozen in liquid nitrogen and stored at −80 °C prior to RNA extraction.

Total RNA extraction
Total RNA of mixed samples of head tissues during different developmental stages was extracted using TRIzol Reagent (Invitrogen, CA, USA) referring to a previous report [79]. The RNA degradation and concentration were first detected by 1% agarose gel electrophoresis and NanoDrop 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA), respectively. Then RNA quality was checked using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) and only those samples with the RNA Integrity Number (RIN) value ≥8.0 were used for RNA library construction.

Transcriptome sequencing and analysis
Fifteen head RNA-Seq libraries from three each of mixed samples during five different developmental stages (1Dph, 3Dph, 5Dph, 15 Dph and 30Dph) were constructed using NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) according to manufacturer's recommendations. Then libraries were subjected to paired-end sequencing of 150 bp on the Illumina HiSeq Xten platform by Novogene Company (Beijing, China). Before assembly, raw sequencing reads were clipped by removing adapter sequences and ambiguous nucleotides. Then all clean reads of the libraries of head tissues from five different development stages were assembled into transcripts by Trinity software, respectively. After assembly, the TGICL clustering software (J. Craig Venter Institute, Rockville, MD, USA) was used to cluster and remove redundant transcripts, and then the remaining sequences were defined as unigenes. All assembled sequences were used for BLAST searches and annotation against the NCBI non-redundant (Nr) protein, Swiss-Prot, gene ontology (GO), clusters of orthologous groups (COG), euKaryotic orthologous groups (KOG), eggNOG, and protein family (Pfam) databases using an e-value of 1e-5, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology results were obtained by comparing with KEGG database using KOBAS 2.0 [80]. After predicting the amino acid sequence of unigenes, the software HMMER [81] was used to compare them with Pfam database to get the annotation information of unigenes.
Gene expression levels were estimated by RSEM (RNA-Seq by Expectation Maximization) software package using a generative model of RNA-Seq reads and the EM algorithm for each sample [82]. Read count statistics of genes were shown in Table S3. The mapped reads were normalized according to fragment per kilobase of exon model per million mapped reads (FPKM) for each unigene between different development groups. Differentially expressed genes (DEGs) between the two groups were identified by the DEseq2 package (samples with three biological replicates) [83]. The false discovery rate (FDR < 0.01) adjusted by Benjamini-Hochberg method was adopted to get DEGs. DEGs were defined as by parameters of FDR < 0.01 and the absolute value of the log2 ratio ≥ 1. FPKM of DEGs from RNA-Seq was shown in Table S4. In order to determine the potential functions and metabolic pathways of these DEGs, DEGs were annotated against the Nr protein, Swiss-Prot, GO, COG, KOG, eggNOG, Pfam, and KEGG databases. KOBAS 2.0 software [80] was used to performed the statistical enrichment of DEGs in KEGG pathways. The principal component analysis (PCA) and heatmap analysis were conducted by R (Version 3.0.3) ggplot2 package and pheatmap package, respctively.

Small RNA library sequencing and analysis
Fifteen head sRNA-Seq libraries from three each of mixed samples during five different developmental stages (1Dph, 3Dph, 5Dph, 15 Dph and 30Dph) were generated using NEBNext ® Multiplex Small RNA Library Prep Set for Illumina ® (NEB, USA.) following manufacturer's recommendations. Then the fifteen libraries were sequenced on an Illumina Hiseq 2500 platform and 50 bp singleend reads were generated. Raw data of fastq format were firstly processed through custom perl and python scripts to obtain clean data. The clean sequence reads were mapped with miRBase 21.0, allowing a mismatch of one or two nucleotide bases, to identify the known miR-NAs. The potential miRNAs and secondary structures were obtained by srna-tools-cli with "srna-tools.pl --tool hp_tool --longSeq --shortSeqs --out" parameter [84]. The available software miREvo [85] and mirdeep2 [86] were integrated to predict novel miRNA through exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the small RNA tags unannotated in the former steps. Hairpin and mature sequences of novel miRNAs from sRNA-Seq was shown in Table  S5. Predicting the target gene of miRNA was performed by miRanda [87]. Read count statistics of miRNAs were shown in Table S6. The expression levels of the identified miRNAs were calculated and normalized to TPM (tags per million). Differential expression analysis of the miR-NAs was performed using the DESeq R package (v1.8.3) with default parameters, and miRNAs with a corrected P value (p-adjust/padj) < 0.05 were considered as differential expressed miRNAs (DEmiRs). TPM of DEmiRs from sRNA-Seq was shown in Table S7. In order to determine the potential functions and metabolic pathways of the target gene candidates of differentially expressed miR-NAs, GOseq [88] and KOBAS [89] were implemented for GO enrichment analysis and KEGG pathway enrichment analysis of the target gene candidates of differentially expressed miRNAs.

Quantitative PCR for miRNA and mRNA expression
In order to examine the reliability of the RNA-Seq and the sRNA-Seq results, 6 DEGs and 6 DEmiRs were randomly selected for validation by qRT-PCR method. Total RNA from 15 samples (three each of five different developmental stages) was extracted individually using TRIzol Reagent (Invitrogen, CA, USA) according to the manufacturer's instruction. 1 μg of total RNA for each sample was used to synthesize cDNA by using PrimeScript ™ RT reagent Kit (TaKaRa, Dalian, China) for DEGs and miRNA first-strand cDNA synthesis kit (Vazyme, Nanjing, China) with stem-loop reverse transcription primers for DEmiRs according to the manufacturer's protocols, respectively. Compared with stem-loop reverse transcription primers of DEmiRs, the reverse transcription primers of U6 reference gene can be directly designed by conventional primer premier 5 software. The downstream primer of U6 can be input to make reverse transcription according to the manufacturer's protocols (Vazyme, Nanjing, China). The PCR primers for DEGs and DEmiRs were list in Table S8. The qRT-PCR reaction solution was 20 μL, containing 10 μL of SYBR Green PCR Master Mix (TaKaRa, Dalian, China), 3 μL of forward and reverse primers (2 μmol/L), 2 μL of cDNA (500 ng) and 2 μL of sterile distilled water. qRT-PCR amplification was performed on a StepOne ™ Real-Time PCR System (Applied Biosystems, USA) according to a previous report [79]. RNA samples of head tissues from five different developmental stages were run in three biological replicates and three technical replicates for qRT-PCR. The expression level of each DEG and DEmiR was normalized to that of the reference gene β-actin and U6 by using the 2 −ΔΔCT method [90] to validate the results of RNA-seq and sRNA-seq, respectively.

Conclusions
This study represents the first application of RNA-Seq and sRNA-Seq in conducting integrated analyses for the head of bighead carp during different development stages, including 1Dph, 3Dph, 5Dph, 15Dph and 30Dph. Firstly, 26 pathways related to head development and bone formation were enriched and identified as the main pathways during early growth. Coupling this to sRNA-Seq data, we picked out six key pathways that may be responsible for head development, namely ECM receptor interaction, TNF signaling pathway, osteoclast differentiation, PI3K-Akt signaling pathway, Neuroactive ligand-receptor interaction and Jak-STAT signaling pathway. Totally, 114 important candidate genes were obtained from these six pathways. Then top 20 key genes were identified according to the degree value by cytohubba, which regulated cell growth, skeletal development and blood homeostasis, such as pik3ca, pik3r1, egfr, vegfa, igf1 and itga2b. Finally, 19 key miRNAs playing multiple roles in the perfection of various tissues in the head (such as brain, eye, muscle, lipids and mouth) and mineralization of head bone system, such as let-7e, miR-142a-5p, miR-144-3p, miR-23a-3p and miR-223 were also acquired. Results of this study will shed lights on further understanding of genetic mechanisms underlying head development and also provide potential candidate targets for the interaction regulation during early growth in bighead carp.