Histochemistry and morphometric analysis
TRAcP staining of scales
The results of TRAcP staining on the scales of juvenile O. mykiss were shown in Fig. 1, the positive TRAcP sites were red and mainly located at the scale margins. As for the CG scales, the area of positive TRAcP staining sites was small. Compared with CG scales, the TRAcP positive staining area and the number of positive reaction sites of the scales in 7D and 21D group were larger.
ALP staining of scales
The results of ALP staining on the scales of juvenile O. mykiss were shown in Fig. 2, and the ALP positive sites were stained purple. A large positive staining area was found in the middle of the scales, and no significant difference was identified in the scales of different groups.
Von Kossa′s staining of scales
The results of Von Kossa′s staining on the scales of juvenile O. mykiss were shown in Fig. 3, and the positive sites of calcium salts were stained brown. A large area of positive Von Kossa′s staining was identified in the middle of the scales. The positive Von Kossa′s staining in 7D scales was stronger than scales of the other three groups. No resorption pit was observed on scales of all the four groups. The location of Von Kossa′s staining area was similar to that of the ALP staining area, and the calcium salt deposits less in the TRAcP positive staining area.
Calcium and phosphorus analysis of scales
Changes in Ca content of the scales throughout the experimental period were shown in Fig. 4 (A), and the result was consistent with the Von Kossa′s staining results. The Ca content first increased and reached the highest level (0.100 µmol/scale) at 7D. The Ca content of the scales at 7D, 14D and 21D were significantly higher than that of the CG scales (P-value < 0.05). The P content of the scales [Fig. 4 (B)] also reached the highest level (0.104 µmol/scale) at 7D, which was significantly higher than that of the CG scales (0.081 µmol/scale) (P-value < 0.05); the 14D scales exhibited the lowest P content (0.067 μmol/scale) (P-value < 0.05).
The molar mass ratio of calcium to phosphorus (Ca/P), indicating the crystalline phase of calcium phosphate (e.g., pure calcium hydroxyapatite, theoretical ratio = 1.67), was also shown in Fig. 4 (C). The Ca/P of 14D scales was significantly higher than scales of the other three groups, and no significant difference was found among Ca/P of CG, 7D and 21D scales.
Illumina sequencing of miRNAs in O. mykiss
To identify miRNAs in the scales of juvenile O. mykiss, twelve miRNA libraries from CG (CG_1, CG_2 and CG_3), 7D (7D_1, 7D_2 and 7D_3), 14D (14D_1, 14D_2 and 14D_3) and 21D (21D_1, 21D_2 and 21D_3) were respectively constructed and sequenced using Illumina sequencing technology. A total number of 11,277,023, 11,052,110, 12,969,211 and 11,929,493 raw reads were respectively obtained from CG, 7D, 14D and 21D scales. After preprocessing steps, a total of 921,415, 824,625, 946,666 and 894,351 clean reads were obtained due to the removal of low-quality reads, adapters et al. The aforementioned reads respectively represented 341,328, 261,252, 245,627 and 236,173 unique sequences (valid reads) (Table 2). Furthermore, the length distribution analysis of these miRNA sequences showed a similar pattern of distribution in length of all libraries (Fig. 5). The length of the miRNAs in all libraries varied from 18 to 26 nt and the majority of read length was 22 and 23 nt, followed by 21, 24, 25 and 26 nt. The read length of 22 nt in the four groups respectively accounted for 23.31%, 24.75%, 33.40% and 35.18% of the total miRNA numbers. The raw reads of the twelve libraries were uploaded into the NCBI database Sequence Read Archive (SRA) and the SRR numbers were from SRR15559621 to SRR15559632.
Table 1 Primer sequences used for qRT-PCR
Primer name
|
Primer sequence (5′~3′)
|
DRE-MIR-193B-3P-F
|
CCCGCAAAGTCCCGCTAAA
|
DRE-MIR-193B-3P-R
|
CAGTGCAGGGTCCGAGGTAT
|
SSA-MIR-203A-3P-F
|
GTGAAATGTTTAGGACCACTTG
|
SSA-MIR-203A-3P-R
|
CAGTGCAGGGTCCGAGGTAT
|
SSA-MIR-1338-3P-F
|
TCAGGTTCGTCAGCCCATG
|
SSA-MIR-1338-3P-R
|
CAGTGCAGGGTCCGAGGTAT
|
SSA-MIR-92A-3P-F
|
TATTGCACTTGTCCCGGCCT
|
SSA-MIR-92A-3P-R
|
CAGTGCAGGGTCCGAGGTAT
|
SSA-MIR-205A-5P-F
|
TTCCTTCATTCCACCGGATC
|
SSA-MIR-205A-5P-R
|
CAGTGCAGGGTCCGAGGTAT
|
ONI-MIR-125A-F
|
CCTGAGACCCTTAACCTGTG
|
ONI-MIR-125A-R
|
CAGTGCAGGGTCCGAGGTAT
|
U6-F
|
CTCGCTTCGGCAGCACATATACT
|
U6-R
|
ACGCTTCACGAATTTGCGTGTC
|
Table 2Analysis of miRNA sequences of juvenile O. mykiss scales collected at different time points during salinity acclimation
Type
|
CG
|
7D
|
14D
|
21D
|
Total
|
%
|
unique
|
%
|
Total
|
%
|
unique
|
%
|
Total
|
%
|
unique
|
%
|
Total
|
%
|
unique
|
%
|
Raw reads
|
11,277,023
|
100.00
|
921,415
|
100.00
|
11,052,110
|
100.00
|
824,625
|
100.00
|
12,969,211
|
100.00
|
946,666
|
100.00
|
11,929,493
|
100.00
|
894,351
|
100.00
|
3ADT& length filter
|
5,251,562
|
46.50
|
567,743
|
61.85
|
5,157,869
|
45.88
|
552,577
|
66.90
|
6,958,124
|
53.57
|
689,174
|
72.84
|
6,294,338
|
52.96
|
644,677
|
72.06
|
Junk reads
|
30,525
|
0.27
|
3,428
|
0.37
|
23,126
|
0.21
|
2,505
|
0.30
|
9,056
|
0.07
|
2,208
|
0.23
|
8,140
|
0.07
|
1,738
|
0.19
|
Rfam
|
542,449
|
4.83
|
7,238
|
0.79
|
564,990
|
5.15
|
7,133
|
0.88
|
509,731
|
3.94
|
8,623
|
0.92
|
679,638
|
5.69
|
10,551
|
1.19
|
mRNA
|
41,515
|
0.37
|
1,721
|
0.18
|
33,869
|
0.32
|
1,167
|
0.14
|
46,402
|
0.36
|
1,114
|
0.12
|
66,780
|
0.55
|
1,308
|
0.15
|
Repeats
|
109,582
|
0.98
|
517
|
0.06
|
154,111
|
1.37
|
481
|
0.06
|
66,809
|
0.51
|
787
|
0.08
|
70,264
|
0.58
|
891
|
0.10
|
valid reads
|
5,412,133
|
48.04
|
341,328
|
36.81
|
5,272,049
|
48.45
|
261,252
|
31.77
|
5,447,228
|
42.07
|
245,627
|
25.90
|
4,880,402
|
40.72
|
236,173
|
26.42
|
Identification and specific expression of miRNAs in O. mykiss scales
To identify miRNAs in the scales of juvenile O. mykiss, all valid sequences were compared with known vertebrate miRNAs and miRNA precursor sequences in miRBase database. As a result, a total of 756 (664 known and 92 putative novel miRNAs) miRNAs were identified (Table S1). The read numbers of these known miRNAs ranged from 1 to 483,390, indicating that there were significantly difference in expression levels of different miRNAs. These miRNAs covered 121 miRNA families, and the most abundant families were let-7 (40 members), miR-10 (29 members) and miR-30 (24 members) (Table S2).
In the 756 miRNAs, 514 miRNAs showed co-expressed in all the four groups (Fig. 6). Eight,18, 29 and 49 miRNAs were uniquely expressed in CG, 7D, 14D and 21D scales, respectively. In all the four groups, 21D scales exhibited the largest number of unique miRNAs. 594 miRNAs were expressed in scales of CG, 607 were expressed in 7D, 636 were expressed in 14D, and 645 were expressed in 21D.
Transcriptome sequencing data
RNA-Seq was performed in juvenile O. mykiss scales of each sample, and a total of 147.29 Gb valid data was obtained. Quality of the RNA-Seq data and comparison of mRNA sequences from different groups were shown in Table 3. The valid data of each sample reached 11.22 Gb, and the minimum percentage of Q30 bases was 97.32%. The valid data of each sample was compared with the O. mykiss genome, and the minimum comparison efficiency reached 82.74%.
Table 3Statistics of RNA-seq data from each juvenile O. mykiss scale sample collected at different time points during salinity acclimation
Sample
|
Raw Read
|
Valid Read
|
Valid data/G
|
Q30/%
|
Mapped reads
|
CG_1
|
89,627,272
|
83,437,990
|
12.52G
|
98.71
|
72,429,195 (86.81%)
|
CG_2
|
90,637,284
|
87,198,854
|
13.08G
|
98.58
|
78,893,393 (90.48%)
|
CG_3
|
80,356,512
|
74,812,354
|
11.22G
|
98.25
|
64,591,320 (86.34%)
|
7D_1
|
92,595,784
|
82,794,096
|
12.42G
|
98.59
|
71,541,252 (86.41%)
|
7D_2
|
95,578,444
|
89,211,058
|
13.38G
|
98.09
|
76,755,268 (86.04%)
|
7D_3
|
90,741,298
|
83,933,042
|
12.59G
|
97.99
|
73,039,003 (87.02%)
|
14D_1
|
91,322,682
|
76,219,276
|
11.43G
|
98.49
|
63,555,106 (83.38%)
|
14D_2
|
94,608,098
|
75,339,848
|
11.30G
|
98.58
|
62,339,088 (82.74%)
|
14D_3
|
90,570,502
|
80,280,090
|
12.04G
|
98.54
|
67,911,295 (84.59%)
|
21D_1
|
94,544,622
|
76,698,318
|
11.50G
|
98.56
|
67,600,190 (88.14%)
|
21D_2
|
91,295,438
|
86,198,836
|
12.93G
|
97.37
|
73,723,516 (85.53%)
|
21D_3
|
88,674,894
|
85,859,784
|
12.88G
|
97.32
|
72,851,690 (84.85%)
|
Differential expression analysis and target genes prediction of the identified miRNAs
The trend of miRNA expression at different timepoints (CG, 7D, 14D and 21D) was analyzed using STEM software. As shown in Fig. 7, the 327 DE miRNAs were clustered into 30 expression patterns, of which five expression patterns were with significant changing trends (Cluster 10, Cluster 11, Cluster 17, Cluster 18 and Cluster 25) (P-value < 0.05). A total of 290 DEC miRNAs (DE miRNAs from clusters with significant trends in the cluster analysis, DEC miRNAs) were clustered into these five patterns (contains 22 DEC miRNAs in Cluster10; Cluster 11 contains 37 DEC miRNAs; Cluster 17 contains 91 DEC miRNAs; Cluster 18 contains 70 DEC miRNAs; Cluster 25 contains 70 DEC miRNAs), and detailed information was listed in Table S3. A total of 22,374 predicted target genes of the DEC miRNAs from the aforementioned five expression modes were obtained; the DEC miRNAs in Cluster 10, Cluster 11, Cluster 17, Cluster 18 and Cluster 25 each had 3,897, 5,756, 9,423, 8,586 and 8,977 predicted mRNAs.
Additionally, trend analysis was performed on 22,374 predicted target genes of DEC miRNAs. The co-clustering was also divided into 30 expression patterns, among which 5,957 DEC mRNAs (DE mRNAs of 22,370 DEC miRNA target genes showed a significant trend in cluster analysis, DEC mRNAs) were clustered into 6 expression patterns with significant changing trends (Cluster 16, Cluster 17, Cluster 22, Cluster 24, Cluster 27 and Cluster 29) (P-value < 0.05) (Fig. 8). Detailed information was listed in Table S4. Cluster 16 contains 682 DEC mRNAs; Cluster 17 contains 1,140 DEC mRNAs; Cluster 22 contains 387 DEC mRNAs; Cluster 24 contains 1,092 DEC mRNAs; Cluster 27 contains 740 DEC mRNAs; Cluster 29 contains 1,916 DEC mRNAs.
GO analysis of DEC miRNA-targeted genes
Gene ontology (GO) enrichment (function) analysis was performed on predicted DEC mRNAs to explain and speculate the function of DEC mRNAs (Fig. 9). Detailed information of the GO analysis was listed in Table S5. As shown in Fig. 9, the predicted DEC mRNAs could be classified into three major categories: biological process (BP), cellular component (CC), and molecular function (MF). The columns in each category presented in the figure were in descending order of percent of genes. The DEC mRNAs in the BP category were enriched in regulation of transcription, DNA-templated, phosphorylation, signal transduction, transport, et al. The DEC mRNAs in the CC category were associated with terms such as membrane, integral component of membrane, cytoplasm, nucleus and cytosol. In the MF category, the DEC mRNAs were associated with terms such as metal ion binding, ATP binding, nucleotide binding and transferase activity.
KEGG pathway analysis of DEC miRNA-targeted genes
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed to analyze the biological pathways of the predicted DEC mRNAs (Fig. 10). Detailed information of the KEGG pathways (P-value< 0.05) was listed in Table S6. predicted DEC mRNAs were significantly enriched in the pathways such as the MAPK signaling pathway, Toll-like receptor signaling pathway, mTOR signaling pathway, Thyroid hormone signaling pathway, Calcium signaling pathway, VEGF signaling pathway, ErbB signaling pathway, Insulin secretion, Wnt signaling pathway, TGF-beta signaling pathway, FoxO signaling pathway, Insulin signaling pathway, Mineral absorption, NF-kappa B (NF-κB) signaling pathway and other bone metabolism related signaling pathways.
Validation of DE miRNAs using qRT-PCR
A total of six randomly selected miRNAs, including DRE-MIR-193B-3P, SSA-MIR-203A-3P, SSA-MIR-1338-3P, SSA-MIR-92A-3P, SSA-MIR-205A-5P and ONI-MIR-125A, were used for qRT-PCR analysis. As shown in Fig. 11, the expression pattern of the six miRNAs were completely identical with the transcriptome sequencing results.