HL produces larger seeds and kernels.
After growing for 21 weeks, the two plant types matured and their seeds were collected. All seeds were subjected to a brief drying period in the sun, and then seed size and kernel number of WT and HL were recorded. As shown in Fig. 2, HL (B) produced larger kernels than WT (A). As Fig. 2C shows, although the seed coat of HL was obviously thicker, the large size of the seeds contributed to larger kernels in HL. It is worth noting that the seed number per plant did not show any difference between the two lines in our data (on average, WT was about 1170 seeds per plant, while HL was 1180), which indicated that the yield difference between these two lines was mainly attributed to seed size.
The difference in fruit length between the two experimental lines became significant at growth stage III.
To further understand the growth regulation of Euryale ferox Salisb., we divided development into five stages (I through V), as in rice studies[23, 24]. At stage I, fruits grow at a slow pace; at stage II, no significant growth in fruit size has occurred; at stage III, fruits show fast growth; at stage IV, fruits are nearly mature, and, at stage V, fruits have entered post-maturity (Fig. S1). We collected all fruits on the 10th week according to the growth pattern, and, more importantly, fruit growth at this phase contributed to the differences between the two lines later in development. As shown in Fig. 3, fruit length in HL was greater than in WT at stages III (p < 0.05), consequently, we selected samples at stages II and III. Further, both HL and WT produced 1–2 fruits above water level per week.
Localization of IAA in HL differed from that in WT.
Phytohormones play important roles in both plant growth and crop yield. IAA is one of the most active effectors in young tissues[25]. We measured IAA content at the initial growth stages and found no difference between fruits of the two lines at stage III (Fig. 4A). However, both concentration and localization of IAA at the cellular level, which must control cell activity and/or its fate, are critical[26]. Therefore, to investigate the localization of IAA in fruits at the same stage, immunofluorescence staining was performed on sectioned fruits. Surprisingly, at stage III, WT allocated free-IAA to its special tissue, pricks, more frequently, whereas the fluorescence signal in fruits of HL concentrated in the ovary and the stamen (Fig. 4B). Similar results were found at stage II (Fig. S2). This most interesting finding suggested a potential role for auxin signal response in regulating ovary size and, subsequently, seed size.
RNA-seq and transcriptome assembly.
A series of de novo assemblies were carried out with Trinity. To obtain useful transcriptomic information for the two lines during development, 12 cDNA libraries were established and paired-end sequence reads were generated using the Illumina Hiseq X10. As shown in Table 1, a total of 630,221,090 clean reads were generated from HL and WT cDNA libraries. The average base fraction with quality scores of Q30 and Q20 were 94.6% and 98.1%, which indicated high quality data. Relative length of assembled sequences is an important evaluation criterion for assembly quality, and the summary statistics of the assembled clean reads are shown in Table 2. Using the Trinity assembler, 197,253 unigenes were generated with an average length of 959 nt, an N50 of 1,574 nt, and a total length of 189,170,382 nt. Total unigenes included 44,511 (22.57%) of less than 300 nt, 151,012 (76.55%) with lengths from 301 to 5,000 nt, and 1,730 (0.88%) with lengths greater than 5,000 nt These results proved high quality sequencing and assembly in our samples.
Table 1
Statistical Analysis of Transcriptome Sequencing Data
Sample | Clean reads number | Clean bases | Clean rate(%) | Q20(%) | Q30(%) |
WT-S2-1 | 49,343,466 | 7,392,663,824 | 95.64 | 97.78 | 93.69 |
WT-S2-2 | 35,660,154 | 5,349,023,100 | 95.3 | 98.03 | 94.77 |
WT-S2-3 | 46,271,622 | 6,946,934,906 | 94.97 | 97.84 | 93.84 |
WT-S3-1 | 38,688,882 | 5,803,332,300 | 95.66 | 98.11 | 94.98 |
WT-S3-2 | 53,603,992 | 8,023,343,356 | 97.29 | 97.69 | 93.01 |
WT-S3-3 | 65,400,866 | 9,791,552,137 | 96.21 | 97.77 | 93.57 |
HL-S2-1 | 106,884,270 | 16,025,803,643 | 97.05 | 97.96 | 93.96 |
HL-S2-2 | 31,283,988 | 4,692,598,200 | 95.41 | 98.07 | 94.89 |
HL-S2-3 | 48,644,712 | 7,296,706,800 | 96.35 | 98.42 | 95.62 |
HL-S3-1 | 42,801,744 | 6,420,261,600 | 96.21 | 98.4 | 95.59 |
HL-S3-2 | 40,785,442 | 6,117,816,300 | 96.05 | 98.43 | 95.66 |
HL-S3-3 | 46,852,120 | 7,027,818,000 | 96.41 | 98.43 | 95.66 |
Table 2
Length distribution of the unigenes in transcriptome
Length range | Unigenes | Percentage (%) |
200–300 | 44,511 | 22.57 |
300–500 | 43,507 | 22.06 |
500-1,000 | 48,729 | 24.70 |
1,000–2,000 | 36,800 | 18.66 |
2,000–5,000 | 21,976 | 11.14 |
5,000+ | 1,730 | 0.88 |
Total number | 197,253 | |
Total length | 189,170,382 | |
N40 length | 1,980 | |
N50 length | 1,574 | |
N60 length | 1,220 | |
N70 length | 910 | |
N80 length | 633 | |
N90 length | 390 | |
Min length | 201 | |
Max length | 17,279 | |
Mean length | 959 | |
Annotation and GO classification.
The annotation of the 197,253 assembled unigenes revealed that 44,363 (22.49%), 18,229 (9.24%), 26,746 (13.56%), 9,633 (4.88%), 3,632 (1.84%), and 2,879 (1.46%) showed significant similarity in the GO, Pathway Pfam, TMhmm, Eggnog, SignalP, and Uniprot databases, respectively. In all, 71,399 unigenes (36.20% of all unigenes) were successfully annotated in at least one database (Table 3). As shown in Fig. 5, total 44363 unigenes annotated in GO were assigned to three principal domains; “cellular component”, “molecular function”, and “biological process”.
Table 3
Unigenes annotated in databases
Database | Number of unigenes | Percentage (%) |
GO | 44,363 | 22.49 |
Pathway Pfam | 18,229 | 9.24 |
TMhmm | 26,746 | 13.56 |
Eggnog | 9,633 | 4.88 |
signalP | 3,632 | 1.84 |
Uniprot | 2,879 | 1.46 |
All database | 68,200 | 34.57 |
At least one database | 71,399 | 36.2 |
Total unigenes | 197,253 | 100 |
SAURs induced differences in the IAA signal transduction pathway between HL and WT.
The KEGG Pathway database records the networks of molecular interactions in cells and species-specific variations[27]. According to the findings on localization of IAA reported above, we focused on genes involved in the auxin transduction pathway. The pathway is shown in Fig. 6A; briefly, auxin promotes an interaction between TIR1/AFB and Aux/IAA proteins, resulting in the degradation of the Aux/IAAs and the release of ARF repression, and down-stream SAURs, GH3 and AUX/IAA are activated subsequently. We compared the transcriptome from HL with that from WT at both stages II and III: only SAURs showed a significant difference in the whole pathway (Fig. 6A). Then, all unigenes annotated with ARF, AUX/IAA, GH3 or SAUR were selected for comparison (Fig. 6B). An obvious difference was observed in SAURs, while the other three families did not show any changes between the two lines at either growth stage. It is worth noting that the unigene number of SAURs is the least among the four families, which indicated less functional redundancy in SAURs of Euryale ferox Salisb. Interestingly, this finding suggests an important role for SAURs in responding to IAA and, probably, up-regulating seed size in HL.
Verification of RNA-Seq by qRT-PCR.
To evaluate the validity of the expression of DEGs, 10 candidate unigenes involved in the auxin transduction pathway (2-ARFs, 3-AUX/IAAs, 2-GH3s, and 4-SAURs) were detected by qRT-PCR (Detailed information about all 10 unigenes is given in Table S1.). Fruits from HL and WT at stage III were selected for qRT-PCR analysis, and β-actin was used as a reference gene in this study. As shown in Fig. 7, most of these DEGs showed a similar expression pattern according to qRT-PCR and RNA-Seq results, which suggested that results from the transcriptome data are reliable.