Origin of reciprocal cross hybrids
We first characterized the divergence of AS between two reciprocal crosses hybrids (BTF3 and TBF3), which were obtained from the self-crossing of respective reciprocal crosses hybrids of M. amblycephala (2n = 48) ´ C. alburnus (2n = 48) [20, 21]. The genotype of chimeric offspring was determined as the allodiploid (2n = 48) with a 1:1 subgenome ratio with chromosome number and 45S rDNA characteristics [20], in which the two types of 45S rRNA were detected and belonged to species-specific sequences of M. amblycephala and C. alburnus, respectively [21] (Additional file 1: Table S1). The expression of mitochondrial genes in the two reciprocal crosses hybrids was considered to be identical to that of the respective inbred female parents based on the mapped reads of transcriptome (Additional file 2: Table S2).
Characteristic differences of the two subgenomes
For the two inbred parental genomes (1.09 Gb in BSB and 1.02 Gb in TC), the distributions of exon numbers and CDS lengths were obtained from 20,131 orthologous gene pairs (Additional file 3: Fig. S1). The average exon number in each gene was 8.83 in BSB and 9.64 in TC, while the average CDS length was 1525 bp in BSB and 1654 bp in TC. Focusing on same characteristic in the two parental genomes, the same exon number was found in 11,414 genes, and the same CDS length was detected in 6832 genes. Analysis of these results showed significant differences in exon number (p < 0.001) and CDS length (p < 0.001). Furthermore, strong associations with exon number and CDS length were detected in BSB (r = 0.7435) and TC (r = 0.7768) (all p < 0.0001 for Pearson correlation coefficients).
Obtaining of long length transcripts and Gene Ontology Analysis
PacBio sequencing was used to detect AS events in reciprocal crosses hybrids. A total of 21.22 Gb and 15.49 Gb data were obtained from TBF3 and BTF3, respectively, and an average of 12 and 13 CCSs and 3080 bp and 2936 bp average insert read lengths were detected in TBF3 and BTF3, respectively (Table 1). To detect the integrity of sequencing data, 663,834 TBF3 and 479,667 BTF3 3ʹ and 5ʹ- untranslated regions were analyzed to determine whether the transcripts were full-length. Then, 586,075 (88.29%) and 431,999 (90.06%) full-length reads were detected in TBF3 and BTF3 (Table 1). After deleting redundant sequences, a total of 316,533 and 268,986 consensus reads were found in TBF3 and BTF3, respectively.
After mapping to the combined genome of the two inbred parents, the 314,298 consensus reads and 76,518 isoforms were obtained from the mapped results of TBF3, while 267,949 consensus reads with 82,083 isoforms were found in BTF3. An average 99.29% and 99.61% of mapping ratios were detected in TBF3 and BTF3, respectively. Then, the sequences of 11,026 genes in TBF3 and 11,448 genes in BTF3 were annotated with KEGG and GO databases (Additional file 4: Fig. S2). The 6071 genes shared between TBF3 and BTF3 were then focused on to help detect differences between the two.
AS between two homoeologs
To better characterize the differences between TBF3 and BTF3, we focused on AS events in BSB- and TC- homoeologs of the two reciprocal cross hybrids. A custom Python script was used to identify 30,007 AS events, and 7029 genes were mapped to the BSB-subgenome detected in 20,131 homoeologous genes of TBF3, while 21,577 AS events and 5286 genes were mapped to the TC-subgenome (Table 2). We also detected 30,305 AS events and mapped 7271 genes to the BSB-subgenome of BTF3, while 30,562 AS events related to 6481 genes were mapped to the TC-subgenome (Table 2). Although the sequencing was performed in a mixture of three tissues, these data suggest a slight BSB-homoeolog expression bias in TBF3 but not in BTF3. Most of the AS events that occurred in hybrids were RI in TC-homoeologs of TBF3 (18.96%) and BTF3 (26.49%), and BSB-homoeologs of BTF3 (26.70%). However, most AS events in BSB-homoeologs of TBF3 were SE (20.28%) (Table 2). Although some error in gene annotation may have led to an increased prediction of RI and SE, these results suggest a clear difference not only between the two reciprocal cross hybrids, but also between BSB- and TC-homoeologous genes. Then, we compared the number of AS events between in each orthologous gene pairs of two subgenomes. Among these, 1290 genes of TBF3 and 2302 genes of BTF3 were shown to possess AS events in both homoeologs. In TBF3, 4862 AS events supported by 6416 isoforms were mapped to the BSB-subgenome, while 4650 AS events supported by 6674 isoforms were mapped to the TC-subgenome (Fig. 1). In addition, we detected the 8054 AS events shared in orthologous gene pairs of BSB-subgenome and TC-subgenome of TBF3, while the 11024 AS events were shared in ones of BTF3.
Expression changes led by maternal effects
Comparison between BTF3 and TBF3, we identified 49 differentially expressed genes (DEGs) in BSB-homoeologous genes of liver, 186 DEGs in muscle, and 348 DEGs in gonads; this compared with 54 DEGs in TC-homoeologous genes of liver, 204 in muscle, and 354 in gonads (Fig. 2, Additional file 5: Table S3). The largest number of DEGs was found in gonads (3.58% in BSB-homoeologs and 3.64% in TC-homoeologs) and the fewest were found in liver (0.50% in BSB-homoeologs and 0.56% in TC-homoeologs) (Fig. 2, Additional file 5: Table S3).
We next focused on DEGs that were shared between BSB- and TC-homoeologs. The same up/down-regulated expression trends of the two homoeologs were exhibited among the three tissues (Fig. 2, Additional file 5: Table S3), indicating that similar differential expression trends occurred in both homoeologs. GO analysis showed that 90, 12, and 51 DEGs (the largest number in GO categories) were involved in the cellular process (GO: 0009987) (level 2) in gonad, liver, and muscle tissues, respectively (Additional file 6: Fig. S3 and S4). Some DEGs were also enriched in other functions, including metabolic process (GO: 0008152), response to stimulus (GO: 0050896), and biological regulation (GO: 0065007), while others were enriched in growth (GO: 0040007), immune system process (GO: 0002376), and reproduction (GO: 0000003) (Additional file 6: Fig. S3 and S4).
Determination of AS in DEGs
To further investigate the maternal effects on expression divergence, AS analysis was performed in homoeologous genes between the reciprocal cross hybrids. ASprofile detected 104 MXE and 3314 SE in gonads, 96 MXE and 2863 SE in liver, and 74 MXE and 2328 SE in muscle (Additional file 7: Table S4). Interestingly, 3103 (90.78%) AS events were found in TC-homoeologous genes of gonads, while 315 (9.22%) AS events were found in BSB-homoeologous genes. Moreover, 2706 (91.45%) AS events were distributed in TC-homoeologous genes of liver, and the remaining 253 (8.55%) AS events occurred in BSB-homoeologous genes. In muscle, 2205 (91.80%) AS events occurred in TC-homoeologous genes compared with 197 (8.20%) AS events in BSB-homoeologous genes (Additional file 7: Table S4). However, no RI, A5SS, or A3SS were identified in Illumina data. A total of 41, 31, and 22 genes were detected as high AS events (number of AS types ≥5 in each gene) from muscle, liver, and gonads, respectively. Among these, the five genes (ptprm, cast, exoc, myo1b, and abi1a) with a high number of AS events were shared among the three tissues (Additional file 8: Fig. S5).
Combining AS and DEG analysis, we identified AS events in 35 DEGs in gonads, 18 DEGs in muscle, and six DEGs in liver. Under different maternal effects, changes to homoeologous gene expression and AS events were found in reciprocal cross hybrids. However, the details of some AS events were inaccurate using Illumina data because of the short length of the reads. Therefore, long length reads of BTF3 PacBio data were used to improve the analysis of AS events in DEGs, including 38 AS events in gonads, 16 in muscle, and two in liver, while TBF3 data improved AS events in 33, 14, and six DEGs in gonads, muscle, and liver, respectively.
AS distribution in bone morphology
In view of the many shared traits between the reciprocal cross hybrids, the observed slight differences in their appearance made us consider their control of bone morphology regulated gene expression. Focusing on the BMP family, 17 orthologous genes were obtained from BSB and TC genomes, which exhibited gene expansion events (Fig. 3A). However, BSB- and TC- homoeologous gene expression in the three tissues studied was only detected in bmpr2a simultaneously. Slight differences in homoeologous gene expression were observed between TBF3 and BTF3 (Fig. 3B), although these were not significant. Therefore, bmpr2a was selected as a model to investigate AS events in BSB- and TC- homoeologs. We identified 12 exons, which was identical to the zebrafish (Fig. 3C) [22]. SE differences of two, one, and zero were detected between TBF3 and BTF3 in gonads, muscle, and liver, respectively, using Illumina data, but fewer AS events were detected in PacBio data. However, longer length transcripts provided more accurate AS prediction in respective BSB- and TC- homoeologs. In TC- homoeologs of bmpr2a, two RI events were observed between exons 9 and 10 and between exons 3 and 4 in TBF3, while three A3SS events and one SE were detected (Fig. 4). Unfortunately, we only detected one isoform because of the potential sequencing bias. In BSB-homoeologs of bmpr2a, one SE event was identified that led to the loss of sequences in exons 7–11 of TBF3. Interestingly, we also identified an SE event with the loss of sequences in exons 2–11 in BTF3, similar to the SE event in the TC-homoeologous gene of TBF3 (Fig. 4).
For further determination of AS events in bmpr2a, the transcripts in the muscle of TB and BT were sequenced by Sanger method. The one RI event between exons 9 and 10 (AS_2) and the two SE events distributed in parts of exon 12 (AS_3 and 4) were detected in TC-homoeologs of BTF3 and TBF3. These phenomena were same to the above results of PacBio data (Fig. 4; Additional file 9: Fig. S6). Furthermore, one SE in exon 7 (AS_1) was observed in BSB-homoeologs of BTF3 and TBF3 (Additional file 9: Fig. S6).