Heat-responsive microRNAs participate in regulating the pollen fertility stability of CMS-D2 restorer line under high-temperature stress

Abstract

Anther development and pollen fertility of cytoplasmic male sterility (CMS) conditioned by Gossypium harknessii cytoplasm (CMS-D2) restorer lines are susceptible to continuous high-temperature (HT) stress in summer, which seriously hinders the large-scale application of "three-line" hybrids in production. Here, integrated small RNA, transcriptome, degradome, and hormone profiling was performed to explore the roles of microRNAs (miRNAs) in regulating fertility stability in mature pollens of isonuclear alloplasmic near-isogenic restorer lines NH and SH under HT stress at two environments. A total of 211 known and 248 novel miRNAs were identified, of which 159 were differentially expressed miRNAs (DEMs). Additionally, 45 DEMs in 39 miRNA clusters (PmCs) were also identified, and most highly expressed miRNAs were significantly induced in SH under extreme HT, especially four MIR482 and six MIR6300 family miRNAs. PmC28 was located in the fine-mapped interval of the Rf₁ gene and contained two DEMs, gra-miR482_L-2R + 2 and gma-miR2118a-3p_R + 1_1ss18TG. Transcriptome sequencing identified 6281 differentially expressed genes, of which heat shock protein (HSP)-related genes, such as HSP70, HSP22, HSP18.5-C, HSP18.2 and HSP17.3-B, presented significantly reduced expression levels in SH under HT stress. Through integrating multi-omics data, we constructed a comprehensive molecular network of miRNA-mRNA-gene-KEGG containing 35 pairs of miRNA/target genes involved in regulating the pollen development in response to HT, among which the mtr-miR167a_R + 1, tcc-miR167c and ghr-miR390a, tcc-miR396c_L-1 and ghr-MIR169b-p3_1ss6AG regulated the pollen fertility by influencing ARF8 responsible for the auxin signal transduction, ascorbate and aldarate metabolism, and the sugar and lipid metabolism and transport pathways, respectively. Further combination with hormone analysis revealed that HT-induced jasmonic acid signaling could activate the expression of downstream auxin synthesis-related genes and cause excessive auxin accumulation, followed by a cascade of auxin signal transduction, ultimately resulting in pollen abortion. The results provide a new understanding of how heat-responsive miRNAs regulate the stability of fertility restoration for CMS-D2 cotton under heat stress.

Keywords: CMS-D2 restorer line; Degradome; High-temperature stress; MiRNA cluster; Plant hormone signal transduction; Pollen fertility stability.

Fig. 1

Comparison of the field performance of flowers and anthers from NH and its isonuclear alloplasmic near-isogenic line SH under mild HT (AY) and extreme HT (JJ) stress. A A representative intact flower from NH (left) and SH (right) under mild HT in AY. B Anthers from NH (left) and SH (right), showing normal pollen release only from NH under mild HT in AY. C, D Pollen grains from NH (C) and SH (D) plants under mild HT stress and stained with Benzidine-α-Naphthol in AY. E A representative intact flower from NH (left) and SH (right) under extreme HT in JJ. F Anthers from NH (left) and SH (right), showing normal pollen release only from NH under extreme HT in JJ. G, H Pollen grains from NH (G) and SH (H) plants under extreme HT stress and stained with Benzidine-α-Naphthol in JJ. More sterile pollen grains from SH were observed under HT stress, especially in JJ spot. Fertile pollen is stained red, part of the vitality shows reddish, and sterile pollen is colorless. Scale bars = 5 mm in A, B, E, F and 100 µm in C, D, G, H. I–K Graphical representations of the percentage of anther dehiscence (I), filament length (J), and exposed length of stigma (K) of NH and SH under mild and extreme HT stress in AY and JJ spots, respectively. AY Anyang, JJ Jiujiang

Fig. 2

Overview of small RNA sequencing data and length distribution. A Statistics of total and unique sRNA reads in 12 sequencing libraries. B, C) Annotation distribution of the total B and unique C reads in each sample. other, other Rfam RNA; snRNA, small nuclear RNA; snoRNA, small nucleolar RNA; tRNA, transfer RNA; rRNA, ribosomal RNA; repeat, repeat associate RNA; mRNA, messenger RNA; VsRNA, valid small RNA. D, E Length distribution of 18–25 nt small RNAs identified in total (D) and unique € reads. The X-axis indicates sRNAs of different lengths, while the Y-axis represents the percentage of sRNAs at a certain length. F Ratio of 24/21 nt sRNAs. The ratio of NH and SH at the JJ spot is significantly lower than that of AY. Three biological replicates are merged into the mean value. AY, Anyang; JJ, Jiujiang. AP_NH, NH under mild HT stress; AP_SH, SH under mild HT stress; JP_NH, NH under extreme HT stress; JP_SH, SH under extreme HT stress

Fig. 3

Statistics of identified miRNAs and analysis of their base biases in the sequencing libraries. A Number of pre-miRNAs and unique miRNAs in different categories. B Length distribution of all identified unique miRNAs. C UpSet Venn diagram showing the number of shared and specific miRNAs in different samples. AP_NH, NH under mild HT stress; AP_SH, SH under mild HT stress; JP_NH, NH under extreme HT stress; JP_SH, SH under extreme HT stress. D The first base preference of mature miRNAs. The X-axis displays the length classification of the miRNAs, and the Y-axis represents the proportion of mature miRNAs with a certain base type as the first base. E Base preference in different positions of mature miRNAs. The X-axis displays the different base positions of the miRNAs, and the Y-axis denotes the proportion of bases at a certain position of mature miRNAs

Fig. 4

Identification of differentially expressed miRNAs in NH and SH under mild and extreme HT stress. A Number of DEMs that were up- or down-regulated under mild and extreme HT stress. B Venn diagram showing the number of unique and shared DEMs. C Frequency percentage of DEMs with different expression levels. D Hierarchical cluster analysis of the identified 126 DEMs using a bending heat map drawn by TBtools software. The relative expression levels in the figure are the normalized expression data displayed with log₁₀ (norm + 1) value. AP_NH, NH under mild HT stress; AP_SH, SH under mild HT stress; JP_NH, NH under extreme HT stress; JP_SH, SH under extreme HT stress

Fig. 5

Identification and expression analysis of miRNA clusters. A Location distribution of all expressed miRNAs and DEMs contained in Pre-miRNA Clusters (PmCs) on different chromosomes of upland cotton. Different colors indicate different categories of miRNAs, namely black, red, purple, green, olive, and cyan belong to the miRNAs of ‘gp1a’, ‘gp1b’ ‘gp2a’, ‘gp2b’, ‘gp3’ and ‘gp4’ categories, respectively. Among them, miRNAs with bold italics and enlarged fonts and underlined are identified DEMs. The blue columns on the left side of each chromosome or scaffold represent the specific genomic position of all 39 PmCs (Inter-distance < = 50,000 nts), and in particular, the blue column on the left side of Chir_D05 chromosome signifies the mapped interval of the fertility restorer gene Rf₁ in our previous study [6]. B Statistics of the number of PmCs and their distribution of miRNAs and DEMs on different chromosomes. The X-axis represents chromosome name; the Y-axis and the numbers above each bar represent the PmC and miRNA numbers on each chromosome. C Venn diagram showing the number of DEMs in PmCs. D Hierarchical cluster analysis of the identified 39 DEMs in PmCs using a heat map drawn by TBtools software. The relative expression levels in the figure are the normalized expression data displayed with log₁₀ (norm + 1) value. AP_NH, NH under mild HT stress; AP_SH, SH under mild HT stress; JP_NH, NH under extreme HT stress; JP_SH, SH under extreme HT stress

Fig. 6

Identification and enrichment analysis of DEGs in NH and SH under mild and extreme HT stress. A Number of DEGs that were up- or down-regulated under mild and extreme HT stress. B Venn diagram showing the number of unique and shared DEGs. C Hierarchical cluster analysis of the top 100 DEGs using a heat map. The abscissa is the sample, and the ordinate is the gene name. The relative expression levels in the figure are the normalized expression data by the Z-value method and displayed with log₁₀ (FPKM + 1) value; red and dark blue indicate high and low expressed genes, respectively, and the colour bar is on the right of the heat map. D, E GO (D) and KEGG (E) enrichment analysis of DEGs. The top 20 GO terms or enriched pathways are plotted based on the significant P-value. The size of the circle represents the number of genes, and the color of the circle signifies the P-value. The X-axis denotes the enrichment factor, which compares the ratio of genes annotated to a GO term or KEGG pathway among the identified DEGs to the ratio of genes annotated to that GO or pathway among all genes, and the Y-axis indicates the GO or pathway name. AP_NH, NH under mild HT stress; AP_SH, SH under mild HT stress; JP_NH, NH under extreme HT stress; JP_SH, SH under extreme HT stress

Fig. 7

Target identification of miRNAs by degradome sequencing and GO and KEGG enrichment. A Number of miRNAs and their target genes. B, C Venn diagrams showing the number of unique and shared miRNAs B and their target genes C in P_NH and P_SH. D–F Categories and statistics of miRNAs (D), cleavage sites (E), and involving transcripts (F). G Target plot (T-plot) showing that gma-miR6300 cleaves the Ghir_A10G015260.1 transcript at the 168th nucleotide position. H T-plot showing that gra-miR482d cleaves the Ghir_D10G019110.1 transcript at the 306th nucleotide position. I T-plot showing that ghr-MIR2948-p3 cleaves the Ghir_A02G011540.1 transcript at the 145th nucleotide position. J, K GO (J) and KEGG (K) enrichment analysis of the miRNA targets. P_NH, NH under HT stress; P_SH, SH under HT stress

Fig. 8

The regulatory miRNA–mRNA interaction pairs in response to HT stress. A Number of four representative regulatory types of miRNA–mRNA interaction pairs. B, C Venn diagrams showing the number of unique and shared miRNAs (B) and their target transcripts (C) in SH compared with NH under HT. D Representative miRNA–mRNA–gene-KEGG regulatory network. Undirected gray lines represent relational pairs, and blue blunt-ended lines represent the inhibitory effect of the miRNAs on the corresponding targets. Yellow star, miRNAs; red vee, target transcripts; green circle, target genes; and purple diamond, KEGG pathways. AP_NH, NH under mild HT stress; AP_SH, SH under mild HT stress; JP_NH, NH under extreme HT stress; JP_SH, SH under extreme HT stress

Fig. 9

Auxin and JA signal transduction pathways are synergistically involved in the regulation of pollen fertility stability under HT stress. A Differentially expressed genes (DEGs) and differential metabolite in the auxin signal transduction pathway. B Relative content of indoleacetic acid. (C, D) Heat maps showing expression levels of DEGs involved in auxin signal transduction. E DEGs and differential metabolite in the JA signal transduction pathway. F Relative contents of JA and methyl jasmonate. G, H Heat maps showing expression levels of DEGs involved in JA signal transduction and indole alkaloid biosynthesis pathways. Green solid circles and red rectangular boxes indicate differential metabolites and DEGs, respectively. AP_NH, NH under mild HT stress; AP_SH, SH under mild HT stress; JP_NH, NH under extreme HT stress; JP_SH, SH under extreme HT stress

Fig. 10

A proposed model showing heat-responsive miRNAs involved in regulating the stability of pollen fertility restoration for CMS-D2 cotton under HT stress. The lines with arrows and blunt ends in the figure signify the promotion and inhibition modes respectively, and the accompanying question marks represent unknown action modes or connections