Increased boll number contributes to the yield heterosis of H318 in different years
H318, as an elite cotton cultivar, showed an obvious yield increase over the control variety, with an average 3,664.5 kg/ha seed cotton and 1,521 kg/ha lint yield production (9.2% higher than control) in two years (2007-2008) of regional testing.
To explore the effects of heterosis on yield, we conducted a four-year field investigation of traits such as the number of bolls, the number of fruit sites, the abscission rate and the number of fruit branches (Figure 1 A-D). At the same time, considering the climate influence on plant growth and production, the climatic changes from May to October (the whole growth period of cotton) in Wuhan from 2010 to 2013 were also obtained from the Statistics Bureau of Hubei Province (Figure S1 A-D).
In 2010, a long rainfall period and a large volume of precipitation in July (Figure S1 C, D) strongly affected the fruit setting efficiency (July and August are the flourishing flowering stage), which caused high abscission rates in 4-5 (paternal line) and B0011 (maternal line) (Figure 1 C). However, a relatively lower abscission rate and more fruit sites were found in H318 (Figure 1 B-C), which might be the reason for its final higher boll numbers (Figure 1A).
In later relatively ideal climate conditions in 2011, there was more water availability in the growing stage (June in Figure S1 C, D); suitable sunshine hours, rain and warming temperatures during the flourishing flowering stage (July and August in Figure S1 A-D); and little difference could be found in the abscission rates between the hybrid and the two parents (Figure 1 C). However, many more fruit sites were produced in hybrid H318 (Figure 1 B), which directly caused the increased boll number of H318 (Figure 1 A).
In 2012 and 2013, long-term high temperature stress that occurred in July and August (Figure S1 A) caused a higher abscission rate than the other years (Figure 1 C), but H318 still showed a much lower abscission rate than its parents, and more bolls were produced (Figure 1 A, C). No obvious differences were found in the fruit site in 2012 and 2013, and fruit branches showed no difference in any year (Figure 1 B, D). Thus, we speculate that the higher yield production of H318 directly derives from the stable higher number of bolls, which is the result of increased fruit sites or decreased abscission rate in different years.
Hybrid H318 shows heterosis of biomass and the growth rate at the seedling stage
In addition to the yield heterosis, H318 also showed growth vigor compared with its parents at the seedling stage (Figure 2 A). The fresh and dry weights (the whole plant) of H318 and its parents were measured at the two-leaf stage, and the results showed that H318 had obviously higher fresh weight than its parents and dry weight than 4-5 (Figure 2 B). Moreover, the fresh weight of B0011 was higher than 4-5, but without obvious changes in dry weight. To evaluate the plant growth status in detail, the cotyledon area was calculated every 2 days after cotyledons spread until 14 days after sowing (DAS) (Figure 2 C). At 6 DAS, hybrid H318 had a slightly larger cotyledon area than two parents. After 8 DAS, the cotyledon area significantly increased and reached its maximum area at 14 DAS. The cotyledon area of H318 remained remarkably larger than its parents after 8 DAS. These results suggested that H318 exhibits obvious heterosis of biomass production and the growth rate at the seedling stage.
Global analysis of differential gene expression in H318 and its parents
To analyze the global gene expression patterns of these three genotypes, a total of 9 RNA sequencing libraries (containing 3 independent biological repeats for each genotype) were constructed for Illumina sequencing using whole seedlings at 8 DAS. In total, at least 43,000,000 clean reads (accounting for 95% of raw reads) were generated from each sample (Table S1). All clean reads were then mapped to the cotton genome (G. hirsutum TM-1 (AD)1). Nearly 89% reads were mapped to the genome, and more than 80% reads were uniquely mapped (Table S2). Most mapped reads (> 85%) were located in exons; reads in intergenic regions occupied 11%-12%, and the rest (2.5%) were located in introns (Figure S2).
The reads mapped to the genome were then used for transcript assembly and gene expression level calculation. The expected number of fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) of all genome annotated genes (70,478) and novel genes (6,564) are listed in Table S3. The correlation analysis showed perfect consistency between biological repeats (Figure S3). All 3 possible comparisons were used for differentially expressed gene (DEG) screening with a P-value < 0.05. Finally, a total of 17,308 DEGs containing 16,047 annotated genes and 1,261 novel genes were found between three genotypes (Table S4), and many few DEGs (306) between H318 and B0011 were found (Figure S4).
To remove the DEGs with low expression or differences, a more stringent criterion (P-value < 0.05 and|log2Ratio (FPKM+1.5)| > 1 in at least one sample) was applied, leaving 3,490 DEGs. The FPKM values of the DEGs mostly ranged from 3-15 indicating that most DEGs were relatively low expressed (Figure S5), which implies low expressed DEGs were the important part for heterosis behaviour at the seedling stages. In addition, 3,472 DEGs were found between H318 and its paternal line (4-5), most DEGs (1,978) were upregulated in H318, and only 64 DEGs were found between H318 and its maternal line (B0011, Figure 3 A). Furthermore, PCA was applied with 3,490 DEGs to show the difference in gene expression between samples (Figure 3 B). The results showed that 3 genotypes were obviously divided from each other, with 4-5 located far from H318 and B0011, indicating that the expression trends of DEGs in 4-5 were much more different from those in H318 and B0011. The DEG distribution in the three genotypes showed that most DEGs (2,682/3,490) were differentially expressed between H318 and 4-5 or 4-5 and B0011; 34 DEGs were shared in the three genotypes (Figure 3 C). These results indicated that the gene expression pattern of H318 is more similar to that of its maternal line and different from that of its paternal line.
To validate the reliability of RNA-Seq, 26 DEGs were selected randomly for qRT-PCR confirmation. A correlation analysis was performed using data from RNA-Seq and qRT-PCR (Figure 3 D). The results show that the R2 value reaches 0.9791, indicating high quality for RNA-Seq.
GO and KEGG analysis
To gain a deeper understanding of the potential functions of these DEGs (3,490) in the three genotypes, GO and KEGG analyses were carried out for gene functional classification (Figure 4, Table S5 and Table S6). For GO analysis, the most significantly enriched biological processes were oxidation-reduction (green), lipid metabolic (blue), photosynthesis (red) and carbohydrate metabolism (yellow)-related processes (Figure 4, Table S5). At the same time, oxidoreductase activity (green) and photosystem-related components (red) were mostly enriched in molecular function and cellular component, respectively. For KEGG analysis, the significantly enriched pathways were various metabolic-related pathways (purple), photosynthesis (red) and carbon metabolism-related pathways (green) (Figure 4, Table S6).
Photosynthesis, lipid metabolic, carbohydrate metabolic and oxidation-reduction pathways are enriched in H318.
From the above GO and KEGG analyses, we found that photosynthesis (65 DEGs), lipid metabolism (176 DEGs), carbohydrate metabolism (178 DEGs) and oxidation-reduction (371 DEGs) pathways were significantly enriched in H318 (Table S7). To be more specific, most DEGs associated with photosynthesis had no changes (64/65) between H318 and B0011 (maternal line), but upregulated (55/65) in H318 comparing to the 4-5 (paternal line). However, most of them were downregulated in 4-5 (51/65) comparing to the B0011. For lipid (176 DEGs) and carbohydrate (178 DEGs) metabolic processes, over 100 DEGs from each pathway were upregulated in H318 comparing to 4-5, but quite a few (6/176 DEGs for lipid and 4/178 for carbohydrate metabolic processes) were found differentially expressed between H318 and B0011. In addition, 252 DEGs involved in oxidation-reduction process were found upregulated between H318 and 4-5 while 119 DEGs were downregulated. However, most of them were unchanged between H318 and B0011. The results indicated that the photosynthesis pathway was enhanced in H318 and B0011 relative to 4-5, and the lipid metabolic, carbohydrate metabolic and oxidation-reduction processes were enhanced in H318, which implied that hybrid H318 might exhibit a stronger capability for photosynthesis, lipid, carbohydrate metabolism and oxidation reduction than its parents, especially to the paternal line. We speculated that the enhanced capability of these functions might contribute to the faster growth rate and enhanced biomass of H318 at the seedling stage (Figure 2).
Higher photosynthesis rates and sugar accumulation were found in H318 at the seedling stage.
To validate the authenticity of this inference, the photosynthesis rate and the content of sucrose and starch at the two-leaf stage were assessed in both H318 and its parental lines (Figure 5). The results showed that the photosynthesis rate was significantly enhanced in H318 relative to its parental lines, especially the paternal line 4-5 (Figure 5 A). Accordingly, the contents of sucrose and starch were also higher in H318 (Figure 5 B). Moreover, the relative expression levels of 6 genes related to starch biosynthesis and photosynthesis processes were also detected and are shown in Figure 5 C, including GRANULE BOUND STARCH SYNTHASE 1 (GBSS1), RUBISCO ACTIVASE (RCA), PROTON GRADIENT REGULATION 5 (PGR5), PHOTOSYSTEM II SUBUNIT X (PSBX), PHOTOSYSTEM I SUBUNIT L (PSAL) and PHOSPHORIBULOKINASE (PRK). The expression trends of these genes were highly consistent with the photosynthesis rate and the content of sucrose and starch, which indicated an enhanced photosynthesis rate and sugar biosynthesis in H318 relative to its parental lines. All these results implied that stronger photosynthesis and metabolic processes contribute to heterosis in H318 at the seedling stage.