Screening and identification of miRNAs regulating Tbx4/5 genes of Pampus argenteus

Sample preparation, sRNA isolation, and cDNA library construction

The samples of P. argenteus used for the experiments were farmed from May to July 2021 in Xixuan Island Breeding Base, Zhejiang Marine Fisheries Research Institute, Zhoushan, Zhejiang Province (28°31′56′′ N, 122°12′24′′ E). The fishes were cultured in indoor cement pools of 25 m² area, and a water depth of 1.4 m. The transcriptomes of P. argenteus were sequenced at days 1, 7, and 13 of larval growth, and two biological replicates were sequenced for each stage. The data represented the mean of two biological replicates for the three individual transcriptomes, and were labelled O_D (A, B), S_D (A, B), and T_D (A, B). The total RNA was extracted in 1.5 ml centrifuge tubes (Guangzhou Jet Bio-Filtration Co., Ltd., Guangzhou, China) by the Trizol method. RNA degradation and contamination were monitored using 1% agarose gels, and the purity, concentration, and integrity of the extracted RNA were evaluated. For each sample, 3 μg of total RNA was used as the input material for constructing the sRNA library. The sRNA library was prepared and sequenced using an Illumina HiSeq 2000 NGS platform. The quality of the library was finally evaluated using DNA High Sensitivity Chips on an Agilent Bioanalyzer 2100 system. The protocol for animal experiments was approved by the Animal Health and Use Committee of Ningbo University (permit number: NO20210706).

Sequence data analyses and qRT-PCR analysis of differentially expressed (DE) miRNAs

The clean reads were obtained by processing the raw reads in FASTQ format using customized Perl and Python scripts. Then, clean reads of certain range of lengths were selected for subsequent analyses. The sRNA tags were mapped to the reference sequence (P. argenteus genomic sequence, GenBank ID: JHEK00000000.1) using Bowtie without any mismatches, and their expression and distribution were analyzed (Langmead et al., 2009). The Rfam10.1 and RepeatMasker databases were used for analyzing the snRNAs, snoRNAs, tRNAs, rRNAs, and repetitive sequences. The miRDeep2 software was used for identifying the known miRNAs and predicting the novel miRNAs of P. argenteus (Friedländer et al., 2012). The dicer cleavage site, secondary structure, and minimum free energy of the miRNAs were determined using the MIREAP software, following which the sRNAs were annotated and summarized.

The transcriptome of the larvae at different developmental stages was compared for identifying the significantly DE miRNAs (DEMs; |log₂(fold change)| > 1, P ≤ 0.05). The TargetScan and MiRanda software were used for analyzing the differences in gene expression among the different development stages. The predicted target genes were subsequently subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses (Mao et al., 2005). In order to verify the sequencing data, eight significantly DE miRNAs were analyzed by qRT-PCR, and the primers used for the specific miRNAs are enlisted in Table 1. U6 was used as the internal control. The qRT-PCR experiments were performed in triplicate, at least three times for each miRNA. The alterations in miRNA expression were determined using the 2^−ΔΔCt method.

Screening of miRNAs that regulate Tbx4/5 genes, and cloning of Tbx4/5

The miRNAs that regulate Tbx4/5 gene expression were screened using the PITA and miRanda software. A Wayne diagram was constructed using the screening results, and the miRNAs in the intersection area were further selected for identifying the miRNAs that regulate the expression of Tbx4/5 genes. The total RNA was extracted by the Trizol method, and subsequently reverse transcribed Primers containing the NotI and XhoI restriction sites, and the GAAT and CCG protective bases, were designed based on the obtained full-length sequence of the Tbx4/5 genes of P. argenteus, and the 3’-UTR fragments of the Tbx4/5 genes were cloned. Site-directed mutagenesis was performed at the corresponding positions based on the 3′-UTR fragments and corresponding miRNA binding sites, using a Mut Express MultiS Fast Mutagenesis Kit V2 (Vazyme Biotech Co., Ltd, Nanjing, China). The primers used for site-directed mutagenesis were designed using the miRNA Design software, and the primers included 4Mut-102n-F/R, 4Mut-301cd-F/R, 4Mut-19bd-F/R, 4Mut-301bd-F/R, 4Mut-113n-F/R, 5Mut-187d-F/R, 5Mut-201n-F/R, 5Mut-219d-F/R, and 5Mut-460d-F/R (Tables 2 and 3). PCR amplification was subsequently performed, and the target fragment was recovered using a Gel Extraction Kit. The pmirGLO plasmid was digested using the Nhel and Xbal restriction enzymes, and the target fragment was subsequently recovered. Following ligation, transformation, bacterial selection, and verification of the bacterial liquid by PCR, the positive bacterial liquid was sequenced.

Construction of dual luciferase reporter gene vector and mutant plasmids

The pMD19-T-TBX4/5-3′-UTR and pSicheck-2 plasmids were digested by XhoI and NotI. and the desired fragments were purified to obtain the 3′-UTR fragments of the Tbx4/5 genes and the linearized pSicheck-2 plasmid containing the digestion sites. The 3′-UTR fragments of the Tbx4/5 genes were ligated into the pSicheck-2 vector, and the resulting recombinant reporter vector was denoted as the pSicheck-3′-UTR-WT plasmid. The reporter vector was transformed into DH5α competent cells (TaKaRa, Zhejiang, PR China), and the positive clones were screened in LB medium containing ampicillin. The plasmids were extracted and verified by enzymatic digestion. Using the psiCHECK-3′-UTR-WT plasmid as the template, site-directed mutagenesis was performed at the binding site of the 3′-UTR domains of the corresponding miRNAs. The protocol was performed according to the instructions in the Mut Express MultiS Fast Mutagenesis Kit V2. The mutant plasmids obtained by site-directed mutagenesis were denoted as 3′-UTR-4Mut-102, 3′-UTR-4Mut-301c, 3′-UTR-4Mut-19b, 3′-UTR-4Mut-301b, 3′-UTR-4Mut-113n, 3′-UTR-4Mut-589n; 3′-UTR-5Mut-187d, 3′-UTR-5Mut-201n, 3′-UTR-5Mut-219, and 3′-UTR-5Mut-460d. The plasmids were extracted using Endo Free Plasmid Mini Kit I (Omega Bio-Tek, Norcross, GA, USA). The experimental protocol was performed according to the instructions provided with the kit, and the mutant plasmids were subsequently sequenced.

Luciferase reporter system

High-confidence DE miRNA-mRNA pairs were selected and targeted for validation. The luciferase activity test confirmed the potential relationship between miRNAs and their target genes. Specifically, the 3′-UTR sequences were synthesized, which contained the miRNA binding sites of the candidate target genes, and subsequently inserted into the multiple cloning sites of the psiCHECK-3′-UTR-WT plasmid. The recombinant psiCHECK-3′-UTR-WT vector and the miRNA mimics or negative control (NC) were co-transfected into DH5α competent cells (TaKaRa, Zhejiang, PR China) at a confluence of 70% using Lipofectamine 3000, according to the manufacturer’s instructions. After 48 h of transfection, the activity of the dual luciferase was detecting using a Dual-Glo^® Luciferase Assay Kit (Promega, Madison, WI, USA).

Characteristics of the sequencing data

Six sRNA libraries were constructed for the three different developmental stages. The raw transcriptome data of P. argenteus were generated using an Illumina sequencing platform, and submitted to the NCBI GenBank database (accession numbers: SAMN20702608–SAMN20702613). The Q20 and Q30 values of all the samples were approximately 98% and 96%, respectively, and the GC content was approximately 48–50%. The length of the 6 sRNA libraries ranged from 11,431,807 to 16,854,046, and the proportion of clean reads to raw reads was above 96%, which indicated that the quality of the sequencing data was highly reliable (Table 4). The lengths of the sRNAs, miRNAs, and siRNAs were 18–35 nt, 21–22 nt, and 24 nt, respectively. Based on the distribution of the sequence lengths of the clean reads, the most abundant size class ranged from 21 to 23 nt (Fig. S1). Additionally, the bias of the first base of known miRNAs was statistically analyzed. The results demonstrated that there was a strong preference for U in the first base of miRNAs in the three different larval stages of P. argenteus, and there was an equal proportion of A and C in the first base of novel miRNAs, without any specific preference (Figs. S2 and S3). Figure 1A indicates that the Pearson correlation coefficients among the groups were greater than 0.8, which indicated a strong correlation among the groups, and also implied that the experimental data were highly reliable.

Identification of DE miRNAs

In this study, we analyzed the expression levels of known and novel miRNAs on days 1, 7, and 13 of the larval stage for determining the significant differences among the 6 sRNA libraries. A total of 662 miRNAs were detected, of which 257 known and 405 novel miRNAs were identified in the 6 sRNA libraries. Compared to day 1, 182 miRNAs were DE on day 7, of which 77 miRNAs were downregulated and 105 miRNAs were upregulated, while 278 miRNAs were DE on day 13, of which 136 miRNAs were downregulated and 142 miRNAs were upregulated. Compared to day 13, 4 miRNAs were DE on day 7, of which 3 miRNAs were downregulated and 1 miRNA was upregulated (Figs. 1B–1E; Table 5). These results indicated that numerous miRNAs are involved in the developmental regulation of P. argenteus.

The results of hierarchical clustering of the miRNAs revealed that the miRNAs had different expression levels among the three different stages of development (Fig. 1F). The genes that were DE between the day 1 (O_D) and day 7 (S_D) larval stages, and between the day 1 (O_D) and day 13 (T_D) larval stages were inversely expressed, indicating that the larvae were in the peak stage of differentiation. However, the number of genes DE between days 7 (S_D) and 13 (T_D) of the larval stages was relatively small, suggesting the initiation of development. Of these, the expression levels of 4 miRNAs (dre-let-7f, dre-miR-129-1-3p, dre-miR-129-3-3p, and novel_154) differed at the three larval developmental periods (Table 6). Compared with those of the O_D larval stage, the expression levels of dre-let-7f, dre-miR-129-1-3p, and dre-miR-129-3-3p miRNAs were upregulated, while those of novel_154 were downregulated in the S_D and T_D larval stages. Compared with those of the T_D larval stages, the expression levels of dre-let-7f were upregulated, while those of dre-miR-129-1-3p, dre-miR-129-3-3p, and novel_154 were downregulated in the O_D and S_D larval stages.

Prediction and functional analysis of target genes

The potential target genes of the DE miRNAs were predicted using the TargetScan webserver. The predicted target genes were classified into three major functional ontologies based on the GO classification system, namely biological process, cellular component, and molecular function. In the biological processes category, the target genes were most enriched in the cellular processes term, followed by the single-cell processes and metabolic processes terms. Additionally, the target genes were most enriched in the cell components category. In the molecular function category, the genes were most enriched in binding and catalytic activity than other terms. The KEGG database (http://www.kegg.jp/kegg/pathway.html) was subsequently used for pathway analysis of the target genes. The 9,440 DE genes in the transcripts of day 1 and day 7 larvae of P. argenteus were subjected to pathway enrichment analyses, which revealed that the main enriched pathways were the MAPK signaling pathway, calcium signaling pathway, tight junction, and the Wnt signaling pathway. Pathway enrichment analyses of the 14,236 DE genes in the transcripts of day 1 and day 13 larvae revealed that the main enriched pathways were the MAPK signaling pathway, calcium signaling pathway, Wnt signaling pathway, adrenergic signaling in cardiomyocytes, and protein processing in endoplasmic reticulum. Similarly, the main enriched pathways of the 45 DE genes in the transcripts of day 7 and day 13 larvae were spliceosome and adrenergic signaling in cardiomyocytes (Fig. 2).

Expression levels of DE miRNAs

In order to verify the reliability of the miRNA-Seq data, the expression levels of the DE miRNAs that regulated the target genes during the development of P. argenteus were detected by qRT-PCR. As depicted in Fig. 3, the expression levels of dre-miR-19b-5p, dre-miR-219-3p, dre-miR-460-5p, novel-212, novel-224, and novel-589 were downregulated. However, graphical depiction of the expression levels of dre-miR-301b-5p and novel-113 showed an “inverted V” curve, and demonstrated that the expression of these miRNAs peaked on day 7. The expression of these miRNAs was consistent with the patterns determined by high-throughput sequencing. The findings confirmed that the miRNA sequencing data were reliable and could be used in further analyses.

Screening of miRNAs regulating Tbx4/5 genes

Similar to the results of previous studies, the expression levels of Tbx4 mRNA increased continually from day 1 to day 13 in this study, and was highest in day 13 larvae, while the expression of Tbx5 mRNA showed a downward trend and was highest when the primordia of the 1-day-old pectoral fins appeared. Of the different tissues, the positive signal of the Tbx4 gene was primarily concentrated in the abdominal epithelium of P. argenteus; however, the expression levels of the Tbx4 protein were low in the abdominal epithelium of day 7 larvae and adult fishes. The positive signal of the Tbx5 gene was primarily concentrated in the pectoral fins, while the expression of the Tbx5 protein was relatively high in the pectoral fins of day 1 larvae and adult fishes (Zhu, 2020). We therefore speculated that the Tbx4/5 genes of P. argenteus could be regulated by miRNAs. In this study, we used the miRanda bioinformatics tool (version 2.2a) for analyzing and predicting the miRNAs that regulate the Tbx4/5 genes of P. argenteus. Of these, a total of 6 miRNAs we identified that met the screening criteria. The results of prediction with miRanda revealed that novel_102, dre-miR-301c-3p, dre-miR-19b-5p, dre-miR-301b-5p, novel_113, and novel_589 could regulate the expression of the Tbx4 gene, while novel_201, dre-miR-187, dre-miR-219-3p, dre-miR-460-5p, novel_212, and novel_224 could regulate the expression of Tbx5. The binding sites of the seed sequence of the miRNAs and the 3′-UTR region of the Tbx4/5 genes were determined, and their distributions are depicted in Figs. S4 and S5 (Tbx4 GenBank ID: MH709128; Tbx5 GenBank ID: MH712458).

Construction of the dual luciferase vector and transfection

The mutated recombinant plasmid was successfully constructed using Mut Express MultiS Fast Mutagenesis Kit V2 for site-directed mutagenesis, and the results were verified by sequencing and comparison. The resulting recombinant plasmid was denoted as the pSicheck-3′-UTR-WT plasmid. The sequences of the plasmids obtained by site-directed mutagenesis are depicted in the Figs. S6 and S7. In order to verify the regulation of Tbx4/5 by the previously identified miRNAs, the psiCHECK-3′-UTR-WT plasmid was co-transfected with the miRNA mimic into HEK293T cells (TaKaRa, Zhejiang, PR China). The cells in the different experimental groups were co-transfected with different plasmids, corresponding to the different predicted regulatory miRNAs, as follows: blank control group (G1): 3′-UTR-NC + miRNA mimic-NC; negative control groups (G2 and G3): (3′-UTR-NC + miRNA mimic-NC) and (3′-UTR-WT + miRNA mimic-NC), respectively; experimental group (G4): 3′-UTR-WT + miRNA mimics (6 miRNA mimics); and point mutation control groups (G5 and G6): (3′-UTR-mut + miRNA mimic-NC and 3′-UTR-mut + miRNA mimic) and (psiCHECK-3′-UTR-mut + miRNA mimic (6 specific miRNA mimics + mimic-NC)). After 48 h of co-transfection, the regulation of Tbx4/5 by the miRNAs was detected according to the instructions provided in the Dual-Glo®Luciferase Assay System kit, and calculated as the ratio of the activity of the firefly luciferase activity to that of the Renilla luciferase, and comparing the ratios calculated for the experimental groups with those of the control group. The results demonstrated that miR-102, miR-301c, and miR-589 had a significant negative regulatory effect on the 3′-UTR of Tbx4 (P < 0.05). The 3′-UTR of Tbx4 was not significantly regulated by miR-113, miR-19b, and miR-301b; however, the 3′-UTR of Tbx5 was significantly negatively regulated by miR-187, miR-201, miR-219, and miR-460. The results also demonstrated that the 3′-UTR of Tbx5 was not significantly regulated by miR-212 and miR-224 (Figs. 4 and 5).