Next Article in Journal
Treatment with Gold Nanoparticles Using Cudrania tricuspidata Root Extract Induced Downregulation of MMP-2/-9 and PLD1 and Inhibited the Invasiveness of Human U87 Glioblastoma Cells
Next Article in Special Issue
Identification and Characterization of NPR1 and PR1 Homologs in Cymbidium orchids in Response to Multiple Hormones, Salinity and Viral Stresses
Previous Article in Journal
Pyramiding Bacterial Blight Resistance Genes in Tainung82 for Broad-Spectrum Resistance Using Marker-Assisted Selection
Previous Article in Special Issue
Overexpression of a Malus baccata NAC Transcription Factor Gene MbNAC25 Increases Cold and Salinity Tolerance in Arabidopsis
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 21
Issue 4
10.3390/ijms21041280
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Different MicroRNA Families Involved in Regulating High Temperature Stress Response during Cotton (Gossypium hirsutum L.) Anther Development

by
Jin Chen
1,
Ao Pan
1,
Shujun He
2,
Pin Su
3,
Xiaoling Yuan
1,
Shengwei Zhu
4,* and
Zhi Liu
1,*
1
College of Bioscience and Biotechnology, Hunan Agricultural University, Changsha 410128, China
2
Cotton Sciences Research Institute of Hunan Province, Changde 415101, China
3
Hunan Plant Protection Institute, Hunan Academy of Agricultural Sciences, Changsha 410125, China
4
Key laboratory of Plant Molecular Physiology, Institute of Botany, Chinese Academy of Sciences, Beijing 100093, China
*
Authors to whom correspondence should be addressed.
Int. J. Mol. Sci. 2020, 21(4), 1280; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms21041280
Submission received: 2 December 2019 / Revised: 26 January 2020 / Accepted: 12 February 2020 / Published: 14 February 2020
(This article belongs to the Special Issue Abiotic Stress and Gene Networks in Plants 2020)
Browse Figures
Versions Notes

Abstract

:
MicroRNAs (miRNAs) are small molecule RNAs widely involved in responses to plant abiotic stresses. We performed small RNA sequencing of cotton anthers at four developmental stages under normal and high temperature (NT and HT, respectively) conditions to investigate the stress response characteristics of miRNA to HT. A total of 77 miRNAs, including 33 known miRNAs and 44 novel miRNAs, were identified, and 41 and 28 miRNAs were differentially expressed under NT and HT stress conditions, respectively. The sporogenous cell proliferation (SCP), meiotic phase (MP), microspore release period (MRP), and pollen maturity (PM) stages had 10 (including 12 miRNAs), four (including six miRNAs), four (including five miRNAs), and seven (including 11 miRNAs) HT stress-responsive miRNA families, respectively, which were identified after removing the changes in genotype-specific miRNAs under NT condition. Seven miRNA families (miR2949, miR167, and miR160 at the SCP stage; miR156 and miR172 at the MP stage; miR156 at the MRP stage; and miR393 and miR3476 at the PM stage), which had expression abundance of more than 10% of the total expression abundance, served as the main regulators responding to HT stress with positive or negative regulation patterns. These miRNAs orchestrated the expression of the corresponding target genes and led to different responses in the HT-tolerant and the HT-sensitive lines. The results revealed that the HT stress response of miRNAs in cotton anthers were stage-specific and differed with the development of anthers. Our study may enhance the understanding of the response of miRNAs to HT stress in cotton anthers and may clarify the mechanism of plant tolerance to HT stress.

Graphical Abstract

1. Introduction

Cotton (Gossypium hirsutum L.) is one of the most economically important crops because its fiber is a raw material used by the textile industry worldwide. During the whole process of growth and development, the cotton plant is sensitive to temperature, and the ideal temperature range is 20–30 °C. At present, the human population growth and global industrialization are the driving forces of the accelerating global warming. The direct manifestation of climate change is the rise in global temperatures, especially the frequent occurrence of extreme high temperature (HT) during summer, which occurs at the same time as the flowering period of cotton. Long-term HT stress may result in anther indehiscence, microspore abortion, and premature degeneration of the tapetum [1], which eventually leads to boll shedding, threatening the global cotton production. HT stress decreases pollen fertility and crop yield, but the molecular mechanism of the effect of HT stress on plant pollen fertility is still unknown.
MicroRNAs (miRNAs) are a class of small endogenous RNA molecules derived from longer single-stranded precursors (pre-miRNA) [2,3]. The length of mature miRNAs range between 20 and 24 nucleotides (nt). The miRNAs ubiquitously exist in the plant kingdom, and are highly conserved in all kinds of organisms like mosses, flowering monocots, and dicots [4]. miRNAs negatively regulate the expression of target genes at the post-transcriptional regulation by degrading target mRNAs or repressing gene translation [5].
For almost a decade, some miRNAs have been categorized into different groups (miRNA families) on the basis of mature miRNA, seed sequence (the positions 2–8 nt from the miRNA 5ʹ-end) [6], and/or structure of pre-miRNAs [7,8]. A miRNA family is an important unit in the classification of miRNA functions because the sequence and function are conserved among the members of miRNAs within the same family [2,9]. Many plant miRNAs are species-specific, however, those functionally important to plant development and environment adaption are conserved across different species [10]. The species-specific miRNAs are not usually in an organ- or tissue-specific pattern, and tend to evolve neutrally [10,11]. Consequently, these species-specific miRNAs obtained stable and specialized functions [12,13].
miRNAs are involved in the plant’s response to abiotic stresses [4], such as drought, salinity, low temperature (LT), and HT [14,15]. Currently, many miRNAs have been detected in plants, including rice, wheat, barley, and Arabidopsis, as a response to heat stress [16,17,18,19]. miR156 was the first miRNA identified in plants and plays a critical role in plant phase change and heat stress memory [20,21]. Under heat stress conditions, miR156 is highly induced [20,22], and plants overexpressing miR156 exhibit enhanced tolerance to heat stress [20,21]. HT stress inhibits the expression of miR156 in cotton [15], the miR157, which belongs to the same miRNA family as that of miR156, shows the same expression trend. Male infertility is observed in miR157-overexpressing plants, suggesting that the target genes of the same family of miRNAs are functionally identical or similar [6,23], and the induced expression of miR157 partially compensates for the loss of function due to the inhibition of miR156 expression. miR172 is involved in the development of anthers in many species, such as Arabidopsis, cotton, and alfalfa [15,24,25]. miR172 regulates floral organ identity genes, such as AP2 and AP2-like genes (TOE, TOE2, and SMZ) [24], and its expression is considerably inhibited under HT stress, whereas its target gene, TOE, is upregulated [16].
At present, only a few miRNAs have been reported to be involved in cotton anther, and the complete picture of miRNAs or miRNA families and their roles during cotton anther development under HT stress need further investigation. In this study, the miRNA profiles of HT-sensitive and HT-tolerant cotton lines were compared, and the miRNAs and the miRNA family at the whole-genome level were analyzed. Moreover, the HT-responsive miRNA families at different developmental stages of cotton anthers were identified. Results showed that the miRNA family responding to HT stress was significantly different at different developmental stages. The miR160, miR167, and miR2949 families are the main HT stress-responsive families at the sporogenous cell proliferation (SCP) stage. The miR156 and miR172 families showed high abundance ratio at the meiotic phase (MP). Only the miR156 family was dominantly expressed at the microspore release period (MRP). Lastly, miR3476 and miR393 were the major families that responded to HT stress at the pollen maturity (PM) stage. The putative target genes of these major HT stress-responsive miRNA families were predicted from the transcriptome sequencing data, and a model of the miRNA-mediated HT stress response during cotton anther development was further established by combining the mRNA-seq and the small RNA (sRNA)-seq information. The identified HT stress-responsive miRNAs and their targets can help illuminate the mechanisms of HT stress tolerance in cotton.

2. Results

2.1. Global Analysis of sRNA Sequencing

Based on the statistical analysis of sequencing results, the 16.8 M clean reads can be measured for each sample on average, and the average ratio of the mapping rate was above 96% (Table S1). The sequences from the anthers constituted four groups according to the four developmental stages (Figure S1), the early stages (SCP and MP) and the late stages (MRP and PM) were clustered together, respectively (Figure S2). The clean and unique reads in the early stages of anther development were 148.2 and 51.56 M, whereas those in the late stages of anther development were 120.4 and 30.81 M, respectively (Figure S3). The read length statistics revealed that all samples showed the highest abundance at 24 nt, followed by 21 nt (Figure S4a). Usually, a 24 nt sRNA is considered a small interfering RNA (siRNA) involved in the RNA-directed DNA methylation [26,27]. The change in the ratio of 24 to 21 nt (24 nt/21 nt) may reflect the degree of DNA methylation. Results showed that the 24 nt/21 nt in the late stages was remarkably lower than in the early stages. T4 had the lowest ratio (2.6), and H1 had the highest ratio (2.9 times that of T4) (Figure S4b).
Clean reads were aligned against the miRBase. Novel miRNAs were predicted for sequences that lacked homology to known miRNAs and the “-novel” was added in the standard name to distinguish from known miRNAs. A total of 77 miRNAs, including 33 known miRNAs and 44 novel miRNAs (seven novel miRNAs cannot be uniformly named), were identified in all samples from the four anther developmental stages (miRNAs with TPM (transcripts per million) < 10 in all samples were removed) (Table S2). The number of miRNAs detected in the early stages of anther development was more than that in the late stages (Figure S5).
All 70 miRNAs (33 known and 37 novel miRNAs) were classified into 41 miRNA families. Among these, 13 families contained two or more members, and the rest were single-member families. The miR482 family contained seven members, followed by the miR3476 and miR156 families, each including five members (Table S2). The seven unnamed miRNAs belonged to seven upland cotton-specific miRNA families, indicating that species-specific miRNAs also exist under the premise of high homology of miRNAs.

2.2. Identification and Expression Abundance Analysis of Differentially Expressed miRNAs

A total of 68 miRNAs (44 upregulated and 24 downregulated) were differentially expressed under normal temperature (NT) condition (Figure 1b). While comparing the HT-tolerant line and the HT-sensitive line at each anther developmental stage, only 38 differentially expressed miRNAs (21 upregulated and 17 downregulated) were found under HT stress (Figure 1a). Redundant miRNAs that appeared at different stages were removed, and 41 and 28 unique differentially expressed miRNAs were identified under NT and HT conditions, respectively.
According to the expression level, 6, 8, and 14 differentially expressed miRNAs had high, moderate, and low abundance, respectively, under HT stress considering all the developmental stages (Figure 2a). All 28 differentially expressed miRNAs under HT stress belonged to 17 miRNA families, which contained six multimember miRNA families (miR156, miR172, miR2949, miR393, miR482, and miR3476). The interesting point is that each multimember family had only one member with high abundance (TPM ≥ 500), and the others had low abundance (TPM < 500) (Figure 2b). Two types of miRNAs generally showed opposite expression trends in the anthers of the HT-sensitive and the HT-tolerant lines under HT stress (Figure 2b, Table S3). Additionally, the single-member families miR160 and miR167 had high expression abundance, and the remaining miRNAs (miR535, miR394, miR396, miR827, miR1530, miR4344, miR6960, miR7486, and miR7495) had low expression abundance.

2.3. Various miRNA Families Responding to HT Stress at Different Developmental Stages

The HT stress-responsive miRNA families at different developmental stages were obtained from the comparative analysis of differentially expressed miRNAs in the HT-tolerant and the HT-sensitive lines under HT treatment after removing the changes in genotype-specific miRNAs under NT condition (Figure 3 and Figure S6). The SCP, MP, MRP, and PM stages had 10 (including 12 miRNAs), four (including six miRNAs), four (including five miRNAs), and seven (including 11 miRNAs) HT stress-responsive miRNA families, respectively (Figure 3). This result illustrated that the SCP stage had more miRNAs or miRNA families that responded to HT stress than the other three stages.
At the SCP stage, the multimember miR172 and miR2949 families showed significantly differential expression patterns in the HT-tolerant and the HT-sensitive lines under HT stress. The high abundance of family members in the two families (ghr-miR172 and ghr-miR2949a-5p) were downregulated in the anther of the HT-tolerant line, whereas the other family members (ghr-novel-miR172b-5p and ghr-miR2949b > ghr-miR2949c) showed opposite expression trends (Table S3). The single-member families, namely, miR160, miR396, miR7495, miR4344, and miR167, were downregulated in the HT-tolerant line after suffering from HT stress, whereas the other single-member families miR156, miR482, and miR3476 were upregulated.
The miR156 family responded to HT stress at the MP and the MRP stages. The expression trends of ghr-novel-miR156e-3p > ghr-novel-miR156f-3p and ghr-miR156a > ghr-miR156b > ghr-miR156d increased and decreased, respectively, at the MP stage. The ghr-novel-miR156e-5p > ghr-novel-miR156f-5p was remarkably induced by HT stress at the MRP stage. The expression trend of ghr-novel-miR156g-3p was suppressed in the HT-tolerant line. The ghr-novel-miR156g-5p, a differentially expressed miRNA shared by the MP and the MRP stages, showed opposite expression trends in the two periods. At the MP stage, ghr-novel-miR156g-5p was substantially induced in the HT-tolerant line, but its expression was inhibited by HT stress at the MRP stage (Table S3). Other families (miR172, miR535, miR827, and miR7486) were suppressed in the HT-tolerant line except for miR2949, which was upregulated at the MP stage.
At the PM stage, the members of the miR3476, miR482, and miR393 families also showed opposite expression patterns in the HT-sensitive and the HT-tolerant lines after HT treatment. ghr-miR393, ghr-miR2948-5p, and ghr-novel-miR3476d were upregulated in the HT-tolerant line, whereas ghr-novel-miR393b-3p > ghr-novel-miR393c-3p, ghr-novel-miR482d-3p, ghr-novel-miR3476b-5p, ghr-novel-miR3476b-3p > ghr-novel-miR3476c-3p, and ghr-novel-miR3476c-5p > ghr-novel-miR3476e were downregulated. Four families (miR394, miR156, miR7495, and miR6960) in the anthers of the HT-tolerant line were induced by HT stress (Table S3).
The identified miRNA families and their diverse expression patterns at different stages implicated the complexity of miRNAs involved in the regulation of the response of anthers to heat stress. Usually, the miRNAs with high abundance are considered the dominant and functional products that respond to stresses, such as HT [28,29]. In our study, seven miRNA families (miR2949, miR167, and miR160 at the SCP stage; miR156 and miR172 at the MP stage; miR156 at the MRP stage; and miR393 and miR3476 at the PM stage) with expression abundance of more than 10% of the total expression abundance served as the main regulators that respond to HT stress with positive or negative regulation patterns (Figure 4a). The expression of the members from these miRNA families at different developmental stages under HT stress were confirmed by qRT-PCR, and the expression trend was highly similar with the results of small RNA sequencing (Figure 5).

2.4. Target Gene Prediction and Functional Enrichment

The main HT stress-responsive miRNA families at different developmental stages were used to predict the target genes to understand the possible biological functions or mechanistic pathways of miRNAs (Table S4). The functions and annotations were analyzed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment (Figure 4b). At the SCP stage, 656 target genes were mainly annotated to cell cycle, fatty acid metabolism, and peroxisomes. The KEGG enrichment results revealed that 562 and 221 predicted target genes at the MP and the MRP stages, respectively, were concentrated on energy metabolism. This finding suggested that energy supply was very important during the early development of anther in response to HT stress. In addition, the genes involved in the metabolism of peroxisomes, arginine, and proline and protein processing at the two critical stages can scavenge the reactive oxygen species (ROS) and denatured proteins generated by HT stress to maintain normal anther development. At the PM stage, the functions of 50 target genes were enriched on cell cycle, plant hormone signal transduction, and p53 signaling pathway, all of which are necessary for next-step pollination and responding to stress stimuli.
All predicted target genes were further searched and aligned against the differentially expressed genes identified from the transcriptome sequencing data at four anther developmental stages in the HT-tolerant and the HT-sensitive lines under HT stress (Table S5). Results showed that ARF6 (Gh_A09G0074 and Gh_D05G3078), ARF8 (Gh_D07G1785), ARF10 (Gh_D03G1293), ARF17 (Gh_D06G0360), and GATA (Gh_D07G1252) were the target genes of miR167, miR160, and miR2949 families at the SCP stage. miR156-target SPLs participated in HT stress response at the MP and MRP stages. However, the different SPL genes at the MP (Gh_A11G0706 and Gh_D01G1503) and the MRP (Gh_A12G1380, Gh_A12G1955, Gh_A11G0706, Gh_D04G1985, and Gh_D13G0874) stages demonstrated diverse expression patterns in the HT-tolerant and the HT-sensitive lines because significantly different expression profiles existed among the three members of the miR156 family. In addition, miR172 was inhibited by HT stress at the MP stage. Its target genes TOE1 (Gh_A03G0292 and Gh_D05G3873) and AP2 (Gh_A02G1495, Gh_A10G0822, and Gh_A11G1795) were dramatically induced in the anthers of the HT-tolerant line. At the PM stage, the target genes TIR1 (Gh_A10G1164, Gh_A13G0392, Gh_D13G0434, Gh_A11G1077, Gh_D08G0477, and Gh_D08G0763) and AFB (Gh_D07G2334) were identified in the miR393 family, and SurE (Gh_A06G1679) was identified in the miR3476 family. The main HT stress-responsive miRNA families and their target genes at different anther developmental stages are summarized in the schematic diagram shown in Figure 6. The expression of 19 predicted genes targeted by the identified HT responsive miRNAs (Figure 5) under HT stress were further confirmed by qRT-PCR (Figure 7), and the expression trends (upregulation or downregulation) of these target genes were highly consistent between the data of transcriptome sequencing and qRT-PCR (Figure 7). The expression and regulation of these miRNA families and target genes fluctuated differently between the HT-tolerant and the HT-sensitive lines during anther development. This fluctuation ultimately led to indehiscent, shrunken, and empty anthers in the HT-sensitive line under HT stress.

3. Discussion

3.1. More miRNAs Participate in the Early Period of Anther Development

The regulation of genomic DNA methylation involves 24 nt siRNAs [26,27]. The changes in the 24 nt/21 nt may reflect changes in the variation of DNA methylation during anther development. The increased frequency of genomic DNA methylation inhibits the expression of miRNAs and mRNAs [27,30]. However, in our study, a high level of 24 nt/21 nt was in parallel with a high number of reads (clean and uniq reads) in the early stages of anthers (Figures S3 and S4b). The same results were detected in miRNAs regardless whether they are HT-tolerant or HT-sensitive lines or in NT or HT conditions (Figure S5). More miRNAs during the early stage of anthers may indicate that a strong post-transcriptional regulation occurred, which may involve nucleic acid and protein syntheses and meiosis [5] These results showed that the early stage of anther development may be an important period for miRNAs to participate in the regulation of anther development.

3.2. The Expression of High- and Low-Abundance miRNAs Changes Diversely in Response to HT Stress

The number of differentially expressed miRNAs under NT condition was remarkably higher than that under HT stress, but no difference was found in the categories of miRNAs detected at both temperature conditions. It is speculated that miRNAs have achieved HT stress response by regulating the expression abundance (upregulation or downregulation) of existing miRNAs in cells rather than inducing the expression of new types of miRNA [31].
Each multimember family, which was derived from a family classification of differentially expressed miRNAs under HT stress, contained one member with high abundance and one or more members with low abundance (Figure 2b). The expression trends of these two types of miRNAs were opposite in the HT-sensitive and the HT-tolerant lines (Figure 2b, Table S3). The miRNAs within the same family have a similar function attributable to their conserved sequence and structural configuration [6,23]. High- and low-abundance miRNAs may play major and minor regulatory roles, respectively [28,29]. The coordinated regulation of these two types of abundance of miRNAs ultimately leads to the regulation of heat stress response [28].

3.3. Different MiRNA Families Temporally and Elaborately Regulate Anther Development under HT Stress

The differentially expressed miRNAs in the HT-sensitive and the HT-tolerant cultivars showed a period of specificity in family classification under HT stress (Figure 3). This discovery implied that the different stages of anther development were regulated by different miRNA families.
The miR160, miR167, and miR2949 families showed high abundance ratio expressions at the SCP stage under HT stress (Figure 4a). The miR160 and miR167 families are involved in HT stress response by regulating their target gene ARF [16]. Overexpressing miR160 increased cotton sensitivity to HT stress, as shown by anther indehiscence, which was associated with the suppression of ARF10 and ARF17 expressions [15]. The KEGG enrichment of target genes showed that the three high-abundance miRNA families were mainly involved in the energy metabolic pathways, such as glycolysis, gluconeogenesis, and fatty acid synthesis (Figure 4b). In addition, some target genes were enriched in diterpenoid biosynthesis and in peroxisome and cell cycle signaling pathways (Figure 4b). At the SCP stage, sporulated cells proliferate in the anthers while the anthers are preserving large amounts of energy [5]. In this strategy, high-abundance miRNAs in the HT-tolerant line, including the members of the miR160, miR167, and miR2949 families (ghr-novel-miR160b-5p > ghr-novel-miR160c-5p, ghr-miR167a > ghr-miR167b, and ghr-miR2949a-5p) were downregulated (Figure 5), contributed to the accumulation of energy substances in cells, and ensured the proliferation of pollen mother cells. However, certain concentrations of secondary metabolites and fatty acids can increase the tolerance of cells to HT [32,33]. The SCP anthers suffered from HT stress, which can disturb cell energy metabolism and redox balance in the HT-sensitive line and eventually lead to male sterility.
The miR156 family had high abundance ratios during the MP and the MRP periods. Furthermore, the miR172 family exhibited high abundance expression during the MP period (Figure 4a). miR156 is involved in cell division and differentiation regulation in the early stage of anther development [34]. miRNA156 can be significantly induced by HT stress, leading to male sterility [21,22,35]. The overexpression of miR156 in Arabidopsis causes ambient temperature-sensitive flowering [22]. In this study, ghr-novel-miR156g-3p and ghr-novel-miR156g-5p were significantly induced in the MRP stage of the HT-sensitive line (Figure 5), which may explain the phenotypic changes in the T-line under HT stress. The miR172 family is involved in multiple species, such as Arabidopsis, cotton, and alfalfa, and well known for its response to HT stress [15,24,25]. Under HT stress, the expression of miRNA172 is inhibited [35,36,37]. The conversion of sepals to petals into carpels and reduction of stamen number were achieved by overexpression of miRNA172 [24]. In this study, ghr-miR172 was significantly downregulated in the anther of the HT-tolerant line (Figure 5).
The target gene prediction and the functional enrichment analysis showed that the target genes of high-abundance miRNA families were involved in energy metabolism at the MP and the MRP stages (Figure 4b). Moreover, a part of the target genes was involved in peroxisome synthesis, arginine and proline metabolism, and biosynthesis of secondary metabolites at the MP and the MRP stages (Figure 4b). The MP and the MRP stages are meiosis and microspore release periods, respectively, during which cells need to consume large amounts of energy to ensure normal meiosis and maturation of pollen grains [5]. However, vigorous respiration produces a large amount of ROS, and the synthesis of peroxisomes, arginine, and proline can effectively remove excess ROS in cells and prevent oxidative damage [38]. Meiosis is the period that is most sensitive to heat stress [39,40,41]. The accumulation of secondary metabolites in cells helps increase the heat tolerance of plants under HT stress [33,42]. Heat stress tolerance is a comprehensive process that requires multiple factors to achieve a dynamic balance. Possibly, the abnormal expression of the miR156 and miR172 families disrupted this balance and eventually caused the T-line to exhibit HT sensitivity.
At the last stage of anther development, the miR393 and miR3476 families showed high abundant expression at the same time (Figure 3). The members of the two families (ghr-miR393 and ghr-novel-miR3476d) were upregulated in the HT-tolerant line (Figure 5). miR393 can be induced by a variety of abiotic stresses, such as HT or LT and drought [35,43,44]. Currently, the function of the miR3476 family has not been reported, suggesting that the response mechanism of anthers to HT stress may exceed our current knowledge.
The functional enrichment results showed that the target genes of miR393 and miR3476 families were involved in cell cycle, plant hormone signal transduction, and the p53 signaling pathways (Figure 4b). The cell cycle pathway was presumed to be related to the transformation of mononuclear pollen grains to binuclear pollen grains. HT stress can considerably affect the expression of p53 pathway genes [45], which are involved in a variety of cellular functions, such as apoptosis, growth inhibition, and cellular stress aging [45,46]. Furthermore, endogenous phytohormones, such as auxin, gibberellin, and jasmonate (JA), have been reported to respond to HT stress during anther development [15,47,48].

3.4. Various MiRNA Families and Their Target Genes Generate Cotton Male Sterility under HT Stress

The miRNAs participate in regulation of plant reproduction development and abiotic stress responses [49]. miR160 and miR167 are responsive to HT stress through their target genes (ARFs), which are involved in the auxin signaling pathway [18,19]. miR172 is a well-studied miRNA [50,51] that targets AP2-like genes, such as TOE1, TOE2, and SMZ [50]. These target genes are important transcription factors during the differentiation of floral organs [52]. TOE1 interacts with a subset of JAZ proteins and promotes the biosynthesis of JA, resulting in male sterility in Arabidopsis [53,54]. Multiple AP2 genes were upregulated in the HT-tolerant line, especially Gh_D05G3873, Gh_A03G0292, Gh_A02G1495, and Gh_A11G1795 (Figure 6 and Figure 7). HT stress may have affected anther development by regulating the expression of miR172.
SPLs, as target genes of miR156, are essential for anther development and male fertility in many plants [55,56]. The upregulation of miR156 in HT-sensitive line suppressed the expression of SPLs, which increased the sensitivity of anthers to HT stress. In Arabidopsis, miR156-mediated SPL downregulation has increased the plant’s response to environmental stresses [20]. In this study, two differentially expressed SPLs (Gh_A11G0706 and Gh_D01G1503) were downregulated in the HT-tolerant line at the MP stage (Figure 6 and Figure 7). At the MRP stage, five SPLs were differentially expressed (two were upregulated, and three were downregulated), and Gh_A11G0706 showed remarkably differential expression in the MP and MRP stages. Interestingly, Gh_A11G0706 was significantly induced in the HT-tolerant line from the MP stage to the MRP stage, whereas the opposite trend was observed in the HT-sensitive line (Figure 6 and Figure 7). However, SPLs act as a transcriptional activator of miR172, which prevents the expression of AP2 [50,51]. HT stress was proposed to inhibit the expression of SPLs in the HT-sensitive line, leading to the accumulation of AP2, which interacted with the JAZ protein. The excessive accumulation of JA ultimately leads to male sterility [53,54].
TIR1 and AFB are important phytohormone auxin receptors involved in the response to HT stress as target genes of miR393 [14,57,58]. At the PM stage, HT stress inhibited the expression of TIR1 and AFB in the HT-sensitive line (Figure 6 and Figure 7) and promoted the excessive accumulation of auxin, which may ultimately lead to male sterility [5].

4. Materials and Methods

4.1. Plant Materials

The cotton (Gossypium hirsutum L.) lines used in this research included an HT-tolerant line HLY11 (H) and an HT-sensitive line TS18 (T) [59]. The cotton seeds were sown and grown in greenhouses with a normal temperature (NT) condition, which was characterized by a 12 h/12 h photoperiod and a 28–34 °C/20–28 °C day/night thermoperiod. For the high-temperature treatment, the temperature condition was set to a 39 ± 2 °C/29 ± 2 °C day/night thermoperiod that lasted for 10 d when the cotton plants were flowering (approximately 60 d after sowing).
The buds at the four developmental stages, namely, SCP, MP, MRP, and PM stages, were harvested under NT and HT conditions [5]. The anther tissue of each sample was isolated, immediately frozen in liquid nitrogen, and stored at −70 °C until use for transcriptomic and sRNA sequencing. After HT stress treatment, the anthers of TS18 were indehiscent, shrunken, and empty, whereas those of HLY11 were split normally with abundant pollen grains. No difference was observed between the two cotton lines under NT condition (Figure S7).

4.2. RNA Extraction and Sequencing

Total RNA was extracted from the materials using the TRIzol reagent kit (Invitrogen, Carlsbad, CA, US) following the manufacturer’s instructions. The RNA yield was determined using the NanoDrop 2100 spectrophotometer (Thermo Scientific, Waltham, MA, USA). Partial RNA for transcriptome sequencing was obtained using the Illumina HiSeqTM2500 (OE biotech company, Shanghai, China). The rest of the RNAs were subjected to gel electrophoresis. The fragments with a certain size were recovered, and 3ʹ- and 5ʹ-end linkers were added. The platform of the Illumina HiSeq XTEN was adopted to obtain the sRNA sequencing (OE biotech company, Shanghai, China). The sequencing data were deposited in the NCBI Sequence Read Archive (http://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/Traces/sra, accessed on 22, August, 2019) with accession number SRP218740.

4.3. Data Analysis

The reference genome and transcriptome link of upland cotton were obtained from ftp://public.genomics.org.cn/BGI/cotton/Gossypium_hirsutum/Gossypium_hirsutum_v1.0.gz, accessed on 1, July, 2018. A comparative analysis of the cotton reference genome was made for clean bases [60]. Clean reads were normalized using the following formula: normalized expression = mapped read count/total reads*1,000,000. The uniform screening conditions for the differential expression of mRNA and miRNA were p ≤ 0.05 and fold change ≤ 0.5 or fold change ≥ 2. The comparison between the sequence in this study and the body sequence of the mature miRNA in miRbase (used as know miRNA) was carried out with the Bowtie software [61]. The miRDeep2 software was adopted to predict the novel miRNA by combining the homologous miRNA sequence of cotton and the RNA secondary structure prediction software like RNA fold [62]. An online prediction website (http://mongrn.org/psRNATarget/, accessed on 18, October, 2019) was used to perform the target gene prediction. The target genes were deciphered via the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis [63]. A heatmap was drawn using the TBtools software (South China Agricultural University, Guangzhou, China) [64].

4.4. Real-Time Quantitative Polymerase Chain Reaction (PCR)

Three biological replicates were obtained from each sample for real-time PCR analysis to detect miRNA and mRNA abundance under NT and HT conditions. The primer pairs (Table S6) used for RT-PCR were designed via the Roche LCPDS2 software and synthesized via the Generay PCR (Generay Biotech, Shanghai, China).
The RT-PCR of mRNA was described in a previous study [65]. The miRNAs quantification was achieved by a two-step reaction method, which consisted of reverse transcription (RT) and PCR. The reaction agents for each RT reaction was 0.5 μg RNA, 2 μL miScript HiSpec buffer, 1 μL Nucleics Mix, and 0.5 μL miScript reverse transcriptase mix (Qiagen, Dusseldor, Germany). The volume for each reaction was 10 μL, and the reactions were performed using the GeneAmp® PCR System 9700 (Applied Biosystems, Foster City, CA, USA) with the following conditions: 60 min at 37 °C, followed by heat inactivation of RT for 5 min at 95 °C. The 10 μL RT reaction mix was then diluted 10 fold in nuclease-free water and stored at −20 °C. The real-time PCR was performed using the Light Cycler® 480 II RT-PCR instrument (Roche, Basel, Swiss) with 10 μL PCR reaction mixture that included 1 μL cDNA, 5 μL 2× QuantiFast® SYBR® Green PCR Master Mix (Qiagen, Dusseldor, Germany), 0.2 μL universal primer (Qiagen, Dusseldor, Germany), 0.2 μL microRNA-specific primer, and 3.6 μL nuclease-free water. Reactions were incubated in a 384-well optical plate (Roche, Swiss) at 95 °C for 5 min, followed by 40 cycles at 95 °C for 10 s and at 60 °C for 30 s. At the end of the PCR cycles, melting curve analysis was performed to validate the specific generation of the expected PCR product.
The actin (NM 001327051) gene was used as reference gene to normalize the expression level of mRNAs, and the 5S genes were used for miRNAs. The expression amounts were calculated using the comparative 2−ΔΔCt method [66].

5. Conclusions

In conclusion, sRNA libraries were established at different developmental stages of cotton anthers under HT stress and NT conditions by using high-throughput sequencing technology. A total of 77 miRNAs, including 33 known miRNAs and 44 novel miRNAs, were identified. A total of 41 and 28 miRNAs were differentially expressed under NT and HT conditions, respectively. After removing the changes in genotype-specific miRNAs under NT condition, SCP, MP, MRP, and PM stages had 10 (including 12 miRNAs), four (including six miRNAs), four (including five miRNAs), and seven (including 11 miRNAs) HT stress-responsive miRNA families, respectively. Seven miRNA families with expression abundance of more than 10% of the total expression abundance at the four stages of cotton anther development were the main and important regulators that respond to HT stress with positive or negative regulation patterns. Overall, our results will enhance the understanding of the response of miRNAs to HT stress in cotton anthers.

Supplementary Materials

Supplementary materials can be found at https://0-www-mdpi-com.brum.beds.ac.uk/1422-0067/21/4/1280/s1. Figure S1. Principal component analysis (PCA) of clean reads in anthers at four developmental stages under NT and HT conditions. The 16 samples include the anthers from the HT-tolerant line (H) and the HT-sensitive line (T) under NT (controls) condition at SCP (HC1 and TC1), MP (HC2 and TC2), MRP (HC3 and TC3), and PM (HC4 and TC4) stages and anthers from the HT-tolerant line and the HT-sensitive line under HT stress at SCP (H1 and T1), MP (H2 and T2), MRP (H3 and T3), and PM (H4 and T4) stages. Figure S2. Cluster analysis of clean reads in anthers at four developmental stages under NT and HT conditions. Figure S3. Statistics of clean and unique reads in the early and late stages of anther development. Figure S4. Length distribution and 24 nt/21 nt sRNAs. (a) Length distribution of 18–26 nt sRNAs in 16 samples. (b) The 24 nt/21 nt sRNAs in 16 samples. Figure S5. Total number of miRNAs in the early and late stages of anther development in the HT-tolerant and the HT-sensitive lines under NT (HC and TC) and HT (H and T) conditions, respectively. Figure S6. Family of differentially expressed miRNAs at four anther developmental stages in the HT-tolerant and the HT-sensitive lines under HT (a) and NT (b) conditions. The red marked miRNA families, including miR1530 at the MP stage, miR3476 and one member of miR156 at the MRP stage, and miR827 at the PM stage, presented the same expression patterns and were considered genotype-specific miRNAs. Figure S7. Morphological observation of cotton anther at the dehiscence period under NT and HT conditions. Anthers from the HT-tolerant (HLY11) and the HT-sensitive (TS18) lines under NT (a and b, respectively) and HT (c and d, respectively) conditions. Table S1. Small RNA sequencing data and miRNA statistics in 16 samples, including the anthers from the HT-tolerant line (H) and the HT-sensitive line (T) under NT (controls) condition at SCP (HC1 and TC1), MP (HC2 and TC2), MRP (HC3 and TC3), and PM (HC4 and TC4) stages and anthers from the HT-tolerant line and the HT-sensitive line under HT stress at SCP (H1 and T1), MP (H2 and T2), MRP (H3 and T3), and PM (H4 and T4) stages. Table S2. All miRNAs detected by sRNA sequencing (miRNAs with TPM < 10 in all samples were removed), miRNA sequences, expression levels (TPM) in each sample, standard name of known and novel miRNAs, and family classification. Table S3. Differentially expressed miRNAs under HT stress. Table S4. List of predicted target genes. Table S5. List of differentially expressed genes. Table S6. List of primers for qRT-PCR.

Author Contributions

Z.L. and S.Z. conceived, organized and planned the research. J.C. drafted the manuscript. J.C., P.S. and A.P. analyzed the data of transcriptomic and microRNA sequencing. J.C., S.H. and X.Y. conducted PCR experiments and HT stress treatments. Z.L. and S.Z. modified the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (Nos. 31771847 and 31701764), and National Key R&D Program of China (2018YFD0100401).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Paupiere, M.J.; Van Heusden, A.W.; Bovy, A.G. The metabolic basis of pollen thermo-tolerance: Perspectives for breeding. Metabolites 2014, 48, 89–920. [Google Scholar]
  2. Kaczkowski, B.; Torarinsson, E.; Reiche, K.; Havgaard, J.H.; Stadler, P.F.; Gorodkin, J. Structural profiles of human miRNA families from pairwise clustering. Bioinformatics 2009, 252, 91–294. [Google Scholar] [CrossRef] [Green Version]
  3. Hertel, J.; Stadler, P.F. Hairpins in a Haystack: Recognizing microRNA precursors in comparative genomics data. Bioinformatics 2006, 22, e197–e202. [Google Scholar] [CrossRef] [PubMed]
  4. Zhang, B.H. MicroRNA: A new target for improving plant tolerance to abiotic stress. J. Exp. Bot. 2015, 661, 749–1761. [Google Scholar] [CrossRef] [PubMed]
  5. Chen, J.; Su, P.; Chen, P.Y.; Li, Q.; Yuan, X.L.; Liu, Z. Insights into the cotton anther development through association analysis of transcriptomic and small RNA sequencing. BMC Plant Biol. 2018, 181, 54–165. [Google Scholar] [CrossRef] [PubMed]
  6. Kamanu, T.K.K.; Radovanovic, A.; Archer, J.A.C.; Bajic, V.B. Exploration of miRNA families for hypotheses generation. Sci. Rep. 2013, 3, 2940. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Kozomara, A.; Griffiths-Jones, S. miRBase: Integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 2011, 39, 152–157. [Google Scholar] [CrossRef] [Green Version]
  8. Ding, J.; Zhou, S.; Guan, J. miRFam: An effective automatic miRNA classification method based on n-grams and a multiclass SVM. BMC Bioinform. 2011, 12, 216. [Google Scholar] [CrossRef] [Green Version]
  9. Sun, G.L. MicroRNAs and their diverse functions in plants. Plant Mol. Biol. 2011, 80, 17–36. [Google Scholar] [CrossRef]
  10. Cuperus, J.; Fahlgren, N.; Carrington, J.C. Evolution and functional diversification of MIRNA genes. Plant Cell 2011, 23, 431–442. [Google Scholar] [CrossRef] [Green Version]
  11. Budak, H.; Akpinar, B.A. Plant miRNAs: Biogenesis, organization and origins. Funct. Integr. Genomics 2015, 15, 523–531. [Google Scholar]
  12. Nozawa, M.; Miura, S.; Nei, M. Origins and evolution of microRNA genes in plant species. Genome Biol. Evol. 2012, 4, 230–239. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Taylor, R.S.; Tarver, J.E.; Hiscock, S.J.; Donoghue, P.C. Evolutionary history of plant microRNAs. Trends Plant Sci. 2014, 19, 175–182. [Google Scholar] [CrossRef] [PubMed]
  14. Sunkar, R.; Zhu, J.K. Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell 2004, 16, 2001–2019. [Google Scholar] [CrossRef] [Green Version]
  15. Ding, Y.H.; Ma, Y.Z.; Liu, N.; Xu, J.; Hu, Q.; Li, Y.Y.; Wu, Y.L.; Xie, S.; Zhu, L.F.; Min, L.; et al. microRNA involved in auxin signaling modulate male sterility under high temperature stress in cotton (Gossypium hirsutum). Plant J. 2017, 91, 977–994. [Google Scholar] [CrossRef] [Green Version]
  16. Zhao, J.G.; He, Q.S.; Chen, G.; Wang, L.; Jin, B. Regulation of non-coding RNAs in heat stress responses of plants. Front. Plant Sci. 2016, 7, 01213. [Google Scholar] [CrossRef] [Green Version]
  17. Wang, Y.; Sun, F.; Cao, H.; Peng, H.R.; Ni, Z.F.; Sun, Q.X.; Yao, Y. TamiR159 directed wheat TaGAMYB cleavage and its involvement in anther development and heat response. PLoS ONE 2012, 7, e48445. [Google Scholar] [CrossRef] [Green Version]
  18. Kruszka, K.; Pacak, A.; Swida-Barteczka, A.; Nuc, P.; Alaba, S.; Wroblewska, Z.; Karlowski, W.; Jarmolowski, A.; Szweykowska-Kulinska, Z. Transcriptionally and post-transcriptionally regulated microRNAs in heat stress response in barley. J. Exp. Bot. 2014, 65, 6123–6135. [Google Scholar] [CrossRef] [Green Version]
  19. Li, S.X.; Liu, J.X.; Liu, Z.Y.; Li, X.R.; Wu, F.J.; He, Y.K. Heat-induced tas1 target1 mediates thermotolerance via heat stress transcription factor A1a-directed pathways in Arabidopsis. Plant Cell 2014, 26, 1764–1780. [Google Scholar] [CrossRef] [Green Version]
  20. Stief, A.; Altmann, S.; Hoffmann, K.; Pant, B.D.; Scheible, W.R.; Baurle, L. Arabidopsis miR156 regulates tolerance to recurring environmental stress through SPL transcription factors. Plant Cell 2014, 26, 1792–1807. [Google Scholar] [CrossRef] [Green Version]
  21. Cui, L.G.; Shan, J.X.; Shi, M.; Gao, J.P.; Lin, H.X. The miR156-SPL9-DFR pathway coordinates the relationship between development and abiotic stress tolerance in plants. Plant J. 2014, 80, 1108–1117. [Google Scholar] [CrossRef] [PubMed]
  22. Kim, J.J.; Lee, J.H.; Kim, W.; Jung, H.S.; Huijser, P.; Ahn, J.H. The microRNA156-SQUAMOSA PROMOTER BINDING PROTEIN-LIKE3 module regulates ambient temperature-responsive flowering via FLOWERING LOCUS T in Arabidopsis. Plant Physiol. 2012, 159, 461–478. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Tang, J.Y.; Chu, C.C. MicroRNAs in crop improvement fine-tuners for complex traits. Nat. Plants 2017, 3, 17077. [Google Scholar] [CrossRef] [PubMed]
  24. Shi, J.; Chen, L.; Zheng, R.; Guan, C.P.; Wang, Y.J.; Liang, W.Y.; Yang, S.J.; Wang, L.J.; Gong, L.; Zheng, G.B.; et al. Comparative phenotype and microRNAome in developing anthers of wild-type and male-sterile Lycium barbarum L. Plant Sci. 2018, 274, 349–359. [Google Scholar] [CrossRef] [PubMed]
  25. Wu, G.; Park, M.Y.; Conway, S.R.; Wang, J.W.; Weigel, D.; Poethig, R.S. The sequential action of miR156 and miR172 regulates developmental timing in Arabidopsis. Cell 2009, 138, 750–759. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Ci, D.; Song, Y.P.; Tian, M.; Zhang, D.Q. Methylation of miRNA genes in the response to temperature stress in Populus simonii. Front Plant Sci. 2015, 6, 921. [Google Scholar] [CrossRef] [Green Version]
  27. Song, Y.P.; Ci, D.; Tian, M.; Zhang, D.Q. Stable methylation of a non-coding RNA gene regulates gene expression in response to abiotic stress in Populus simonii. J. Exp. Bot. 2016, 67, 1477–1492. [Google Scholar] [CrossRef] [Green Version]
  28. Ravichandran, S.; Raqupathy, R.; Edwards, T.; Domaratzki, M.; Cloutier, S. MicroRNA-guided regulation of heat stress response in wheat. BMC Genom. 2019, 20, 488. [Google Scholar] [CrossRef] [Green Version]
  29. Hutvagner, G. Small RNA asymmetry in RNA: Function in RISC assembly and gene regulation. FEBS Lett. 2005, 579, 5850–5857. [Google Scholar] [CrossRef] [Green Version]
  30. Matzke, M.A.; Mosher, R.A. RNA-directed DNA methylation: An epigenetic pathway of increasing complexity. Nat. Rev. Genet. 2014, 15, 394–408. [Google Scholar] [CrossRef]
  31. Ragupathy, R.; Ravichandran, S.; Mahdi, M.S.R.; Huang, D.; Reimer, E.; Domaratzki, M.; Cloutier, S. Deep sequencing of wheat sRNA transcriptome reveals distinct temporal expression pattern of miRNAs in response to heat, light and UV. Sci Rep. 2016, 6, 39373. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Wahid, A.; Gelani, S.; Ashraf, M.; Foolad, M.R. Heat tolerance in plants: An overview. Environ. Exp. Bot. 2007, 61, 199–223. [Google Scholar] [CrossRef]
  33. Castander-Olarieta, A.; Montalban, I.A.; Oliveira, E.D.M.; Aversana, E.D.; Amelia, L.D.; Carillo, P.; Steiner, N.; Fraga, H.P.D.F.; Guerra, M.P.; Goicoa, T.; et al. Effect of thermal stress on tissue ultrastructure and metabolite profiles during initiation of radiata pine somatic embryogenesis. Front. Plant Sci. 2019, 9, 2004. [Google Scholar] [CrossRef] [PubMed]
  34. Xing, S.P.; Salinas, M.; Hohmann, S.; Berndtgen, R.; Huijser, P. miR156-targeted and nontargeted SBP-box transcription factors act in concert to secure male fertility in Arabidopsis. Plant Cell 2010, 22, 3935–3950. [Google Scholar] [CrossRef] [Green Version]
  35. Xin, M.M.; Wang, Y.; Yao, Y.Y.; Xie, C.J.; Peng, H.R.; Ni, Z.F.; Sun, Q.X. Diverse set of microRNAs are responsive to powdery mildew infection and heat stress in wheat (Triticum aestivum L.). BMC Plant Biol. 2010, 10, 123. [Google Scholar] [CrossRef] [Green Version]
  36. May, P.; Liao, W.; Wu, Y.; Shuai, B.; McCombie, W.R.; Zhang, M.Q.; Liu, Q.A. The effects of carbon dioxide and temperature on microRNA expression in Arabidopsis development. Nat. Commun. 2013, 4, 2145. [Google Scholar] [CrossRef] [Green Version]
  37. Khaksefidi, R.E.; Mirlohi, S.; Khalaji, F.; Fakhari, Z.; Shiran, B.; Fallahi, H.; Rafiei, F.; Budak, H.; Ebrahimie, E. Differential expression of seven conserved microRNAs in response to abiotic stress and their regulatory network in Helianthus annuus. Front. Plant Sci. 2015, 17, 741. [Google Scholar]
  38. Selinski, J.; Scheibe, R. Pollen tube growth where does the energy come from? Plant Signal. Behav. 2014, 9, e977200. [Google Scholar] [CrossRef] [Green Version]
  39. Song, L.R.; Liu, Z.Q.; Tong, J.H.; Xiao, L.T.; Ma, H.; Zhang, H.Q. Comparative proteomics analysis reveals the mechanism of fertility alternation of thermo-sensitive genic male sterile rice lines under low temperature inducement. Proteomics 2015, 15, 1884–1905. [Google Scholar] [CrossRef]
  40. Morgan, C.H.; Zhang, H.K.; Bomblies, K. Are the effects of elevated temperature on meiotic recombination and thermotolerance linked via the axis and synaptonemal complex? Philos. Trans. R. Soc. Lond. B Biol. Sci. 2017, 372, 20160470. [Google Scholar] [CrossRef] [Green Version]
  41. Draeger, T.; Moore, G. Short periods of high temperature during meiosis prevent normal meiotic progression and reduce grain number in hexaploid wheat (Triticum aestivum L.). Theor. Appl. Genet. 2017, 130, 1785–1800. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Shamloo, M.; Babawale, E.A.; Furtado, A.; Henry, R.J.; Eck, P.K.; Jones, P.J.H. Effects of genotype and temperature on accumulation of plant secondary metabolites in Canadian and Australian wheat grown under controlled environments. Sci. Rep. 2017, 7, 9133. [Google Scholar] [CrossRef] [PubMed]
  43. Zhou, X.F.; Wang, G.D.; Sutoh, K.; Zhu, J.K.; Zhang, W.X. Identification of cold-inducible microRNAs in plants by transcriptome analysis. Biochim. Biophys. Acta. 2008, 1779, 780–788. [Google Scholar] [CrossRef] [PubMed]
  44. Zhao, Y.; Srivastava, D. A developmental view of microRNA function. Trends Biochem. Sci. 2007, 32, 189–197. [Google Scholar] [CrossRef]
  45. Cheng, C.H.; Yang, F.F.; Liao, S.A.; Miao, Y.T.; Ye, C.X.; Wang, A.L.; Tan, J.W.; Chen, X.Y. High temperature induces apoptosis and oxidative stress in pufferfish (Takifugu obscurus) blood cells. J. Therm Biol. 2015, 53, 172–179. [Google Scholar] [CrossRef]
  46. Smardova, J.; Liskova, K.; Ravcukova, B.; Kubiczkova, L.; Sevcikova, S.; Michalek, J.; Svitakova, M.; Vybihal, V.; Kren, L.; Smarda, J. High frequency of temperature-sensitive mutants of p53 in glioblastoma. Pathol. Oncol. Res. 2013, 19, 421–428. [Google Scholar] [CrossRef]
  47. Chhun, T.; Aya, K.; Asano, K.; Yamamoto, E.; Morinaka, Y.; Watanabe, M.; Kitano, H.; Ashikari, M.; Matsuoka, M.; Ueguchi-Tanaka, M. Gibberellin regulates pollen viability and pollen tube growth in rice. Plant Cell 2007, 19, 3876–3888. [Google Scholar] [CrossRef] [Green Version]
  48. Sakata, T.; Oshino, T.; Miura, S.; Tomabechi, M.; Tsunaga, Y.; Higashitani, N.; Miyazawa, Y.; Takahashi, H.; Watanabe, M.; Higashitani, A. Auxins reverse plant male sterility caused by high temperatures. Proc. Natl. Acad. Sci. USA 2010, 107, 8569–8574. [Google Scholar] [CrossRef] [Green Version]
  49. Song, X.W.; Li, Y.; Cao, X.F.; Qi, Y.J. MicroRNAs and their regulatory roles in plant-environment interactions. Annu. Rev. Plant Biol. 2019, 70, 489–525. [Google Scholar] [CrossRef]
  50. Aukerman, M.J.; Sakai, H. Regulation of flowering time and floral organ identity by a MicroRNA and its APETALA2-like target genes. Plant Cell 2003, 15, 2730–2741. [Google Scholar] [CrossRef] [Green Version]
  51. Zhu, Q.H.; Upadhyaya, N.M.; Gubler, F.; Helliwell, C.A. Over-expression of miR172 causes loss of spikelet determinacy and floral organ abnormalities in rice (Oryza sativa). BMC Plant Biol. 2009, 9, 149. [Google Scholar] [CrossRef] [Green Version]
  52. Huang, Z.; Shi, T.; Zheng, B.; Yumul, R.E.; Liu, X.; You, C.; Gao, Z.; Xiao, L.; Chen, X. APETALA2 antagonizes the transcriptional activity of AGAMOUS in regulating floral stem cells in Arabidopsis thaliana. New Phytol. 2017, 215, 1197–1209. [Google Scholar] [CrossRef] [PubMed]
  53. Zhai, Q.; Zhang, X.; Wu, F.; Feng, H.; Deng, L.; Xu, L.; Zhang, M.; Wang, Q.; Li, C. Transcriptional mechanism of jasmonate receptor COI1-mediated delay of flowering time in Arabidopsis. Plant Cell 2015, 27, 2814–2828. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Wang, Y.; Li, Y.; He, S.P.; Gao, Y.; Wang, N.N.; Lu, R.; Li, X.B. A cotton (Gossypium hirsutum) WRKY transcription factor (GhWRKY22) participates in regulating anther/pollen development. Plant Physiol. Biochem. 2019, 141, 231–239. [Google Scholar] [CrossRef] [PubMed]
  55. Wang, Z.S.; Wang, Y.; Kohalmi, S.E.; Amyot, L.; Hannoufa, A. SQUAMOSA PROMOTER BINDING PROTEIN-LIKE 2 controls floral organ development and plant fertility by activating ASYMMETRIC LEAVES 2 in Arabidopsis thaliana. Plant Mol. Biol. 2016, 92, 661–674. [Google Scholar] [CrossRef] [PubMed]
  56. Xing, S.P.; Salinas, M.; Garcia-Molina, A.; Hohmann, S.; Berndtgen, R.; Huijser, P. SPL8 and miR156-targeted SPL genes redundantly regulate Arabidopsis gynoecium differential patterning. Plant J. 2013, 75, 566–577. [Google Scholar] [CrossRef]
  57. Dharmasiri, N.; Dharmasiri, S.; Estelle, M. The F-box protein TIR1 is an auxin receptor. Nature 2005, 435, 441–445. [Google Scholar] [CrossRef]
  58. Budak, H.; Kantar, M.; Bulut, R.; Akpinar, B.A. Stress responsive miRNAs and isomiRs in cereals. Plant Sci. 2015, 235, 1–13. [Google Scholar] [CrossRef]
  59. Liu, Z.; Yuan, Y.L.; Liu, S.Q.; Yu, X.N.; Rao, L.Q. Screening for high-temperature tolerant cotton cultivars by testing in vitro pollen germination, pollen tube growth and boll retention. J. Integr. Plant Biol. 2006, 48, 706–714. [Google Scholar] [CrossRef]
  60. Kim, D.; Langmead, B.; Salzberg, S. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef] [Green Version]
  61. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  62. Friedländer, M.R.; Mackowiak, S.D.; Li, N.; Chen, W.; Rajewsky, N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012, 40, 37–52. [Google Scholar] [CrossRef] [PubMed]
  63. Kanehisa, M.; Araki, M.; Goto, S.; Hattori, M.; Hirakawa, M.; Itoh, M.; Katayama, T.; Kawashima, S.; Okuda, S.; Tokimatsu, T.; et al. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008, 36, D480–D484. [Google Scholar] [CrossRef] [PubMed]
  64. Chen, C.J.; Xia, R.; Chen, H.; He, Y.H. TBtools, a Toolkit for Biologists integrating various HTS-data handling tools with a user-friendly interface. bioRxiv 2018. [Google Scholar] [CrossRef]
  65. Huang, Y.; Mo, Y.J.; Chen, P.Y.; Yuan, X.L.; Meng, F.N.; Zhu, S.W.; Liu, Z. Identification of SET domain-containing proteins in Gossypium raimondii and their response to high temperature stress. Sci. Rep. 2016, 6, 32729. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  66. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−△△CT Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
Ijms 21 01280 g001 550
Figure 1. Histograms of differentially expressed microRNAs (miRNAs) under normal temperature (NT) and high temperature (HT) conditions. (a) The number of differentially expressed miRNAs in the anthers between HT-tolerant line and the HT-sensitive line under HT stress at sporogenous cell proliferation (SCP) (H1 vs. T1), meiotic phase MP (H2 vs. T2), microspore release period MRP (H3 vs. T3), and pollen maturity (PM) (H4 vs. T4) stages. (b) The number of differentially expressed miRNAs in the anthers between HT-tolerant line and the HT-sensitive line under NT condition at SCP (HC1 vs. TC1), MP (HC2 vs. TC2), MRP (HC3 vs. TC3), and PM (HC4 vs. TC4) stages.
Figure 1. Histograms of differentially expressed microRNAs (miRNAs) under normal temperature (NT) and high temperature (HT) conditions. (a) The number of differentially expressed miRNAs in the anthers between HT-tolerant line and the HT-sensitive line under HT stress at sporogenous cell proliferation (SCP) (H1 vs. T1), meiotic phase MP (H2 vs. T2), microspore release period MRP (H3 vs. T3), and pollen maturity (PM) (H4 vs. T4) stages. (b) The number of differentially expressed miRNAs in the anthers between HT-tolerant line and the HT-sensitive line under NT condition at SCP (HC1 vs. TC1), MP (HC2 vs. TC2), MRP (HC3 vs. TC3), and PM (HC4 vs. TC4) stages.
Ijms 21 01280 g001
Ijms 21 01280 g002 550
Figure 2. Differential expression of miRNAs at the four stages of anther development under HT stress. (a) Heatmap of the 28 differentially expressed miRNAs responding to HT stress. miRNAs in red-, black-, and green-dotted boxes have high, medium, and low expression levels, respectively. (b) Heatmap of six major miRNA families responding to HT stress.
Figure 2. Differential expression of miRNAs at the four stages of anther development under HT stress. (a) Heatmap of the 28 differentially expressed miRNAs responding to HT stress. miRNAs in red-, black-, and green-dotted boxes have high, medium, and low expression levels, respectively. (b) Heatmap of six major miRNA families responding to HT stress.
Ijms 21 01280 g002
Ijms 21 01280 g003 550
Figure 3. Anther stage-specific miRNA families responding to HT stress in the HT-tolerant and the HT-sensitive lines.
Figure 3. Anther stage-specific miRNA families responding to HT stress in the HT-tolerant and the HT-sensitive lines.
Ijms 21 01280 g003
Ijms 21 01280 g004 550
Figure 4. Heatmap of the expression abundance of differentially expressed miRNA families and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the predicted target genes. (a) The top five expression abundance rates of miRNA families responding to HT stress in the HT-tolerant and the HT-sensitive lines. The expression abundance rate is determined as the expression abundance (in percentage) of a miRNA family over the total expression abundance of miRNA families at a certain stage (sum of family member’s expression abundance/the sum of expression abundance of all differentially expressed miRNAs at a certain stage). (b) KEGG enrichment of target genes of miRNA families with high expression abundance.
Figure 4. Heatmap of the expression abundance of differentially expressed miRNA families and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the predicted target genes. (a) The top five expression abundance rates of miRNA families responding to HT stress in the HT-tolerant and the HT-sensitive lines. The expression abundance rate is determined as the expression abundance (in percentage) of a miRNA family over the total expression abundance of miRNA families at a certain stage (sum of family member’s expression abundance/the sum of expression abundance of all differentially expressed miRNAs at a certain stage). (b) KEGG enrichment of target genes of miRNA families with high expression abundance.
Ijms 21 01280 g004
Ijms 21 01280 g005 550
Figure 5. qRT-PCR validation of the main HT responsive miRNAs identified by small RNA sequencing. (a) Fold-change of miRNAs at the four developmental stages under HT stress based on small RNA sequencing data. (b) The relative expression of the same miRNA members under HT stress detected by qRT-PCR and represented as means of three replicates (n = 3) ± standard error. Significant differences indicated with letters above each bar (p ≤ 0.05).
Figure 5. qRT-PCR validation of the main HT responsive miRNAs identified by small RNA sequencing. (a) Fold-change of miRNAs at the four developmental stages under HT stress based on small RNA sequencing data. (b) The relative expression of the same miRNA members under HT stress detected by qRT-PCR and represented as means of three replicates (n = 3) ± standard error. Significant differences indicated with letters above each bar (p ≤ 0.05).
Ijms 21 01280 g005
Ijms 21 01280 g006 550
Figure 6. A schematic diagram showing various miRNA families and their target genes involved in cotton male sterility caused by HT stress at different stages of anther development.
Figure 6. A schematic diagram showing various miRNA families and their target genes involved in cotton male sterility caused by HT stress at different stages of anther development.
Ijms 21 01280 g006
Ijms 21 01280 g007 550
Figure 7. qRT-PCR validation of 19 predicted target genes identified by transcriptome sequencing. (a) Fold-change of the predicted genes at the four developmental stages under HT stress based on transcriptome sequencing data. (b) The relative expression of the same genes under HT stress detected by qRT-PCR and represented as means of three replicates (n = 3) ± standard error. Significant differences indicated with letters above each bar (p ≤ 0.05).
Figure 7. qRT-PCR validation of 19 predicted target genes identified by transcriptome sequencing. (a) Fold-change of the predicted genes at the four developmental stages under HT stress based on transcriptome sequencing data. (b) The relative expression of the same genes under HT stress detected by qRT-PCR and represented as means of three replicates (n = 3) ± standard error. Significant differences indicated with letters above each bar (p ≤ 0.05).
Ijms 21 01280 g007

Share and Cite

MDPI and ACS Style

Chen, J.; Pan, A.; He, S.; Su, P.; Yuan, X.; Zhu, S.; Liu, Z. Different MicroRNA Families Involved in Regulating High Temperature Stress Response during Cotton (Gossypium hirsutum L.) Anther Development. Int. J. Mol. Sci. 2020, 21, 1280. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms21041280

AMA Style

Chen J, Pan A, He S, Su P, Yuan X, Zhu S, Liu Z. Different MicroRNA Families Involved in Regulating High Temperature Stress Response during Cotton (Gossypium hirsutum L.) Anther Development. International Journal of Molecular Sciences. 2020; 21(4):1280. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms21041280

Chicago/Turabian Style

Chen, Jin, Ao Pan, Shujun He, Pin Su, Xiaoling Yuan, Shengwei Zhu, and Zhi Liu. 2020. "Different MicroRNA Families Involved in Regulating High Temperature Stress Response during Cotton (Gossypium hirsutum L.) Anther Development" International Journal of Molecular Sciences 21, no. 4: 1280. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms21041280

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop