Yulong Chen1,2†
Shiming Wan1,2†Qing Li3Xiaoru Dong1,2Jinghan Diao1,2Qing Liao1,2Gui-Ying Wang3
Ze-Xia Gao1,2,4*- 1Key Lab of Agricultural Animal Genetics, Breeding and Reproduction of Ministry of Education, Key Lab of Freshwater Animal Breeding, Ministry of Agriculture, College of Fisheries, Huazhong Agricultural University, Wuhan, China
- 2Engineering Research Center of Green Development for Conventional Aquatic Biological Industry in the Yangtze River Economic Belt, Ministry of Education, Wuhan, China
- 3Fisheries Research Institute, Wuhan Academy of Agricultural Sciences, Wuhan Xianfeng Aquaculture Technology Co. Ltd, Wuhan, China
- 4Engineering Technology Research Center for Fish Breeding and Culture in Hubei Province, Wuhan, China
Intermuscular bone (IB) occurs in the myosepta of teleosts. Its existence has an adverse influence on the edible and economic value of fish, especially for aquaculture species belonging to Cypriniformes. The growth mechanism of IBs is quite lacking. In this study, we firstly used single molecular real-time sequencing (SMRT) technology to improve the draft genome annotation and full characterization of the transcriptome for one typical aquaculture species, blunt snout bream (Megalobrama amblycephala). The long non-coding RNA (lncRNA), microRNA (miRNA), and messenger RNA (mRNA) expression profiles in two IB growth stages (1 and 3 years old) were compared through transcriptome and degradome analyses. A total of 126 miRNAs, 403 mRNAs, and 353 lncRNAs were found to be differentially expressed between the two stages. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that the significantly upregulated map2k6 and cytc in the MAPK/p53 signaling pathway and the significantly downregulated lama3 and thbs4b in the extracellular matrix (ECM)–receptor pathway may play a key regulatory role in IB growth. Bioinformatics analysis subsequently revealed 14 competing endogenous RNA (ceRNA) pairs related to the growth of IBs, consisting of 10 lncRNAs, 7 miRNAs, and 10 mRNAs. Of these, dre-miR-24b-3p and dre-miR-193b-3p are core regulatory factors interacting with four lncRNAs and three mRNAs, the interaction mechanism of which was also revealed by subsequent experiments at the cellular level. In conclusion, our data showed that IBs had higher activity of cell apoptosis and lower mineralization activity in IB_III compared to IB_I via interaction of MAPK/p53 and ECM–receptor signaling pathways. The downregulated zip1 interacted with miR-24a-3p and lnc017705, decreased osteoblast differentiation and Ca2+ deposition in the IB_III stage. Our identified functional mRNAs, lncRNAs, and miRNAs provide a data basis for in-depth elucidation of the growth mechanism of teleost IB.
Introduction
Intermuscular bone (IB), which only occurs in the myosepta of lower teleosts, is a small spicule-like bone derived from tendons (Danos and Ward, 2012). Most freshwater-farmed fishes have a substantial amount of IB, which has an extremely negative impact on their edible and economic value. Earlier studies have revealed the number and complex morphological characteristics of IB in different species (Dong et al., 2006; Wan et al., 2014). Subsequent research also revealed the possibility of reducing IB numbers based on ploidy change (Li et al., 2013), distant hybridization (Jiang et al., 2016), and genetic breeding (Cao et al., 2015; Xiong et al., 2018). It is worth noting that an IB-deficient grass carp (Ctenopharyngodon idellus) mutant and some specimens of tambaqui (Colossoma macropomum) without IB (normal individuals possessing significant numbers of IB) were found in an artificial gynogenetic group and a culture population (Xu et al., 2015; Perazza et al., 2016), respectively, which strongly indicated the feasibility of genetic improvement of IB numbers. Therefore, clarification of the regulatory mechanism of IB growth will contribute to trait improvement in aquaculture.
The biological phenotypic diversity was dominated by complex molecular mechanisms involving DNA modification and the regulation of transcription and translation, in which long non-coding RNA (lncRNA) and microRNA (miRNA) also play important roles through regulating gene expression. Numerous studies have paid attention to the biological function of lncRNA and revealed its important roles in cell differentiation, transcription regulation, and development (Yang et al., 2014; Larsson et al., 2015; Liu et al., 2015; Ali and Soudeh, 2017). In fish, lncRNA has also been reported to participate in many biological regulatory processes, including early sex differentiation (Cai et al., 2019; Yuan et al., 2019), skin color regulation (Luo et al., 2019), and immune regulation (Wang M. et al., 2018). miRNA can regulate transcript levels by binding the 3′ untranslated region (UTR) of the target mRNA to participate in various biological processes, such as osteoblast differentiation, development, disease, gene transcription, and translation (Hobert, 2008). So far, the miRNA-related studies had been reported in several fishes including zebrafish (Danio rerio), blunt snout bream (Megalobrama amblycephala), Botia superciliaris, and Pelteobagrus vachelli. These studies mainly pay attention to the regulatory role of miRNA in the growth (Zhao et al., 2017) and development of organs including the skeleton, gonad, skin, and liver (Gan et al., 2016; Zhang et al., 2016; Lan et al., 2019; Zhou et al., 2019). Nonetheless, there are few miRNA studies on fish IB development. In our previous studies, the dynamic expression patterns and potential function of mRNA and miRNA within IB development were revealed in the typical aquaculture species blunt snout bream (Wan et al., 2016). Meanwhile, miRNA transcriptome features of IB and surrounding connective tissues had also been investigated in M. amblycephala (Wan et al., 2015), which provided theoretical support for revealing the potential molecular regulatory mechanisms of IB development. Massive evidences had shown that lncRNA and miRNA play an important regulatory role in animal bone development, but the specific regulatory function of lncRNA and miRNA in fish IB development is still unclear.
The single molecular real-time sequencing (SMRT) technology can provide important information for improving draft genome annotation and the full characterization of the transcriptome for corresponding species (Sharon et al., 2013; Gao et al., 2016). This technology has also been widely used in fishes, by which Gong et al. (2018) constructed the first high-quality chromosome-level genome assembly in yellow catfish (Pelteobagrus fulvidraco), which offered a valuable reference for functional genomics studies of yellow catfish to decipher its economic traits and sex determination. The genome assembly and gene evolution analysis of goldfish (Carassius auratus) also provide an important resource for understanding the causes of goldfish variants (Chen et al., 2017). With the development of high-throughput sequencing technology, transcriptome sequencing has been utilized in various fields and species. However, the regulatory networks that determine specific traits are often extremely complex; it is difficult to accurately discover key regulatory factors from a single omics data. Fortunately, the development of multi-omics research has overcome this shortcoming, and the competing endogenous RNA (ceRNA) research has also become a new hot spot in revealing the mechanism of RNA interaction.
In this study, we firstly used SMRT technology to improve the draft genome annotation and full characterization of the transcriptome for M. amblycephala. Then, the lncRNA, mRNA, and miRNA expression profiles of IB samples in 1- and 3-year old specimens were revealed, respectively. Based on conjoint analysis of multi-omics at different periods, the lncRNA–miRNA–mRNA ceRNA regulatory network related to IB growth was detected and verified. This research aims to reveal the key lncRNA, miRNA, and target mRNA involved in teleost IB growth and provides an important reference for further revealing the molecular regulatory mechanism of IB growth (Figure 1).
Figure 1. Phenotypic characteristics of IB and research scheme diagram.
Materials and Methods
Experiment Materials and Sample Preparation
Before sampling, experimental M. amblycephala were anesthetized in well-aerated water containing the 100 mg/L concentration of tricaine methanesulfonate. IB samples of three 1-year-old (IB_I, IBs growing rapidly at this stage) and three 3-year-old (IB_III, IBs growing quite slowly at this stage) fishes were collected and pooled, respectively, with three biological replicates. Total RNA was isolated using RNAiso Plus Reagent (Takara, Japan) following the manufacturer's protocol. RNA degradation and contamination were monitored on 1% agarose gels. RNA concentration was measured using the NanoDrop 2000 spectrophotometer. Qualified RNA samples are stored at −80°C and then used for RNA sequencing. All animals and experiments were conducted in accordance with the “Guidelines for Experimental Animals” of the Ministry of Science and Technology (Beijing, China). The study was approved by the Institutional Animal Care and Use Ethics Committee of Huazhong Agricultural University. All efforts were made to minimize animal suffering.
Library Construction and Sequencing
An equal mixture of RNA from all samples of two periods was used for the Pacific Bioscience (PacBio) library construction. The mRNA containing polyA was enriched with Oligo(dT), which was then reversely transcribed into cDNA using the Clontech SMARTer PCR cDNA synthesis kit. After PCR amplification, DNA damage repair, end repair, SMRT adapter ligation, size fractionation, and selection (<4 and >4 kb), two SMRT bell libraries were constructed, and the combined SMRT bell library was then sequenced on the PacBio Sequel System.
lncRNA and small RNA sequencing libraries were generated using NEBNext® Ultra™ Directional RNA Library Prep Kit (NEB, USA) and NEBNext® Multiplex Small RNA Library Prep Set (NEB, USA), respectively, following the manufacturer's recommendations. The PCR products were purified (AMPure XP system), and library quality was assessed on the Agilent Bioanalyzer 2100 system. After cluster generation of the index-coded samples using TruSeq PE Cluster Kit v3-cBot-HS (Illumina), the libraries were sequenced on an Illumina HiSeq 2500 platform. Then 125 bp paired-end reads and 50 bp single-end reads were generated in lncRNA and small RNA libraries, respectively.
Correction, Annotation, and Structural Analysis of Full-Length Transcripts
The adaptor and raw reads with lengths <50 bp were removed to generate subreads, which subsequently were processed using the SMRT Link software to generate a circular consensus sequence (CCS). Subsequently, the full-length non-chimeric (FLNC) reads with complete 5′ primer, 3′ primer, and poly(A) tails of the same transcript were clustered using the iterative clustering for error correction (ICE) algorithm to remove redundancy and obtain consensus reads, and the non-full-length (NFL) sequences were further used to correct consensus reads using the arrow algorithm to obtain polished consensus reads, which were further corrected and improved via the Illumina RNA-Seq data (NCBI BioProject: PRJNA329421) generated by our lab and LoRDEC software (Salmela and Rivals, 2014). Later, corrected polished consensus sequences were aligned to the reference genome (Liu et al., 2017) using the genome mapping and alignment program (GMAP) (Wu and Watanabe, 2005). According to the genome mapping results, the polished consensus sequence was further corrected and clustered to remove redundancy and obtain the final high-quality isoforms using the TAPIS pipeline (Abdel-Ghany et al., 2016).
SUPPA (Alamancos et al., 2015) was used to identify seven types of alternative splicing (AS) events: skipped exon (SE), mutually exclusive exon (MX), alternative 5′ splice site (A5), alternative 3′ splice site (A3), retained intron (RI); alternative first exon (AF); and alternative last exon (AL). In order to identify poly(A) sequence signals and transcription factors (TFs) in our PacBio Iso-Seq data, MEME analysis (Zhang et al., 2015) and Animal Transcription Factor Database (AnimalTFDB) analysis were performed, respectively. In order to predict open reading frames (ORFs), the ANGEL pipeline (Kana et al., 2006) was used to find potential coding sequences for corrected isoforms. Meanwhile, PLEK (Lei et al., 2007), CNCI (Li et al., 2014), CPC (Sun et al., 2013) software, and Pfam (Finn et al., 2016) databases were used to predict the coding potential of PacBio sequencing data to identify lncRNA. At last, unmapped transcripts and novel gene transcript functions were annotated based on the following databases: NCBI non-redundant protein sequences (NR); NCBI non-redundant nucleotide sequences (NT); Pfam (protein family), KOG/COG (Clusters of Orthologous Groups of proteins (KOG/COG); Swiss-Prot (a manually annotated and reviewed protein sequence database); Kyoto Encyclopedia of Genes and Genomes (KEGG); and Gene Ontology (GO).
Assembly and Annotation of lncRNA and miRNA Transcriptomes
After quality control, paired-end clean reads were aligned to the reference genome using HISAT2 (Langmead and Salzberg, 2012). The mapped reads of each sample were assembled by StringTie (Mihaela et al., 2016). After transcripts that overlapped with known mRNAs or shorter than 200 bp were discarded, the coding potential analysis was performed on the remaining transcripts to define candidate lncRNA by CPC2 (Kang et al., 2017), CNCI (Li et al., 2014), and Pfam (Finn et al., 2016). Target gene prediction of lncRNA was then performed, and the expression level of both lncRNAs and coding transcripts were quantified.
The small RNA tags were mapped to the reference sequence by Bowtie (Langmead et al., 2009) without mismatch to analyze their expression and distribution on the reference. Mapped small RNA tags were used to look for known miRNAs based on the miRBase 20.0 database and miRDeep2 (Friedländer et al., 2012) and to predict novel miRNAs by miREvo (Wen et al., 2012) and miRDeep2. The target gene prediction of miRNA was then performed by miRanda (Enright et al., 2003), and miRNA expression levels were estimated.
Differential Expression and Functional Enrichment Analyses
A comparison of expression levels of mRNA, lncRNA, and miRNA between IB_I and IB_III was performed using DESeq or edgeR (Robinson et al., 2010). Differentially expressed mRNAs/lncRNAs/miRNAs were screened with a threshold q < 0.05. The lncRNA target genes were predicted by co-location analysis within 100 kb upstream and downstream of the lncRNA. And co-expression analysis was also used for lncRNA target prediction based on the Pearson correlation coefficient method with a threshold P < 0.05 and a correlation absolute value (R2) >0.95. Besides, GO and KEGG pathway enrichment analyses of differentially expressed genes or lncRNA target genes were implemented by the GO-seq R package (Young et al., 2010) and KOBAS software (Mao et al., 2005), respectively.
ceRNA Analysis
Pearson correlation coefficients of lncRNA and mRNA were calculated using the corresponding matrix data, and then a correlation test was performed. According to the ceRNA mechanism, significant positively correlated co-expression of lncRNA–mRNA was screened under conditions of correlation coefficient r > 0.95 and P < 0.001. The lncRNA–miRNA and mRNA–miRNA targeted pairs were constructed based on the targeting relationship between differential expressions of mRNA, lncRNA, and miRNA. The ceRNA regulatory network (lncRNA–miRNA–mRNA) was constructed based on the co-targeted lncRNA–miRNA–mRNA pairs and co-expressed lncRNA–mRNA pairs, followed by visualization of the ceRNA regulatory network using the Cytoscape software.
Validation of Transcriptome Sequencing
For validation of AS isoforms, five isoform genes were randomly selected, and gene-specific primers were designed using NCBI Primer-BLAST, and the expression profiles of nine lncRNAs and nine miRNAs were verified using qPCR technology. The PCR products of AS isoforms mixed with 10 × loading buffer and GelRed fluorescent nucleic acid dyes were electrophoresed in 1.4% agarose gel for 1 h. Real-time qPCR was performed to validate gene expression using TB Green® Fast qPCR Mix on the QuantStudio™ 6 Flex qRT-PCR system (ABI, Germany). The specific primer pairs of the validated genes were shown in Supplementary Table 1. The relative expression of the genes was calculated based on the comparative CT (2−ΔΔCT) method. Statistical analysis was performed using SPSS software. Data were presented as mean ± standard deviation (SD), with a statistical significance of P < 0.05.
Vector Construction
The whole lengths of lnc017705 and the dre-miR-24b-3p gene were cloned into the overexpression vector pcs2(+). The ZIP1 3′ UTR including the binding site of miR-24b-3p was amplified and inserted into the pmir-GLO vector (Promega, Madison, WI, USA) at the 3′ end of the luc2 (firefly luciferase) gene (pmir-GLO-ZIP1-3′ UTR). Similarly, the vectors of pmir-GLO-lnc017705 were obtained using the same method. Primer sequences were shown in Supplementary Table 1. All constructs were verified by sequencing.
Overexpression Assay
When the cell confluence reached about 80%, the pcs2+-miR-24b-3p and pcs2+-lnc017705 were transfected into cells, which were derived from the connective tissue surrounding IB in M. amblycephala and cultured in our laboratory using Lipofectamine 2000 (Invitrogen, USA). The relative expressions of miR-24b-3p, lnc017705, and ZIP1 were tested following the above method.
Luciferase Activity Assay
When the cell confluence reached about 80%, the pcs2+-miR-24b-3p, pcs2+-lnc017705, and pmir-GLO-ZIP1-3′ UTR were co-transfected into HEK293T cells. Similarly, the pcs2+-miR-24b-3p and pmir-GLO-lnc017705 were co-transfected into cells. After incubation for 24 h, the cells were collected and lysed, and dual luciferase activity was measured using a dual luciferase assay kit (Vazyme, Nanjing, China) and an automatic microplate reader (Molecular Devices, Sunnyvale, USA). The firefly luciferase activity was normalized against Renilla luciferase activity.
Results
Correction and Mapping of SMRT Sequencing
Of 326,088 generated CCS, FLNC reads from SMRT sequencing have a total number of 278,126 with a mean length of 2,543 bp (Table 1). Through the rigorous screening process, a total of 179,990 correct consensus reads with a mean length of 2,481 bp and N50 of 3,310 bp were generated. According to the mapping results, the sequences are classified into four types: unmapped (2.34%), multiple mapped (9.06%), reads mapped to a sense strand “+” (46.39%), reads mapped to an antisense strand “–” (42.21%) (Supplementary Figure 1A). The mapping coverage and identity of each transcript and density of FLNC sequences in chromosomes showed that transcripts with 98–100% coverage accounted for 75%, the mapping identification of which was 100% (Supplementary Figures 1B,C).
Table 1. Summary of reads from PacBio Iso-Seq.
Gene Structure Analysis and Identification of LncRNA From SMRT Sequencing
A total of 34,544 isoforms with a mean length of 2,537 bp and a N50 of 3,258 bp were obtained, including 2,209 isoforms of known genes, 26,302 novel isoform of known genes, and 6,033 isoforms of novel genes (Supplementary Figure 2A). The results also indicated that Iso-Seq generates longer transcripts and more gene isoforms than the reference genome (Supplementary Figures 2B–D). The density distribution of new genes and transcripts on chromosomes were analyzed and shown in Figure 2A. The SUPPA analysis revealed that SE is the most abundant AS event followed by A3, A5, and RI events (Figures 2A,B). In order to identify poly(A) sequence signals in our PacBio Iso-Seq data, a MEME analysis was performed on the sequence of 50 nucleotides upstream of the 8,541 poly(A) sites in 5,639 genes (Figures 2A,C). The conserved motif (AAUAAA) was identified upstream of poly(A) cleavage sites (Figure 2D). According to the Pfam file of the TF family, a total of 1,462 genes belonging to 59 TF families were identified using hmmsearch (Figure 2E). Besides, a total of 31,186 complete ORFs as well as 17,901 5′ and 27,602 3′ UTR sequences of the isoforms were predicted. And a total of 4,458 fusion transcripts were identified in the present study, whose linkage on intra-chromosomes and inter-chromosomes was shown in Figure 2A. After prediction of coding potential, transcripts without coding potential are defined as our candidate set of lncRNAs (Supplementary Figure 3A), which had a total number of 5,199 and were divided into four groups based on their position in the genome: long intergenic non-coding RNA (lincRNA), antisense, sense intronic, and sense overlapping (Supplementary Figure 3B). Besides, the results showed that the mean length of lncRNAs was longer than that of mRNAs, and most lncRNAs only had one exon (Supplementary Figures 3C,D).
Figure 2. Gene structure analysis. (A) Circos visualization of Iso-Seq data at the genome-wide level. From the outside to the inside: 1. Chromosome sequence. 2. Accumulation histogram of alternative splicing events (light blue indicates RI, green indicates A3, yellow indicates A5, purple indicates SE, red indicates MX, brown indicates AF, and dark blue indicates AL). 3. Alternative polyadenylation sites. 4. Distribution of novel transcripts from Iso-Seq. The color closer to red represents higher density, and the color closer to blue represents lower density. 5. Distribution of novel genes from Iso-Seq. The color closer to red represents higher density. 6. LncRNA density. The dot closer to the center represents lower density. 7. Linkage of fusion transcripts. Purple indicates intrachromosomal, and yellow indicates interchromosomal. (B) Statistics of alternative splicing events. (C) Statistics of alternative polyadenylation events. (D) MEME analysis of motif upstream of the poly(A) site. (E) Transcription factor analysis.
Expression Profiles of IB Growth-Dependent lncRNA and mRNA
In total, clean reads ranging from 83,654,272 to 91,891,080 and 87,038,458 to 95,445,674 were generated from IB_I and IB_III libraries, respectively, with three replicates. Of these, 84.59–87.58% reads in IB_I and 85.69–87.03% reads in IB_III can be mapped to the reference genome (Supplementary Table 2). The proportion of clean reads aligning uniquely to the reference sequence ranged from 78.83 to 81.92% in six libraries. According to the workflow as shown in Figure 3A, a total of 21,969 candidates lncRNAs were identified in IB libraries (Figure 3B), while the percentages of lincRNA, antisense lncRNAs, and intronic lncRNAs were 48.0, 7.0, and 45.0%, respectively (Figure 3C). The distribution pattern of differently expressed lncRNAs on different chromosomes was shown in Figure 3D. As illustrated in Figure 3E, the mean length of lncRNAs was shorter than that of protein-coding transcripts. In order to analyze the expression of different transcripts including lncRNAs, mRNAs, and transcripts of uncertain coding potential (TUCPs), StringTieeB software was used, and the result showed that the expression level of mRNA was higher than the expression level of lncRNA and TUCP (Figure 3F). Based on the expression level comparison, we found that 353 lncRNAs and 403 mRNAs were significantly differentially expressed between IB_I and IB_III (P < 0.05), respectively (Figures 4A,B; Supplementary Tables 3, 4). The hierarchical clustering analysis of their expression levels along with the Pearson correlation analysis was conducted and shown in Supplementary Figures 4A–C. Compared with those in IB_I, 170 and 183 lncRNAs were upregulated and downregulated in IB_III, respectively, while 219 and 184 mRNAs were upregulated and downregulated with an expression level fold change of more than 2 (Figures 4A,B). Besides, LNC_020948 (XLOC_153030) and MamblycephalaGene17723 (act2) were the most significantly differentially expressed lncRNA and mRNA between IB_III and IB_I, respectively. Among the lncRNAs, LNC_000059 (XLOC_007806) was the most highly expressed in IB_I, followed by LNC_020433 (XLOC_14908) and LNC_02798 (XLOC_159148), while the LNC_021798 (XLOC_159148) showed the highest expression in IB_III, followed by LNC_006232 (XLOC_052479) and LNC_004556 (XLOC_052479).
Figure 3. Identification of lncRNA in intermuscular bones. (A) Workflow for the preparation and analysis of lncRNA libraries. (B) Venn diagram analysis of lncRNAs identified from three software. (C) The classification of lncRNAs. (D) The distribution of lncRNAs in different chromosomes: the outermost circle is the chromosome; the second circle is the FPKM value of the samples on the corresponding chromosome; the third and fourth circles are the distribution of the significantly upregulated genes and the significantly downregulated genes in IB_III compared to IB_I on the chromosome, respectively. (E) Density distribution of lncRNA and mRNA with different ORF lengths. (F) Violin diagram of the expression of lncRNA, mRNA, and TUCP transcripts.
Figure 4. The analysis of differential expression and KEGG enrichment of lncRNA, mRNA, and miRNA in IB_I and IB_III. (A–C) Volcano map of differentially expressed lncRNA, mRNA, and miRNA. (D) The KEGG pathways enriched by lncRNA, mRNAand miRNA, relating to bone development. The color and size of shape represent p-value and enrichment numbers, respectively; the circle, triangle and square represent lncRNA, mRNAand miRNA, respectively. (E) The schematic diagram of MAPK/p53 signaling pathway. The black arrow represents the promotor; the red arrow represents the up-regulated (IB_III vs IB_I); the dotted line represents the intermediate process.
Expression Profiles of IB Growth-Dependent miRNA
After discarding low-quality reads, 3′ and 5′ adaptors and sequences with <18 nt, 19,553,745, 21,510,084, 25,513,637, 118,755,553, 22,297,421, and 19,168,240 clean reads were obtained from six libraries (Supplementary Table 5). Little differences in the length distribution of the sequences were identified among the six libraries. Most of the small RNAs were 21–23 nt in length, with 22 nt being the most abundant length (>40%) in all libraries. The total number of unique reads in IB_I and IB_III libraries ranged from 18,289,440 to 23,841,205 and 17,744,474 to 20,789,160, respectively. After reference sequence mapping, 16,999,090 (92.94%), 18,418,061 (92.53%), 21,712,011 (91.07%), 16,645,236 (93.81%), 19,471,908 (93.66%), and 16,694,103 (93.10%) mapped unique sRNAs were then used for miRNA identification and ncRNA annotation (Supplementary Table 5). The identification of the known miRNAs and novel miRNAs was accomplished by mapping the unique sequences in miRBase 20. The results showed that a total of 304 mapped unique miRNAs were detected and that 204 miRNAs were identified as novel in this study (Supplemental Table 6). After differential analysis, 126 miRNAs were identified with differential expression (P < 0.05) in IB_I and IB_III, including 68 upregulated and 58 downregulated miRNAs (Figure 4C). Of these differentially expressed miRNAs, the expression level of dre-miR-190b has the most significant difference and dre-miR-152 and dre-miR-99 had the highest expression levels in IB_I and IB_III, respectively (Supplementary Table 7).
GO and KEGG Analyses
We performed GO enrichment analysis to further investigate the potential functions of differentially expressed lncRNAs, miRNAs, and mRNAs in regulating IB development. As a result, GO enrichment of target genes of differentially expressed lncRNAs was categorized into 4,081 function groups (P < 0.05, Supplementary Table 8), while the GO enrichment of differentially expressed mRNAs was categorized into 1,520 function groups (P < 0.05, Supplementary Table 9). The GO function annotations of lncRNAs are shown in Supplementary Figure 5A. In addition, our data showed that the targets of differentially expressed lncRNAs and differentially expressed mRNAs were enriched into 158 and 92 pathways (P < 0.05, Supplementary Tables 10, 11), respectively. The top 20 enriched KEGG pathways of lncRNAs and mRNAs were presented in Supplementary Figures 6A,B. The miRNAs can regulate the expression of target genes; thus, differentially expressed target mRNA gene enrichment analysis can substitute for enrichment analysis of differentially expressed miRNA. A total of 3,517 GO term groups and 150 KEGG term groups of mRNAs (P < 0.05, Supplementary Tables 12, 13) targeted by differentially expressed miRNAs were enriched, while the enrichments and the top 20 KEGG pathways were shown in Supplementary Figures 5B,C, respectively. Based on the pathway analysis, the differentially expressed lncRNAs/miRNAs/mRNAs were commonly involved in eight signaling pathways related to bone development, including the Wnt signaling pathway, MAPK signaling pathway, p53 signaling pathway, mTOR signaling pathway, extracellular matrix (ECM)–receptor, ErbB signaling pathway, calcium signaling pathway, and Notch signaling pathway (Figure 4D). In these pathways, the differentially expressed genes (map2k6 and cytc) in the MAPK/p53 signaling pathway were significantly upregulated, and the differentially expressed genes (lama3 and thbs4b) in the ECM–receptor pathway were significantly downregulated. The pattern diagram of the MAPK/p53 signaling pathway was shown in Figure 4E.
Construction of the ceRNA Network
To further explore the potential lncRNA–miRNA–mRNA regulatory network, the co-expression analysis of differentially expressed lncRNAs and mRNAs between IB_I and IB_III was conducted firstly. A total of 3,197 co-expression lncRNA–mRNA pairs were obtained (r > 0.95 and P < 0.001) (Supplementary Table 14), and 2,053 lncRNA–miRNA–mRNA targeted pairs were generated (Supplementary Table 15). A total of 14 ceRNA pairs was obtained by merging the co-expression pairs and lncRNA–miRNA–mRNA targeted pairs, in which six upregulated lncRNAs and eight downregulated lncRNAs were involved in the ceRNA pairs. The detailed prediction process was shown in Figure 5, and the ceRNA regulatory network was visualized using Cytoscape software. Of these, we found that dre-miR-24b-3p and dre-miR-193b-3p were the core elements of the ceRNA pairs, interacting with four lncRNAs and three mRNAs, including LNC_007210, LNC_011298, LNC_001774, LNC_017705, MamblycephalaGene01083 (ZIP1), MamblycephalaGene20857 (C6), and MamblycephalaGene23275 (Figures 5B,C).
Figure 5. Identification of the ceRNA regulatory network. (A) The schematic diagram of lncRNA–miRNA–mRNA axes. “–” represents the negative correlation of the expression. “+” represents the positive correlation of the expression. (B) The filtering and categories of lncRNA–miRNA–mRNA interactions. (C) The cytoscape diagram of lncRNA–miRNA–mRNA interactions.
Validation of Transcriptome Sequencing
The gel electrophoresis separation results showed that all five genes have different amounts of isoforms (Supplementary Figure 2E), which indicated the reliability of the data generated in the present study and the different expression levels of different transcript isoforms. For example, in Gene15055, the expression of isoform3 is obviously higher than that of isoform1 and isoform2 (Supplementary Figure 2E). The qPCR results showed that the expression patterns between IB_I and IB_III in M. amblycephala were consistent with the RNA-Seq results (Figure 6). The qPCR and RNA-Seq results of nine lncRNAs and nine miRNAs in the two periods also showed similar expression changes (Figure 6).
Figure 6. The expression comparison of lncRNAs and miRNAs in RNA-Seq and qRT-PCR. The significant difference was presented as *P < 0.05, **P < 0.01, and ***P < 0.001.
Lnc017705 Acts as a ceRNA for miR-24b-3p
Given that ZIP1 can promote osteogenic differentiation by regulating runx2 and osterix expressions (Fu et al., 2018), we further explored the interaction mechanism of the regulatory axis lnc017705–miR-24b-3p–ZIP1 identified from the ceRNA network (Figures 5B,C). We transfected pcs2+-lnc017705 into connective tissue cells and found that lncRNA overexpression led to a 3-fold increase of lnc017705 expression (Figure 7A). The increasing lnc017705 could induce the expression of endogenous miR-24b-3p, and the pcs2+-miR-24b-3p transcription was promoted when co-transfecting pcs2+-lnc017705 and pcs2+-miR-24b-3p (Figure 7A). On the other hand, miR-24b-3p overexpression caused a significant increase in the expressions of lnc017705 and ZIP1 (Figure 7A). After co-transfection of pcs2+-lnc017705 and pcs2+-miR-24b-3p, the lncRNA expression increased compared with that only transfected by miR-24b-3p but lower than that only transfected by pcs2+-lnc017705 (Figure 7A). To determine the targeted relationship between lnc017705 and miR-24b-3p, we generated a miR-24b-3p sensor by inserting the lnc017705 sequence downstream of the luc2 gene in the pmir-GLO vector (Figure 7C). We found that miR-24b-3p markedly decreased the luc2 activity of pmir-GLO-lnc017705 in HEK293 cells (Figure 7D), which determined the target relationship between lnc017705 and miR-24b-3p.
Figure 7. The results of overexpression experiments and dual luciferase reporter gene of the lnc017705/miR-24b-3p/ZIP1 axis. (A) The result of lnc017705 and miR-24b-3p overexpression. (B) Targeting relationship prediction between miR-24b-3p and lnc017705 as well as ZIP1. (C) The construction of the dual luciferase vector. (D) The target verification of lnc017705 or ZIP1 with miRNA. (E) A regulatory schematic diagram of the lnc017705/miR-24a-3p-ZIP1 axis in osteogenic differentiation. *p < 0.05, **p < 0.01, ***p < 0.001.
Our bioinformatics analysis results show that ZIP1 is a potential target gene of miR-24b-3p (Figure 7B). Meanwhile, we found that ZIP1 mRNA expression could be rescued by co-transfection of pcs2+-lnc017705 and pcs2+-miR-24b-3p, compared with that only transfected by pcs2+-miR-24b-3p (Figure 7B). To eliminate the influence of miR-24b-3p on ZIP1, a dual luciferase activity assay was performed, referring to lnc017705 and miR-24b-3p. The results showed that miR-24b-3p could significantly increase luc2 activity when co-transfected with pcs2+-miR-24b-3p and pmir-GLO-ZIP1-3′ UTR, but luc2 activity was decreased after transfecting with additional pcs2+-lnc017705 (Figure 7D).
Discussion
The multi-omics analysis is an effective way to explore the molecular mechanism underlying complex trait development. In the present study, IB samples were isolated from M. amblycephala and used for PacBio library construction, which was then utilized to improve the draft genome annotation and full characterization of the transcriptome. Meanwhile, an integrative analysis of the transcriptome and degradome related to IB growth was presented. Among the KEGG pathways enriched by differentially expressed genes from the transcriptome and degradome of IB, eight KEGG pathways, such as the MAPK signaling pathway, p53 signaling pathway, Wnt signaling pathway, and ECM–receptor interaction pathway, were commonly enriched by differentially expressed lncRNAs, miRNAs, and mRNAs (Figure 4D). This is the first comprehensive overview of multi-omics about fish IB growth.
The MAPK–extracellular signal regulated kinase (ERK) pathway is an important link between the cell surface and nucleus to regulate proliferation and differentiation migration, as well as cell death (Lu et al., 2016). It also plays a critical role in bone formation (Schindeler and Little, 2006). MKK, identified as a mitogen-activated protein kinase (MAPK), can activate p38 protein through phosphorylating a critical -Thr-Gly-Tyr- motif (Goedert et al., 1997; Ono and Han, 2000). Then, phosphorylated p38 can act on the downstream target gene or participate in the p53 signaling pathway to promote cell apoptosis. In the present study, the gene map2k6 (MamblycephalaGene13655), belonging to MKK families, was upregulated in IB_III and might enhance the activity of p38, followed by promoting cell apoptosis. In the p53 signaling pathway, the electron carrier cytochrome c (cytc) was released under certain extreme situations (Kalpage et al., 2020). Once cytc is released from the mitochondria into the cytosol, it can interact with the protein apoptosis protease activating factor-1 (Apaf-1), which results in the formation of the apoptosome, leading to downstream caspase activation and cell death (Kalpage et al., 2020). Therefore, the significantly upregulated cytc (MamblycephalaGene13678) might be the reason for the slow growth of IBs in IB_III. Namely, the activity of cell apoptosis in IB_III was enhanced via upregulation of genes (map2k6 and cytc) involved in the MAPK/p53 signaling pathway.
The ECM consists of a complex mixture of structural and functional macromolecules, including glycosaminoglycans and fibrous proteins (e.g., collagen, elastin, fibronectin, and lammin) (Timpl, 1996; Van der Flier and Sonnenberg, 2001; Mariman and Wang, 2010), which are usually linked to the cytoskeleton through integrins to control various signaling pathways (Van der Flier and Sonnenberg, 2001). The previous studies have shown that the ECM proteins are associated with biomineralization in the bone tissue (Ravindran and George, 2014; Murshed, 2018) and provide structural and mechanical support for cells (Baroncelli et al., 2018). Therefore, downregulated ECM-related genes (e.g., MamblycephalaGene08032: lama3 and MamblycephalaGene22942: thbs4b) may explain the lower biomineralization activity in IB_III.
ceRNAs can regulate each other at the post-transcription level by competing for shared miRNAs (Qi et al., 2015). For example, lncRNA can indirectly regulate mRNA expression by absorbing endogenous miRNA, which has been introduced in various research fields (Wang R. et al., 2018; Zhang and Lu, 2018). In this study, we had constructed the ceRNA regulator network (lncRNA–miRNA–mRNA) to discover the core regulatory factors related to IB development. A total of 14 ceRNA pairs were generated, composed of 10 lncRNAs, 7 miRNAs, and 10 mRNAs. Among the ceRNA pairs, we found that dre-miR-24b-3p and dre-miR-193b-3p play the role of core elements, interacting with five lncRNAs and three mRNAs, such as lnc017705, MamblycephalaGene01083 (ZIP1), and MamblycephalaGene20857 (C6). As is known to all, bone homeostasis depends on bone resorption and formation by osteoclasts and osteoblasts (Chen et al., 2018), respectively, which are the two main cells participating in bone development (Matsuo and Irie, 2008). Meanwhile, osteoblast and osteoclast formation, differentiation, or apoptosis was affected by functional genes, such as runx2, osterix, rank, opg, β-catenin, and zip1 (Anderson et al., 1997; Komori et al., 1997; Nakashima et al., 2002; Boyce and Xing, 2008). Among them, ZIP1, as an essential trace element, is involved in diverse metabolic and signaling pathways (Mohammed et al., 2005) and promotes osteogenic differentiation by forming a zinc–Runx2/Osterix–ZIP1 regulation axis (Fu et al., 2018). The majority of Zn2+ in the body was involved in osteogenesis as an activator or coactivator of a variety of proteins such as runx2 and osterix with zinc finger motifs (Yamaguchi et al., 2008; Zhao et al., 2015). Besides, Mohammed et al. (2005) had determined that intracellular zinc content gradually accumulated and further found that ZIP1 was one dominant zinc transporter required for zinc uptake during osteogenesis of mesenchymal stem cells (MSCs). In addition, addition of Zn2+ or overexpression of ZIP1 can enhance MSC differentiation and promote the deposition of citrate and Ca2+ in mineralized MSCs (Fu et al., 2018). Therefore, it is reasonable to speculate that slow-growth IBs at the IB_III stage might result from decreased osteoblast differentiation and Ca2+ deposition caused by the downregulated ZIP1. Meanwhile, miR-24a-3p is the downstream targets of ZIP1 and lnc017705. Considering that miR-24 can regulate osteoblast differentiation of MSCs by targeting Tcf-1 (Nakashima et al., 2002), the lnc017705–miR-24a-3p–ZIP1 regulatory axis is likely to play an important role in the development of IB.
It is well-known that miRNAs inhibit the translation of mRNAs into protein and promote mRNA degradation. In the targeting relationship verification experiment, we found that lnc017705 can bind to miR-24b-3p to influence ZIP1 expression. However, when we singly transfected pcs2+-miR-24b-3p, the expressions of ZIP1 and lnc017705 were increased, which was contrary to expectation. This is because miR-24b-3p produced from pcs2+-miR-24b-3p could relieve the repression of target genes of endogenous miRNAs by available RNA-induced silencing complex (RISC) competition, resulting in upregulation of the corresponding mRNAs and proteins (Khan et al., 2009). This phenomenon was called machinery saturation (Castanotto et al., 2007). Besides, more miR-24b-3ps from pcs2+-miR-24b-3p were produced through co-overexpression of lnc017705 and miR-24b-3p, and ZIP1 expression will decrease since the machinery saturation was destroyed. Thereby, the mRNA level of ZIP1 and luc2 activity of pmir-GLO-ZIP1-3′ UTR vector were decreased by miRNA oversaturation state. Among numerous ceRNA researches, using miRNA mimics, but not pre-miRNAs (which rely on the nuclear export machinery), becomes the first choice because it can have highly effective inhibition and avoid saturation-related effects (McBride et al., 2008). Despite that, the higher off-targeted rate of miRNA mimics compared with miRNA expression vector is also a problem. So we had provided ideals by using pre-miRNA vector to investigate the ceRNA mechanism in our study. Finally, a summary of the lnc017705–miR-24a-3p–ZIP1 axis in osteogenesis difference is illustrated (Figure 7E), which offers novel clues to elucidate in depth the growth mechanism of teleost IB in future studies.
Conclusion
This is the first study to reveal the characteristics of the transcriptome and degradome associated with the growth of IBs in fish species. Through pathway enrichment analysis, we found that the expression change of key functional genes between IB_III and IB_I might be the reason for the slow growth of IBs in IB_III. In a word, the activity of cell apoptosis in IB_III was enhanced via upregulated genes (map2k6 and cytc) involved in the MAPK/p53 signaling pathway, and the lower mineralization activity in IB_III may be caused by downregulated genes (lama3 and thbs4b) in the ECM–receptor. Meanwhile, based on targeting relationship verification of lncRNA–miRNA–mRNA ceRNA and expression analyses, it is reasonable to believe that the slow IB growth in IB_III resulted from decreased osteoblast differentiation and Ca2+ deposition regulated by the lnc017705–miR-24a-3p–ZIP1 ceRNA network.
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 at: NCBI, PRJNA544738, and PRJNA640807.
Ethics Statement
All animals and experiments were conducted in accordance with the Guidelines for Experimental Animals of the Ministry of Science and Technology (Beijing, China). The study was approved by the Institutional Animal Care and Use Ethics Committee of Huazhong Agricultural University. All efforts were made to minimize suffering.
Author Contributions
Z-XG designed and supervised the study. SW, XD, JD, QLi, QLia, and G-YW performed the experiments. YC and SW analyzed the data and wrote the manuscript. All authors reviewed and approved the manuscript.
Funding
This work was financially supported by the National Natural Science Foundation of China (Grant No. 31872559), National Key Research and Development Program (Grant No. 2018YFD0900102), Modern Agriculture Industry Technology System Construction Projects of China titled as Staple Freshwater Fishes Industry Technology System (Grant No. CARS-46-08), Fundamental Research Funds for the Central Universities (Grant No. 2662018PY035), and Wuhan Applied Foundational Frontier Project (2020020601012253 and 2019020702011375).
Conflict of Interest
QLi and G-YW were employed by the company Wuhan Xianfeng Aquaculture Technology Co. Ltd, Wuhan, China.
The remaining 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.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.603815/full#supplementary-material
References
Abdel-Ghany, S. E., Hamilton, M., Jacobi, J. L., Ngam, P., Devitt, N., Schilkey, F., et al. (2016). A survey of the sorghum transcriptome using single-molecule long reads. Nat. Commun. 7:11706. doi: 10.1038/ncomms11706
Alamancos, G. P., Pagès, A., Trincado, J. L., Bellora, N., and Eyras, E. (2015). Leveraging transcript quantification for fast computation of alternative splicing profiles. RNA 21, 1521–1531. doi: 10.1261/rna.051557.115
Ali, D., and Soudeh, G. F. (2017). Long non-coding RNA expression intersecting cancer and spermatogenesis: a systematic review. Asian Pac. J. Cancer Prev. 18, 2601–2610. doi: 10.22034/APJCP.2017.18.10.2601
Anderson, D. M., Maraskovsky, E., Billingsley, W. L., Dougall, W. C., Tometsko, M. E., Roux, E. R., et al. (1997). A homologue of the TNF receptor and its ligand enhance T-cell growth and dendritic-cell function. Nature 390, 175–179. doi: 10.1038/36593
Baroncelli, M., Van der Eerden, B. C., Kan, Y. Y., Alves, R. D., Demmers, J. A., Van de Peppel, J., et al. (2018). Comparative proteomic profiling of human osteoblast-derived extracellular matrices identifies proteins involved in mesenchymal stromal cell osteogenic differentiation and mineralization. J. Cell Physiol. 233, 387–395. doi: 10.1002/jcp.25898
Boyce, B. F., and Xing, L. (2008). Functions of RANKL/RANK/OPG in bone modeling and remodeling. Arch. Biochem. Biophys. 473, 139–146. doi: 10.1016/j.abb.2008.03.018
Cai, J., Li, L., Song, L., Xie, L., Luo, F., Sun, S., et al. (2019). Effects of long term antiprogestine mifepristone (RU486) exposure on sexually dimorphic lncRNA expression and gonadal masculinization in Nile tilapia (Oreochromis niloticus). Aquat. Toxicol. 215:105289. doi: 10.1016/j.aquatox.2019.105289
Cao, D. C., Kuang, Y. Y., Zheng, X. H., Tong, G. X., Li, C. T., and Sun, X. W. (2015). Comparative analysis of intermuscular bones in three strains of common carp. J. Appl. Ichthyol. 31, 32–36. doi: 10.1111/jai.12483
Castanotto, D., Sakurai, K., Lingeman, R., Li, H., Shively, L., Aagaard, L., et al. (2007). Combinatorial delivery of small interfering RNAs reduces RNAi efficacy by selective incorporation into RISC. Nucleic Acids. Res. 35, 5154–5164. doi: 10.1093/nar/gkm543
Chen, X., Wang, Z., Duan, N., Zhu, G., Schwarz, E. M., and Xie, C. (2018). Osteoblast-osteoclast interactions. Connect Tissue Res. 59, 99–107. doi: 10.1080/03008207.2017.1290085
Chen, Z., Omori, Y., Koren, S., Shirokiya, T., Kuroda, T., Miyamoto, A., et al. (2017). De novo assembly of the goldfish (Carassius auratus) genome and the evolution of genes after whole-genome duplication. Sci. Adv. 5:eaav0547. doi: 10.1126/sciadv.aav0547
Danos, N., and Ward, A. B. (2012). The homology and origins of intermuscular bones in fishes: phylogenetic or biomechanical determinants. Biol. J. Linn. Soc. 106, 607–622. doi: 10.1111/j.1095-8312.2012.01893.x
Dong, Z., Zhong, H., Li, L., Yuan, X., Miao, W., Chen, Q., et al. (2006). Preliminary study on intermuscular bones of several cultured cyprinids. J. Ocean U. China 15, 425–429. doi: 10.1360/aps050066
Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA targets in Drosophila. Genome Biol. 5:R1. doi: 10.1186/gb-2003-5-1-r1
Finn, R. D., Coggill, P., Eberhardt, R. Y., Eddy, S. R., Mistry, J., Mitchell, A. L., et al. (2016). The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 44, D279–285. doi: 10.1093/nar/gkv1344
Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2012). MiRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi: 10.1093/nar/gkr688
Fu, X., Li, Y., Huang, T., Yu, Z., Ma, K., Yang, M., et al. (2018). Runx2/osterix and zinc uptake synergize to orchestrate osteogenic differentiation and citrate containing bone apatite formation. Adv. Sci. 5:1700755. doi: 10.1002/advs.201700755
Gan, S., Huang, Z., Liu, N., Su, R., Xie, G., Zhong, B., et al. (2016). MicroRNA-140-5p impairs zebrafish embryonic bone development via targeting BMP-2. FEBS Lett. 590, 1438–1446. doi: 10.1002/1873-3468.12190
Gao, S., Ren, Y., Sun, Y., Wu, Z., Ruan, J., He, B., et al. (2016). PacBio full-length transcriptome profiling of insect mitochondrial gene expression. RNA Biol. 13, 820–825. doi: 10.1080/15476286.2016.1197481
Goedert, M., Cuenda, A., Craxton, M., Jakes, R., and Cohen, P. (1997). Activation of the novel stress-activated protein kinase SAPK4 by cytokines and cellular stresses is mediated by SKK3 (MKK6); comparison of its substrate specificity with that of other SAP kinases. EMBO J. 16, 3563–3571. doi: 10.1093/emboj/16.12.3563
Gong, G., Dan, C., Xiao, S., Guo, W., Huang, P., Xiong, Y., et al. (2018). Chromosomal-level assembly of yellow catfish genome using third-generation DNA sequencing and Hi-C analysis. Gigascience 7:giy120. doi: 10.1093/gigascience/giy120
Hobert, O. (2008). Gene regulation by transcription factors and microRNAs. Science 319, 1785–1786. doi: 10.1126/science.1151651
Jiang, W., Jia, Y., Liu, S., Li, Q., Li, T., and Gu, Z. (2016). Comparative analysis of intermuscular bones in hybrid F1, F2 of (C. alburnus) (♀) × (M. amblycephala) (♂). Acta. Hydrobiol. Sin. 40, 277–286. doi: 10.7541/2016.38
Kalpage, H. A., Wan, J., Morse, P. T., Zurek, M. P., Turner, A. A., Khobeir, A., et al. (2020). Cytochrome c phosphorylation: control of mitochondrial electron transport chain flux and apoptosis. Int. J. Biochem. Cell Biol. 121:105704. doi: 10.1016/j.biocel.2020.105704
Kana, S., Jun, A., and Yoichi, M. (2006). ANGLE: a sequencing errors resistant program for predicting protein coding regions in unfinished cDNA. J. Bioinf. Comput. Biol. 4, 649–664. doi: 10.1142/S0219720006002260
Kang, Y., Yang, D., Kong, L., Hou, M., Meng, Y., Wei, L., et al. (2017). CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 45, W12–W16. doi: 10.1093/nar/gkx428
Khan, A. A., Betel, D., Miller, M. L., Sander, C., Leslie, C. S., and Marks, D. S. (2009), Transfection of small RNAs globally perturbs gene regulation by endogenous microRNAs. Nat. Biotechnol. 27, 549–555. doi: 10.1038/nbt.1543
Komori, T., Yagi, H., Nomura, S., Yamaguchi, A., Sasaki, K., Deguchi, K., et al. (1997). Targeted disruption of Cbfa1 results in a complete lack of bone formation owing to maturational arrest of osteoblasts. Cell 89, 755–764. doi: 10.1016/S0092-8674(00)80258-5
Lan, T., Chen, Y. L., Gul, Y., Zhao, B. W., and Gao, Z. X. (2019). Comparative expression analysis of let-7 microRNAs during ovary development in Megalobrama Amblycephala. Fish Physiol. Biochem. 45, 1101–1115. doi: 10.1007/s10695-019-00624-7
Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25
Larsson, L., Castilho, R. M., and Giannobile, W. V. (2015). Epigenetics and its role in periodontal diseases: a state-of-the-art review. J. Periodontol. 86, 556–568. doi: 10.1902/jop.2014.140559
Lei, K., Yong, Z., Zhi-Qiang, Y., Xiao-Qiao, L., Shu-Qi, Z., Liping, W., et al. (2007). CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–349. doi: 10.1093/nar/gkm391
Li, A., Zhang, J., and Zhou, Z. (2014). PLEK: a tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics 15:311. doi: 10.1186/1471-2105-15-311
Li, L., Zhong, Z., Zeng, M., Liu, S., Zhou, Y., Xiao, J., et al. (2013). Comparative analysis of intermuscular bones in fish of different ploidies. Sci. China Life Sci. 56, 341–350. doi: 10.1007/s11427-013-4465-5
Liu, H., Chen, C., Gao, Z., Min, J., Gu, Y., Jian, J., et al. (2017). The draft genome of blunt snout bream (Megalobrama amblycephala) reveals the development of intermuscular bone and adaptation to herbivorous diet. Gigascience 6, 1–13. doi: 10.1093/gigascience/gix039
Liu, Y., Zhang, R., and Ying, K. (2015). Long non-coding RNAs: novel links in respiratory diseases (review). Mol. Med. Rep. 11, 4025–4031. doi: 10.3892/mmr.2015.3290
Lu, D., Qu, J., Lei, J., Wang, P., Sun, C., Huang, L., et al. (2016). LPS-stimulated inflammation inhibits BMP-9-induced osteoblastic differentiation through crosstalk between BMP/MAPK and Smad signaling. Exp. Cell Res. 341, 54–60. doi: 10.1016/j.yexcr.2016.01.009
Luo, M., Wang, L., Yin, H., Zhu, W., Fu, J., and Dong, Z. (2019). Integrated analysis of long non-coding RNA and mRNA expression in different colored skin of koi carp. BMC genomics 20:515. doi: 10.1186/s12864-019-5894-8
Mao, X., Cai, T., Olyarchuk, J. G., and Wei, L. (2005). Automated genome annotation and pathway identification using the KEGG Orthology (KO) as a controlled vocabulary. Bioinformatics 21, 3787–3793. doi: 10.1093/bioinformatics/bti430
Mariman, E. C., and Wang, P. (2010). Adipocyte extracellular matrix composition, dynamics and role in obesity. Cell Mol. Life Sci. 67, 1277–1292. doi: 10.1007/s00018-010-0263-4
Matsuo, K., and Irie, N. (2008). Osteoclast-osteoblast communication. Arch. Biochem. Biophys. 473, 201–209. doi: 10.1016/j.abb.2008.03.027
McBride, J. L., Boudreau, R. L., Harper, S. Q., Staber, P. D., Monteys, A. M., Martins, I., et al. (2008). Artificial miRNAs mitigate shRNA-mediated toxicity in the brain: implications for the therapeutic development of RNAi. Proc. Natl. Acad. Sci. U.S.A. 105 5868–5873. doi: 10.1073/pnas.0801775105
Mihaela, P., Daehwan, K., Geo, M. P., Jeffrey, T. L., and Steven, L. S. (2016). Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11, 1650–1667. doi: 10.1038/nprot.2016.095
Mohammed, A. K., Surasri, N. S., Guang, B., Sunia, A., and Anandarup, G. (2005). Expression of the zinc transporter ZIP1 in osteoclasts. Bone 37, 296–304. doi: 10.1016/j.bone.2005.04.035
Murshed, M. (2018). Mechanism of bone mineralization. Cold Spring Harb. Perspect Med. 8:a031229. doi: 10.1101/cshperspect.a031229
Nakashima, K., Zhou, X., Kunkel, G., Zhang, Z., Deng, J. M., Behringer, R. R., et al. (2002). The novel zinc finger-containing transcription factor osterix is required for osteoblast differentiation and bone formation. Cell 108, 17–29. doi: 10.1016/S0092-8674(01)00622-5
Ono, K., and Han, J. (2000). The p38 signal transduction pathway: activation and function. Cell Signal. 12, 1–13. doi: 10.1016/S0898-6568(99)00071-6
Perazza, C. A., de Menezes, J. T. B., Ferraz, J. B. S., Pinaffi, F. L. V., Silva, L. A., and Hilsdorf, A. W. S. (2016). Lack of intermuscular bones in specimens of Colossoma macropomum: An unusual phenotype to be incorporated into genetic improvement programs. Aquaculture 472, 57–60. doi: 10.1016/j.aquaculture.2016.05.014
Qi, X., Zhang, D. H., Wu, N., Xiao, J. H., Wang, X., and Ma, W. (2015). CeRNA in cancer: possible functions and clinical implications. J. Med. Genet. 52, 710–718. doi: 10.1136/jmedgenet-2015-103334
Ravindran, S., and George, A. (2014). Multifunctional ECM proteins in bone and teeth. Exp. Cell Res. 325, 148–154. doi: 10.1016/j.yexcr.2014.01.018
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). EdgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Salmela, L., and Rivals, E. (2014). LoRDEC: accurate and efficient long read error correction. Bioinformatics 30, 3506–3514. doi: 10.1093/bioinformatics/btu538
Schindeler, A., and Little, D. G. (2006). Ras-MAPK signaling in osteogenic differentiation: Friend or foe? J. Cell Physiol. 21, 1331–1338. doi: 10.1359/jbmr.060603
Sharon, D., Tilgner, H., Grubert, F., and Snyder, M. (2013). A single-molecule long-read survey of the human transcriptome. Nat. Biotechnol. 31, 1009–1014. doi: 10.1038/nbt.2705
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41:e166. doi: 10.1093/nar/gkt646
Timpl, R. (1996). Macromolecular organization of basement membranes. Curr. Opin. Cell Biol. 8, 618–624. doi: 10.1016/S0955-0674(96)80102-5
Van der Flier, A., and Sonnenberg, A. (2001). Function and interactions of integrins. Cell Tissue Res. 305, 285–298. doi: 10.1007/s004410100417
Wan, S., Yi, S., Zhong, J., Wang, W., Jiang, E., Chen, B., et al. (2014). Development and morphological observation of intermuscular bones in megalobrama amblycephala. Acta. Hydrobiol. Sin. 38, 1143–1151. doi: 10.7541/2014.166
Wan, S. M., Yi, S. K., Zhong, J., Nie, C. H., Guan, N. N., and Chen, B. X., (2015). Identification of microRNA for intermuscular bone development in blunt snout bream (Megalobrama amblycephala). Int. J. Mol. Sci. 16, 10686–10703. doi: 10.3390/ijms160510686
Wan, S. M., Yi, S. K., Zhong, J., Nie, C. H., Guan, N. N., Zhang, W. Z., et al. (2016). Dynamic mRNA and miRNA expression analysis in response to intermuscular bone development of blunt snout bream (Megalobrama amblycephala). Sci. Rep. 6:31050. doi: 10.1038/srep31050
Wang, M., Jiang, S., Wu, W., Yu, F., Chang, W., Li, P., et al. (2018). Non-coding RNAs function as immune regulators in teleost fish. Front. Immuno. 9:2801. doi: 10.3389/fimmu.2018.02801
Wang, R., Zhang, S., Chen, X., Li, N., Li, J., Jia, R., et al. (2018). CircNT5E acts as a sponge of miR-422a to promote glioblastoma tumorigenesis. Cancer Res. 78, 4812–4825. doi: 10.1158/0008-5472.CAN-18-0532
Wen, M., Shen, Y., Shi, S., and Tang, T. (2012). MiREvo: an integrative microRNA evolutionary analysis platform for next-generation sequencing experiments. BMC Bioinformatics 13:140. doi: 10.1186/1471-2105-13-140
Wu, T. D., and Watanabe, C. K. (2005). GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics 21, 1859–1875. doi: 10.1093/bioinformatics/bti310
Xiong, X. M., Robinson, N. A., Zhou, J. J., Chen, Y. L., Wang, W., Wang, X. B., et al. (2018). Genetic parameter estimates for intermuscular bone in blunt snout bream (Megalobrama amblycephala) based on a microsatellite-based pedigree. Aquaculture 502, 371–377. doi: 10.1016/j.aquaculture.2018.12.072
Xu, X., Zheng, J., Qian, Y., and Luo, C. (2015). Normally grown and developed intermuscular bone-deficient mutant in grass carp, Ctenopharyngodon Idellus. Chinese Sci. Bull. 60, 52–57. doi: 10.1360/N972014-00637
Yamaguchi, M., Goto, M., Uchiyama, S., and Nakagawa, T. (2008). Effect of zinc on gene expression in osteoblastic MC3T3-E1 cells: enhancement of Runx2, OPG, and regucalcin mRNA expressions. Mol. Cell Biochem. 312, 157–166. doi: 10.1007/s11010-008-9731-7
Yang, G., Lu, X., and Yuan, L. (2014). LncRNA: a link between RNA and cancer. BBA-Gen. Subjects 1839, 1097–1109. doi: 10.1016/j.bbagrm.2014.08.012
Young, M. D., Wakefield, M. J., Smyth, G. K., and Oshlack, A. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11:R14. doi: 10.1186/gb-2010-11-2-r14
Yuan, W., Jiang, S., Sun, D., Wu, Z., Wei, C., Dai, C., et al. (2019). Transcriptome profiling analysis of sex-based differentially expressed mRNAs and lncRNAs in the brains of mature zebrafish (Danio rerio). BMC genomics 20:830. doi: 10.1186/s12864-019-6197-9
Zhang, G., Yin, S., Mao, J., Liang, F., Zhao, C., Li, P., et al. (2016). Integrated analysis of mRNA-seq and miRNA-seq in the liver of Pelteobagrus vachelli in response to hypoxia. Sci. Rep. 6:22907. doi: 10.1038/srep22907
Zhang, H., and Lu, W. (2018). LncRNA SNHG12 regulates gastric cancer progression by acting as a molecular sponge of miR-320. Mol. Med. Rep. 17, 4812–4825. doi: 10.3892/mmr.2017.8143
Zhang, H. M., Liu, T., Liu, C. J., Song, S., Zhang, X., Liu, W., et al. (2015). AnimalTFDB 2.0: a resource for expression, prediction and functional study of animal transcription factors. Nucleic Acids Res. 43, D76–81. doi: 10.1093/nar/gku887
Zhao, B. W., Zhou, L. F., Liu, Y. L., Wan, S. M., and Gao, Z. X. (2017). Evolution of fish let-7 microRNAs and their expression correlated to growth development in blunt snout bream. Int. J. Mol. Sci. 18:646. doi: 10.3390/ijms18030646
Zhao, W., Wu, C., Dong, Y., Ma, Y., Jin, Y., and Ji, Y. (2015). MicroRNA-24 regulates osteogenic differentiation via targeting T-Cell factor-1, Int. J. Mol. Sci. 16, 11699–11712. doi: 10.3390/ijms160511699
Keywords: teleost, intermuscular bones, growth, transcriptome profiling, RNA interaction
Citation: Chen Y, Wan S, Li Q, Dong X, Diao J, Liao Q, Wang G-Y and Gao Z-X (2021) Genome-Wide Integrated Analysis Revealed Functions of lncRNA–miRNA–mRNA Interaction in Growth of Intermuscular Bones in Megalobrama amblycephala. Front. Cell Dev. Biol. 8:603815. doi: 10.3389/fcell.2020.603815
Received: 08 September 2020; Accepted: 23 December 2020;
Published: 04 February 2021.
Edited by:
Kang-le Lu, Jimei University, ChinaReviewed by:
Wanbo Li, Jimei University, ChinaChangwei Shao, Yellow Sea Fisheries Research Institute (CAFS), China
Qiang Lin, Chinese Academy of Sciences, China
Copyright © 2021 Chen, Wan, Li, Dong, Diao, Liao, Wang and Gao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ze-Xia Gao, gaozexia@hotmail.com; gaozx@mail.hzau.edu.cn
†These authors have contributed equally to this work