Phenotypic comparison of male fertility in cotton restorer lines with different cytoplasm under HT stress

Our previous studies have found that the male fertility of cotton restorer line with sterile cytoplasm is susceptible to continuous external HT stress [8]. To further explore the potential mechanism of pollen fertility instability under mild and extreme HT stress, two isonuclear alloplasmic near-isogenic cotton restorer lines with obvious differences in anther phenotypes under external HT stress were used, namely, NH (HT-tolerant) and SH (HT-sensitive) (Fig. 1). At the two ecological spots, AY and JJ, there was generally no obvious difference in the external morphology and size of the intact flowers of NH and SH under mild or extreme HT stress (Fig. 1A, E). However, the anther and pollen fertility of the restorer line SH with sterile cytoplasm was significantly lower than that of normal upland cotton cytoplasm NH, especially under extreme HT stress in JJ spot where the summer temperature is higher (Fig. 1B-D, F–H). In addition, whether in the AY with a mild HT or the JJ with an extreme HT, it was observed that compared with NH, SH showed a significant decrease in the normal anther dehiscence ratio and filament length, and a significant increase in the exposed length of the stigma (Fig. 1I–K). This indicates that the male fertility of the sterile cytoplasmic restorer line SH was significantly affected under HT stress, which is finally manifested in the non-cracking of the anthers and the significantly reduced pollen fertility, whereas the NH performed normally (Fig. 1).

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

Global small RNA sequencing data analysis in cotton pollen grains under HT stress

To explore whether miRNAs are involved in male fertility instability under HT in CMS system cotton, mature pollen grains from NH and SH were collected separately with three biological repeats under mild and extreme HT stress, respectively, and used for small RNA (sRNA) library construction and sequencing to identify changes in miRNA abundance in response to HT. Pearson correlation analysis showed that the correlation coefficients among the three replicates of each sample exceeded 0.95, indicating that the sequencing results were relatively accurate and reliable (Additional file 1: Fig. S1). More than 120 million total raw reads were obtained from the 12 small RNA sequencing (sRNA-seq) libraries, and the average total and unique reads for each sample were approximately 10 and two million, respectively (Fig. 2A, Additional file 2: Table S1). After removing 3′adapter sequences, poly-A tags, small tags less than 18-nt and poor-quality reads from raw data, the clean reads obtained were filtered using mRNA, Rfam and Repbase databases, and the number and percentage of different sRNA reads were classified and counted (Figs. 2, Additional file 1: Figs. S2, S3, Additional file 2: Tables S1, S2). The proportion of valid small RNA (VsRNA) reads in SH was significantly less than that of NH, especially in JJ spot with extreme HT (Fig. 2B, C).

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

Afterward, sRNAs within a certain length range (18–25 nucleotide (nt)) in total and unique reads were selected and counted (Fig. 2D, E, Additional file 2: Table S3). All samples showed the highest abundance at 24-nt, and sRNAs of 21–24 nt in length accounted for more than 66% of all reads in each library, representing the majority of small RNAs. Additionally, no significant difference was found in the abundance of 21 and 24 nt sRNAs between NH and SH under HT (Fig. 2D, E). Considering 24-nt, small interfering RNA (siRNA) is a class of double-stranded small RNA molecules associated with RNA-directed DNA methylation (RdDM) process [35]; the changes in the ratio of 24/21 nt might reflect changes in DNA methylation levels during pollen development. Hence, the ratios of 24/21 nt sRNAs in each sample were compared. The results showed that NH had the highest ratios (3.91 for total and 6.55 for unique reads) and the lowest ratios (2.12 for total and 4.19 for unique reads) under mild and extreme HT, respectively, and the ratios of NH and SH at JJ ecological spot were significantly lower than that of AY (Fig. 2F), so it is speculated that extensive DNA demethylation occurred during pollen development under HT stress [8].

Systematic identification of miRNAs in pollen grains of cotton restorer lines

A total of 404 pre-miRNAs and 459 (358 unique) miRNAs were identified in all samples, including 211 known and 248 newly identified miRNAs (Fig. 3A, Additional file 2: Tables S4, S5). Among the five different categories, the number of pre-miRNAs and unique miRNAs in 'gp2b' was the largest, and 'gp3' was the least, especially in SH (Fig. 3A, Additional file 1: Fig. S4A, B, Additional file 2: Table S5). Almost half of the identified miRNAs were 21-nt in length (228), and 22-nt was the second (93), accounting for about 20.26% (Fig. 3B, Additional file 2: Table S6). Moreover, 104 shared miRNAs were found in all samples, and specific miRNAs accounted for only a small portion (Fig. 3C). Subsequently, the detected miRNAs were analyzed for family and conservation. The results showed that the MIR166 family had the most miRNA members (19), followed by MIR156 and MIR482 with 14 and 13 miRNA members, respectively, and miRNAs were highly conserved among different species (Additional file 1: Fig. S5, Additional file 2: Tables S7–S9).

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

For all identified miRNAs, the first base preference of the shortest 18-nt (57.14%) and longest 25-nt (81.82%) mature miRNAs was mainly ‘G’; mature 19–22 nt miRNAs mostly started with ‘U’ as the first base (37.50–65.59%), while 23 and 24 nt miRNAs mostly started with ‘C’ (42.42 and 45.16%, respectively) (Fig. 3D, Additional file 2: Table S10). The mature miRNAs of the 'gp1_3' categories were found to have the same first base preference as the above, but the first base of the 'gp4' miRNAs was mainly 'U' or 'C' (Additional file 1: Fig. S4C, D, Additional file 2: Tables S4, S10). Furthermore, the base preference of mature miRNAs at different positions was also counted. Clearly, the 5'start and 3'end bases were mainly 'U' and 'C', respectively, whereas there was no obvious preference for the base distribution in most of the middle positions (Fig. 3E, Additional file 1: Fig. S4E, F, Additional file 2: Tables S4, S11).

Expression changes of miRNAs in pollen grains in response to HT stress

To determine the miRNAs that were differentially expressed between NH and SH with different cytoplasm under HT stress, four comparisons were performed and among all 459 expressed miRNAs (200 located in At and 225 in Dt sub-genome), totally 159 differentially expressed miRNAs (DEMs; 126 non-redundant) were identified (Fig. 4, Additional file 2: Table S12), including 74 and 78 that were distributed across the 13 chromosomes of At and Dt sub-genomes of upland cotton, respectively, and seven DEMs on four different scaffolds (Additional file 1: Fig. S6). Remarkably, Chir_D05 (with restorer gene Rf1; eight DEMs out of 26 miRNAs) and its homologous Chir_A05 (10 DEMs out of 27 miRNAs) carried relatively more DEMs than the other chromosomes, which was consistent with our previous sequencing results [8, 36]. Among them, a total of 27 and 12 miRNAs were up- and down-regulated, respectively, in JP_NH versus AP_NH, whereas 34 up- and 54 down-regulated miRNAs were identified in JP_SH versus AP_SH. Under mild HT stress, only 26 up- and six down-regulated miRNAs were identified for the comparison of AP_SH versus AP_NH; comparatively, a total of 72 miRNAs were differentially expressed (35 up- and 37 down-regulated) between NH and SH under extreme HT stress (Fig. 4A). The Venn diagram displayed that the combinations AP_SH versus AP_NH and JP_SH versus JP_NH had only 20 DEMs in common. However, the DEMs only under extreme HT stress (52) accounted for 61.90% of the total DEMs (84) of these two combinations. This indicated that most DEMs between NH and SH had differential expression changes in response to HT (Fig. 4B). The expression levels of more than half of DEMs in all four comparison combinations were moderately expressed, especially under extreme HT stress (Fig. 4C). Hierarchical clustering analysis of genome-wide differential expression levels in all 12 samples is shown in Fig. 4D with a bending heat map. Apparently, most miRNAs between NH and SH were found to have a wide range of expression levels in response to HT stress.

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 log10 (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

Whole-genome identification and expression level analysis of miRNA clusters

Some miRNA genes are distributed in clusters on chromosomes and co-transcribed through a promoter and exist in the form of polycistrons [37]. Based on this, we co-localized all 459 known and newly predicted miRNAs on the genome, aiming to find miRNAs that may be simultaneously transcribed in the form of gene clusters (Additional file 2: Table S4). Totally 39 miRNA clusters or Pre-miRNA Clusters (PmCs) containing 106 expressed miRNAs were identified based on 50-kb cluster spacing, including 18 and 20 that were distributed across the eight and 12 chromosomes of At and Dt sub-genomes of upland cotton, and contained 45 and 55 miRNAs, respectively, and one miRNA cluster PmC39 with six miRNAs on Scaffold897 (Fig. 5A, B). Specifically, there were at least two PmCs distributed on ten chromosomes, namely Ghir_A03 (5), Ghir_A04 (2), Ghir_A05 (2), Ghir_A07 (3), Ghir_A08 (2), Ghir_A13 (2), Ghir_D02 (4), Ghir_D05 (3), Ghir_D11 (3) and Ghir_D12 (2), and the remaining chromosomes only carried one PmC. Among them, 45 miRNAs (39 specific) were differentially expressed, including 21 and 22 located in the At and Dt sub-genomes, respectively, and two DEMs on Scaffold897 (Fig. 5C). Heat map analysis of the expression levels of DEMs in the miRNA clusters above found that most highly expressed miRNAs were significantly induced in SH under extreme HT stress compared with NH, especially four MIR482 and six MIR6300 family miRNAs (Fig. 5D). Interestingly, Chir_D05 with restorer gene Rf1 carried three PmCs containing 11 expressed miRNAs, five of which were differentially expressed. It is worth noting that a miRNA cluster, PmC28, was located in the fine-mapped interval of the Rf1 gene, and contained three expressed miRNAs, of which gra-miR482_L-2R + 2 and gma-miR2118a-3p_R + 1_1ss18TG were differentially expressed (Fig. 5), indicating that some HT-responsive miRNAs may be involved in the process of pollen development and fertility restoration in cotton.

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 Rf1 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 log10 (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

Changes in gene transcript levels are involved in response to HT stress

To explore whether the changes in target gene expression in response to temperature variation in pollen grains, transcriptome profiles were performed on the same materials used for the sRNA-seq, and a total of 6281 DEGs were identified in the above four comparisons (Additional file 2: Table S13). Among them, 340 and 541 genes were up- and down-regulated, respectively, in JP_NH versus AP_NH, while 1547 up-regulated and 1417 down-regulated genes were identified in JP_SH versus AP_SH. Under mild HT stress, only 1143 genes were differentially expressed (1053 up-regulated and 90 down-regulated genes) for the comparison of AP_SH versus AP_NH; comparatively, a total of 1711 up-regulated and 147 down-regulated were identified between NH and SH under extreme HT stress (Fig. 6A). The Venn diagram showed that the combinations AP_SH versus AP_NH and JP_SH versus JP_NH had only 311 DEMs in common. However, the DEGs only under extreme HT stress (1547) accounted for 57.51% of the total DEMs (2690) of these two combinations, indicating that most DEGs between NH and SH had differential expression changes in response to HT (Fig. 6B). Subsequently, hierarchical clustering analysis of genome-wide differential expression levels for the top 100 DEGs in all 12 samples is shown in Fig. 6C with a heat map. Notably, most genes were up-regulated in SH compared to NH under extreme HT stress, but heat shock protein (HSP)-related genes, such as HSP70, HSP22, HSP18.5-C, HSP18.2 and HSP17.3-B (Fig. 6C), showed significantly reduced expression levels, suggesting that a large number of genes in response to HT may be involved in pollen development.

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 log10 (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

To gain insight into the potential role of changes in gene transcript levels in cotton pollen development under HT stress, the GO functional enrichment analysis of all 6281 DEGs was performed. Totally 5986 DEGs were annotated to 3571 functional GO categories, of which 682 were significantly enriched (P-value < 0.05), including 341 biological processes (BPs), 108 cellular components (CCs), and 233 molecular functions (MFs) (Additional file 1: Fig. S7, Additional file 2: Table S14). This mainly included BPs such as ‘response to heat’, ‘response to hydrogen peroxide’, and ‘response to cadmium ion’. Besides, ‘vacuole’, ‘cytosol’, and ‘vacuolar membrane’ were the three most significant enriched CC terms, while MFs related to ‘magnesium ion binding’, ‘phosphoprotein phosphatase activity’, and ‘UDP-glucuronate decarboxylase activity’ were the three important functional GO categories that were involved (Fig. 6D, Additional file 2: Table S14). To better understand biological functions and regulatory networks, a KEGG pathway enrichment analysis of all 6281 DEGs was also carried out, and only 2553 DEGs were enriched in 137 pathways, of which 28 were significantly enriched (P-value < 0.05) (Additional file 2: Table S15). Among these, ‘Oxidative phosphorylation’ was the most significantly enriched pathway, followed by ‘Citrate cycle (TCA cycle)’, ‘Protein processing in endoplasmic reticulum’, and’Glycolysis/Gluconeogenesis’ (Fig. 6E). Additionally, the biosynthesis or metabolism of amino acid and fatty acid pathways, such as ‘Valine, leucine and isoleucine biosynthesis’, ‘Cysteine and methionine metabolism’, ‘Alanine, aspartate and glutamate metabolism’, ‘Arginine and proline metabolism’, ‘Lysine degradation’ and ‘Fatty acid degradation’, were also significantly enriched (Fig. 6E, Additional file 2: Table S15). Based on these results, it is therefore inferred that the homeostasis of energy metabolism-related pathways among carbohydrates, fatty acids and proteins may play a crucial role in pollen development and fertility restoration under HT stress.

Identification of miRNA targets by degradome sequencing

To determine the target genes of identified heat-responsive miRNAs in pollen grains of cotton CMS-D2 restorer, two degradome sequencing libraries (P_NH and P_SH) for the same pollen samples were constructed to validate the cleavage sites of miRNAs, and more than 15 million raw reads per library were obtained, namely 15,368,314 and 18,907,383 for P_NH and P_SH, respectively. After removing the 3' adaptor and the number of sequences with shorter splicing site tags (< 15-nt), the ratio of unique mappable reads and transcript mapped reads were all larger than 99% and 70%, respectively. Also, the number of covered transcripts ranged from 79,971 to 91,566, which accounted for approximately 55.78% and 63.87% of the total number of input transcripts for P_NH and P_SH, respectively (Table S16), indicating that the degradome sequencing produced a relatively high coverage of degradation fragments. In total, 196 miRNAs and 5,133 cleavage events involving 3,740 candidate target genes were identified (Fig. 7A, Additional file 2: Table S17), of which 180 miRNAs and 1160 target genes were shared between P_NH and P_SH (Fig. 7B, C). Based on the signature number and abundance at each occupied transcript position, these cleaved transcripts were categorized into five categories: 0, 1, 2, 3 and 4. Obviously, category 2 had the most miRNAs, while category 4 exhibited the most cleavage events and transcripts in these two degradome libraries (Fig. 7D–F). Specifically, category 0 contained 23 and 28 miRNAs corresponding to 33 and 40 cleavage events and transcripts in P_NH and P_SH, respectively, and we used the target plot (T-plot) to display the detailed information of three representative miRNAs cleaving the corresponding target genes (Fig. 7G–I). Besides, the other four categories 1, 2, 3 and 4 were counted, and one miRNA and its cleaved target gene were also selected for T-plot in each category (Fig. 7D–F, Additional file 1: Fig. S8).

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

Based on GO functional enrichment analysis, the targets of the miRNAs described above were significantly enriched in ‘regulation of pollen tube growth’, ‘pollen tube growth’, ‘pollen wall assembly’, ‘pollen germination’, ‘pollen sperm cell differentiation’, and ‘pollen development’ (Fig. 7J, Additional file 2: Table S18). Moreover, KEGG enrichment analysis revealed that 724 of all target genes were enriched into 104 different pathways, of which 16 were significantly enriched (P-value < 0.05), including ‘Protein processing in endoplasmic reticulum’, ‘Spliceosome’, ‘Valine, leucine and isoleucine degradation’, ‘Oxidative phosphorylation’, and ‘Terpenoid backbone biosynthesis’ (Fig. 7K, Additional file 2: Table S19). These results suggest that some miRNAs did act in concert with each other to regulate the expression of genes involved in protein metabolism and mitochondrial energy metabolism by cleaving target mRNAs during pollen development and fertility restoration under HT stress.

Comprehensive analysis of miRNA expression profiles and target genes in response to HT in cotton restorer lines

To further determine the regulatory role of miRNAs during pollen development in response to HT, we integrated combined sRNA, transcriptome, and transcriptome sequencing data to compressively analyze the expression profiles of miRNAs and their corresponding target genes. In total, seven and 51 differently expressed miRNA–target mRNA pairs were identified between SH and NH under mild and extreme HT stress, respectively (Fig. 8A). However, only three up-regulated miRNAs, namely ghr-miR390a, tcc-miR396c, and ptc-miR319a, were shared under either mild or extreme HT stress, but no common target transcripts were found in the two comparative combinations (Fig. 8B, C). Given the inhibitory effects of miRNAs on potential targets in general [38], three up–down and one down–up miRNA–mRNA interaction pairs were identified under mild HT; correspondingly, ten and 21 miRNA–mRNA pairs were found to exhibit opposite up–down and down–up expression patterns under extreme HT, respectively (Fig. 8A, Table 1). Among these 35 negative regulatory pairs, seven miRNAs were found to correspond to multiple targets respectively, and the remaining 16 miRNA–target mRNA pairs were expressed one-to-one (Table 1). Further KEGG enrichment annotation analysis found that these miRNA targets were significantly enriched in ‘Ascorbate and aldarate metabolism’, ‘Oxidative phosphorylation’, ‘Phenylpropanoid biosynthesis’, ‘Plant hormone signal transduction’, and carbohydrate metabolism such as ‘Galactose metabolism’, ‘Other glycan degradation’, ‘Pentose and glucuronate interconversions’ (Additional file 2: Table S20, Fig. 8D).

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

Table 1 The negative regulatory miRNA–mRNA interaction pairs in SH compared with NH under HT stress

Among these targets, a down-regulated target of mtr-miR167a_R + 1, Ghir_A07G018870.1, is a transcription factor encoding auxin response factor 8-like isoform X1 (ARF8) and responsible for the auxin signaling transduction. In contrast, two receptor-like protein kinase (RLK) genes PERK12 and RBK1 encoding proline-rich RLK and receptor-like cytosolic serine/threonine-protein kinase, respectively, were up-regulated in HT-sensitive SH, while ghr-MIR156d-p3 and gra-MIR8643a-p3 targeting them were inhibited under HT. In the ‘Ascorbate and aldarate metabolism’, the target Ghir_A13G022250.1 of tcc-miR167c and the target Ghir_A13G022220.2 of ghr-miR390a are L-ascorbate oxidase homolog encoding pollen-specific protein (Bp10), and two other targets of ghr-miR390a, Ghir_A11G033170.1 (Receptor-like kinase, LIP1) and Ghir_D08G026350.1 (Aldehyde oxidase, GLOX1), were also down-regulated in SH compared with NH under HT stress. Moreover, the target gene Ghir_A09G001320.1 (Pectinesterase 31-like, PME31) of tcc-miR396c_L-1 involved in the ‘Pentose and glucuronate interconversions’ pathway and the target genes Ghir_A01G016050.1 (Beta-galactosidase 13-like precursor, BGAL13) and Ghir_A04G005510.1 (Putative phosphatidylglycerol/phosphatidylinositol transfer protein, DDB_G0282179) of ghr-MIR169b-p3_1ss6AG were also significantly different, indicating that the sugar and lipid metabolism and transport pathways are involved in the regulation of cotton male fertility under HT. Based on the above results, we constructed a comprehensive molecular network of miRNA–mRNA–gene-KEGG involved in the regulation of pollen development in response to HT stress (Fig. 8D).

Excessive auxin accumulation leads to pollen abortion in CMS-D2 restorer line under HT stress by disrupting the homeostasis of JA metabolism

To better demonstrate the involvement of predicted miRNAs and their target genes in pollen fertility stability under HT stress, we also performed metabolite profiling on the same pollen samples [39], and primarily compared the expression levels of related DEGs involved in ‘Plant hormone signal transduction’ pathway (Fig. 9A, E). The relative contents of indoleacetic acid (IAA/auxin) and JA were significantly higher in SH than that in NH, but the active MeJA in SH was obviously down-regulated under HT stress (Fig. 9B, F). Most of the DEGs involved in the auxin signal transduction pathway were significantly up-regulated in SH, among which especially the nine auxin response factor (ARF) genes were more pronounced under extreme HT stress in JJ (Fig. 9C, D), suggesting that heat-activated auxin signal transduction may result in pollen abortion. Similarly, the majority of DEGs involved in the JA signal transduction pathway, especially the two receptor genes COI1, were significantly up-regulated in SH under both mild and extreme HT stress (Fig. 9G). Interestingly, DEGs related to the ‘indole alkaloid synthesis’ pathway downstream of JA signal transduction, including YUC2, YUC6 and AAEs were also up-regulated in SH (Fig. 9H). In summary, we thus preliminarily infer that HT may disturb the dynamic equilibrium of JA and MeJA synthesis, and that induced JA signaling may activate the expression of downstream auxin synthesis-related genes and cause excessive auxin accumulation, followed by a cascade of auxin signal transduction that finally lead to pollen abortion in HT-sensitive restorer line SH under heat stress (Fig. 1).

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