Uncovering Male Fertility Transition Responsive miRNA in a Wheat Photo-Thermosensitive Genic Male Sterile Line by Deep Sequencing and Degradome Analysis

MicroRNAs (miRNAs) are endogenous small RNAs which play important negative regulatory roles at both the transcriptional and post-transcriptional levels in plants. Wheat is the most commonly cultivated plant species worldwide. In this study, RNA-seq analysis was used to examine the expression profiles of miRNA in the spikelets of photo-thermosenisitive genic male sterile (PTGMS) wheat line BS366 during male fertility transition. Through mapping on their corresponding precursors, 917–7,762 novel miRNAs were found in six libraries. Six novel miRNAs were selected for examination of their secondary structures and confirmation by stem-loop RT-PCR. In a differential expression analysis, 20, 22, and 58 known miRNAs exhibited significant differential expression between developmental stages 1 (secondary sporogenous cells had formed), 2 (all cells layers were present and mitosis had ceased), and 3 (meiotic division stage), respectively, of fertile and sterile plants. Some of these differential expressed miRNAs, such as tae-miR156, tae-miR164, tae-miR171, and tae-miR172, were shown to be associated with their targets. These targets were previously reported to be related to pollen development and/or male sterility, indicating that these miRNAs and their targets may be involved in the regulation of male fertility transition in the PTGMS wheat line BS366. Furthermore, target genes of miRNA cleavage sites were validated by degradome sequencing. In this study, a possible signal model for the miRNA-mediated signaling pathway during the process of male fertility transition in the PTGMS wheat line BS366 was developed. This study provides a new perspective for understanding the roles of miRNAs in male fertility in PTGMS lines of wheat.


INTRODUCTION
MicroRNAs (miRNAs) are an extensive class of non-coding small endogenous RNAs approximately 21 nucleotides (nts) in length that are crucial regulators of gene expression in plants (Bartel, 2004). MiRNAs were first described by Lee et al. (1993) in C. elegans (Lee et al., 1993). To date, 46,231 mature miRNAs have been identified in the Sanger miRNA database (miRBase, version 21.0). In plants, these miRNAs are usually transcribed into primary miRNAs (pri-miRNAs) by RNA polymerase (Pol II) (Lee et al., 2004). The pri-miRNA is then processed into a stem-loop precursor miRNA (pre-miRNA) and further processed by Dicerlike 1 (DCL1), or occasionally by DCL4, into a mature miRNA duplex (miRNA/miRNA * ) with 20-22 nt in the cell nucleus (Rajagopalan et al., 2006). Subsequently, the miRNA/miRNA * is exported into the cytoplasm by Hasty (the plant ortholog of exportin 5) (Yang et al., 2006). The miRNA/miRNA * interacts with argonaute (AGO) and is then incorporated into the RNAinduced silencing complex (RISC) which mediates the unwinding of the miRNA/miRNA * (Baumberger and Baulcombe, 2005). The miRNA serves as a guide for recognition of its target transcripts and interacts with its target genes by perfect or near-perfect complementary base pairing, whereas the miRNA * is often degraded (Bartel, 2004;Park et al., 2005;Jung et al., 2009). Finally, RISC suppresses the expression of mRNA targets in a sequence specific manner by degrading target mRNAs or repressing the translation process (Bartel, 2004). Extensive studies have demonstrated that miRNA plays a critical role in many biological processes such as organ senescence (Lim et al., 2010;Pei et al., 2013), flowering time (Wollmann et al., 2010), reproductive development (Peng et al., 2012;Yan et al., 2015), leaf development (Scarpella et al., 2010), fertility transition Ding et al., 2016), and root morphology (Moreno-Risueno et al., 2010). miRNAs are also reported to be associated with plant responses to biotic and abiotic stresses (Leung and Sharp, 2010;Sunkar, 2010;Todesco et al., 2010;Feng et al., 2014;Shu et al., 2016).
Hybrid cultivars have been widely used in plant production to improve yield and have other favorable agronomic traits. For wheat, Triticum aestivum L. (2n = 6x = 42), one of the main grain crops for humans, hybrid cytoplasmic male sterility (CMS) systems or environment-sensitive sterility systems are the basis of hybridization systems in China and India (Longin et al., 2012). Currently, cytoplasmic male sterility lines (CMS) with a maintenance line and a restoration line are mainly used in a three-line breeding system, whereas environment (temperature and photoperiod)-sensitive sterility systems (PTGMS) are mainly used in a two-line breeding system. Plant fertility is controlled by temperature and/or photoperiod in environment-sensitive sterility systems. For example, the critical temperature for anther/pollen (male) fertility transition of PTGMS line BS20 is in the range of 10-13 • C for temperature and 11.5-12 h for light period, during the spikelet differentiation stage. At the day average temperature of 12 • C, the fertility of BS20 increased and the critical temperature decreased with longer day length (Ru et al., 2015). Thus, the PTGMS line can be used as both a maintenance line and a male-sterile line for exploiting heterosis that can improve crop yield, enhance resistance, and improve crop quality (Jordaan, 1996;Virmani and Ilyas-Ahmed, 2001;Yang et al., 2007;Tang et al., 2011;Zhou et al., 2014;Wang et al., 2017).
To date, many miRNAs have been identified in anther development of male sterile lines, such as maize (Shen et al., 2011), cotton (Wei et al., 2013), Brassica juncea (Yang et al., 2013), wheat , soybean , tomato (Omidvar et al., 2015), and rice (Yan et al., 2015). Some of these miRNAs may be involved in plant fertility. Overexpression of osa-miR159 caused male sterility in rice, while mutating the target of osa-miR159 also resulted in abnormalities in anther and pollen development (Kaneko et al., 2004;Tsuji et al., 2006). Recently, osa-miR164d and osa-miR444f were identified and confirmed to be connected to pollen development by interacting with their targets in a rice PTGMS line . In Arabidopsis, miRNA172, interacting with its target genes APETALA2-like, can control flowering time (Chen, 2004). Moreover, overexpression of miRNA167 can lead to male fertility defects (Siré et al., 2009). Furthermore, some studies have reported that MYB33/MYB65 which could be down-regulated by miRNA159a, accelerates anther development in Arabidopsis (Millar and Gubler, 2005).  found that miR167 and transacting small interfering RNA (tasiRNA)-ARF played roles in cold-induced male sterility in wheat . However, few studies have demonstrated the relationship between miRNA and male fertility transition in PTGMS wheat lines. The genetic mechanism is complex and the sterility mechanism is still not clear, both of which are obstacles that restrict its application. In this study, we identified conserved and novel miRNAs which may be involved in male fertility transition during pollen and anther development in PTGMS wheat line BS366 using RNA-seq.

Plant Materials, Growth Conditions, and Sample Collection
The wheat (Triticum aestivum L.) PTGMS line BS366 was used in this study. Plants were grown in soil in plastic pots and vernalized naturally in the field. After the four-leaf stage (when the tip of the fourth leaf had emerged), plants were randomly divided into two groups. One group (fertile condition) was transferred to an artificial climate incubator at 20 • C (with 12-h day/12-h night) and relative humidity of 60-80% for the entire reproductive period. The other group (sterile condition) was held at 10 • C (with 12-h day/12-h night). The plants received cold stress for 5 days. Based on leaf age index and anther length, spikelets were collected from each treatment group at three different developmental stages identified using a dissecting microscope (Olympus SZX12,Tokyo, Japan): stage 1: secondary sporogenous cells had formed (anther length: 1.5 mm); stage 2: all cell layers were present and mitosis had ceased (anther length: 2.2 mm), and stage 3: meiotic division stage (anther length: 3.0 mm) (Tang et al., 2011;. Therefore, samples from these two groups, fertile plants and sterile plants, were designated FS1, FS2, FS3 and SS1, SS2, SS3 for each developmental stage, respectively. Small RNA and degradome sequencing was previously carried out for these six sample sets  and further analyzed here under Sections Small RNA Sequence Data Analysis. For the current work we also collected anther samples from the same stages.

Small RNA Sequence Data Analysis
After filtering the low-quality reads (reads less than 18 nt and reads whose adaptors were null) and removing adaptors/adaptor sequences, clean reads in the range of 18-30 nt were retained for further analysis. rRNAs, tRNAs, small nuclear (sn)RNAs, and small nucleolar (sno)RNAs were identified and eliminated by comparing clean reads to GenBank (http://www.ncbi.nlm. nih.gov) and Rfam (12.0) (http://rfam.xfam.org/) RNA family databases (Gardner et al., 2009). Known miRNAs in the samples were identified by comparison with the specified range in the miRNA database (http://www.mirbase.org). Unannotated small RNA sequences were used for predicting novel miRNAs based on the characteristic hairpin structure of miRNA precursors. Novel miRNAs were predicted using miRdeep2 (https://www.mdcberlin.de/rajewsky/miRDeep). The secondary structures were predicted using the MFOLD3.2 web server (http://mfold.rna.albany.edu/) (Zuker, 2003). Wheat (Triticum aestivum L.) transcripts, downloaded from phytozome v11 (www. phytozome.net), were used to determine potential target mRNA candidates for miRNAs using target gene prediction software psRNATarget with default parameters (http://plantgrn.noble.org/ psRNATarget/) (Dai and Zhao, 2011) according to the rules proposed by Allen et al. (2005) and Schwab et al. (2005). Functional annotations of predicted targets were analyzed using the Blast2GO program (P ≤ 0.05) to obtain the GO annotations from the unigene database using BLASTX hits against the available Nr database in NCBI (Conesa et al., 2005). The KEGG pathway was analyzed using BLASTX (P ≤ 0.05) through the KEGG database.

Degradome Sequence Data Analysis
We used the CleaveLand pipeline to find sliced miRNA targets using wheat unigenes, Viridiplantae miRBase 21.0 mature miRNAs, and miRNAs identified in our study as inputs (Addo-Quaye et al., 2008). To allow for inaccurate target cleavage or variation in miRNA 5 ′ ends, the pipeline was modified to recognize targets cleaved between the 10th and 11th nucleotides of mature miRNA.

RNA Extraction and qRT-PCR Validation of miRNAs and Target Gene Expression
Total RNA was obtained for each anther sample using TRIzol Reagent (Invitrogen, Nottingham, UK) according to the manufacturer's instructions. The purified RNA was dissolved in RNase-free water and stored at −80 • C until subsequent analysis.
RNAs isolated from anthers were transcribed into cDNA using M-MLV Reverse Transcriptase according to the manufacturer's instructions (Takara, Japan). Quantitative real-time PCR (qRT-PCR) was performed with a SYBR Premix Ex Taq TM kit (Takara, Japan) on a Roche LightCycler 480 real-time PCR machine, following the manufacturer's instructions. The protocol used to perform qRT-PCRs was 10 min at 95 • C, followed by 40 cycles of 15 s at 95 • C and 60 s at 61 • C. Ubiquitin was used as a reference gene for miRNA targets analysis, and U6 was used as a reference gene for miRNA analysis. The reactions were performed with three biological replicates and at least two technical replicates per sample. The data were analyzed using the 2 − Ct method (Livak and Schmittgen, 2001), and means ± standard errors (SE) of three biological replicates are presented. The primers for the miRNA qRT-PCR are listed in Table S1.

Overview of Small RNA Sequences Data
Photo-thermosensitive genic male sterile (PTGMS) line wheat had abnormal fertility when treated with low temperature during the meiosis stage (Tang et al., 2011). In this study, the goals were to identify new miRNAs which may be involved in the regulation of wheat male-sterility from a PTGMS wheat line and to illustrate the regulation mechanism of miRNAs in the low temperaturemale fertility transition network. Small RNA sequencing data were generated by Solexa technology and yielded 11.4-13.5 million raw reads of each sample from all six small RNA libraries. Approximately 86% were clean reads and were kept for further analysis after eliminating the unavailable sequences (low quality reads, 3 ′ adapter null, 5 ′ adapter null, insert null, reads shorter than 18 nt, and polyA reads) ( Table 1). The most abundant size classes of sRNA ranged from 18 to 21 nt, and the most abundant size sRNA in all six libraries was 21 nt (45.56-61.89%) (Figure 1). Thus, the sRNA transcriptome in wheat is dominated by 21 nt sRNAs. This is in contrast with other plant species, including Arabidopsis (Fahlgren et al., 2007), rice (Sunkar et al., 2008), peanut (Chi et al., 2011), radish , and cucumber (Martínez et al., 2011;Mao et al., 2012), where 24 nt sRNAs are more abundant.
The clean small RNAs reads were mapped to the wheat genome in order to analyze their expression and distribution in the genome. The total and unique sRNAs were classified into different RNA categories ( Table 2). More than approximately 91% of unique small RNAs corresponding to more than 92% of total small RNAs from the different libraries (except for FS3 25.13%/35.98%) were mapped onto the wheat genome ( Table 2). After removing non-coding RNAs including rRNA, tRNA, snRNA, and snoRNA by blasting sRNA sequences against Genbank and Rfam, we annotated and obtained 2378 (0.06%),

Identification of Known miRNAs
To identify known miRNAs from three pollen developmental stages, the small RNAs mapped on the genome were analyzed with a homology-based method of mapping the unique sequences using miRBase Release 21.0 (July 3, 2014). Currently, approximately 5554 miRNAs derived from 6992 pre-miRNAs are published in miRBase. In our study, a total of 3304 known miRNAs in the spikes of BS366 at the different developmental stages evaluated (1736in SS1, 1705in SS2, 1822in SS3, 1677in FS1, 1699in FS2, and 1980 were identified as orthologs of known miRNAs in other species (Table 3), of which 108 miRNAs were from wheat (  Figure 2). Among these miRNAs, 9 miRNAs in SS1/FS1, 14 miRNAs in SS2/FS2, and 49 miRNAs in SS3/FS3 were up-regulated, and the others were downregulated (Figure 2). Not all differentially expressed miRNAs expressed at different stage. For example, tae-miR171a showed differential expression in FS1/SS1 and FS2/SS2, not in FS3/SS3. However, many differentially expressed miRNAs such as tae-miR167a, tae-miR156 and tae-miR164, showed differentially expressed miRNAs peaked at stage 3 (Figure 2), suggesting that the transition of male fertility may occur at this time. Furthermore, these miRNAs different in abundance at different stages and conditions. For example, tae-miR167, tae-miR5048, tae-miR5062, tae-miR9666, and tae-miR9672 were found over 10,000 times (Table 4). Remarkably, few or no copies of all known miRNAs from wheat were found in the FS3 library ( Table 4 and Table S2), suggesting that expression levels varied markedly among different miRNA families, probably owing to developmental stage-specific expression. This result also indicates that stage 3 (meiotic division) is perhaps the important phase for male fertility transition. In plants, the determination of which strand of the miRNA:miRNA * duplex is incorporated into the RNA induced silencing complex (RISC) is largely based on the identity of the first nucleotide. In Arabidopsis, the AGO family consists of 10 members. AGO1 tends to recruit miRNAs beginning with a uracil (U), while AGO2 and AGO4 display preferences for miRNAs with a 5 ′ terminal adenosine (A), and AGO5 recruits miRNAs initiated with a 5 ′ cytosine (C) (Mi et al., 2008). In this study, the base preference of 21 nt known miRNAs was analyzed. The results revealed that the order of the bias for the first nucleotide of 21 nt known miRNA is A>G>U>C in SS1, SS2, SS3, FS1, and FS2 ( Figure S1). Notably, in FS3, the bias for the first nucleotide of 21 nt miRNA is usually A, suggesting that miRNAs of different stages and different conditions might be recruited by different AGO proteins and enter distinct regulatory pathways.

Identification of Novel miRNAs
Novel miRNAs could be distinguished from other small RNAs by mapping to genomes of target species to predict secondary structures (Zuker, 2003;Yan et al., 2015). Through mapping onto their corresponding precursors, thousands of novel miRNAs were found, ranging from 917 to 7,762 in six libraries (6825 in SS1, 6700 in SS2, 2483 in SS3, 7762 in FS1, 6277 in FS2, and 917 in FS3) ( Table 5); many of these miRNAs had low expression levels (Table S3). To validate the predicted novel miRNAs, 6 novel miRNAs including novel-miR-2, novel-miR-964, novel-miR-4990, novel-miR-2186, novel-miR-5020, and novel-miR-1221, were selected randomly for confirmation by stem-loop RT-PCR. DNA fragments of approximately 60 bp were amplified ( Figure 3B), suggesting that these miRNAs were all expressed in the young panicles of BS366. Moreover, the folding structures for these novel miRNA precursors were predicted using MFOLD (http://mfold.rna.albany.edu/). Six structures of the novel miRNAs with perfect secondary structures are shown in Figure 3A. All of these six novel miRNAs were expressed in both fertile plants and sterile plants ( Figure 3C). Furthermore, all of these novel miRNAs were up-regulated in sterile plants at the three developmental stages ( Figure 3C). The novel miRNAs were named sequentially, according to the chromosome locations, and are described as novel-miR-number.

miRNA Target Prediction and Functional Analysis
Identifying the candidate genes targeted by the miRNAs is crucial to understanding the biological functions of these miRNAs. To further investigate the regulatory effect of miRNAs, the sequences of all known and novel miRNAs were aligned to the wheat genome (downloaded from phytozome v11) to predict the potential targets genes (Allen et al., 2005). We found 3,286 targets for unique known miRNAs (Table S4-1) and 17,016 targets for novel miRNAs (Table S4-2). The targets of differentially expressed miRNAs of wheat were the focus of this study. These predicted targets of known miRNAs, such as miR156, miR159, miR164, miR1120, and miR167, are described as growth-regulating factors, MYB family transcription factors, F-box domain-containing proteins, MADS-box family proteins, and SBP-box gene family members (Table S5). These proteins are known to play various roles during plant growth and development processes. To more thoroughly investigate miRNA functions during male fertility transition in PTGMS wheat, gene ontology analysis of the predicted targets was performed. The most frequent term for "biological process" was "regulation of transcription" [such as the targets of tae-miR1120c-5p (Traes_4BS_35F8C45F6), tae-miR1130b-3p (Traes_2AL_EE1350B36) and tae-miR164 (Traes_2BL_6AEE8AC28)] followed by "transporter" and "auxin-activated signaling pathway" [such as the targets of tae-miR167a (Traes_2AL_A7941CB12)] in three comparison groups (FS1/SS1, FS2/SS2 and FS3/SS3) of targets for known miRNAs (Figure 4 and Table S6). In addition, the terms "regulation of flower development" [such as the targets of tae-miR167a (Traes_2DL_8434C0251), tae-miR1130b-3p (Traes_2DL_DA1A74C0D) and tae-miR156 (Traes_2BS_186EA570A)] and "recognition of pollen" [such as the targets of tae-miR1122b-3p (Traes_1DS_7868656E4.1) and tae-miR5049-3p (Traes_2AL_2A541092D. 2)] were annotated in targets of differentially expressed miRNAs of FS1/SS1and FS3/SS3, respectively, but not in FS2/SS2. The most frequent "molecular function" term was "double−stranded DNA binding" and followed by "DNA binding, " "ATP binding, " and "protein binding" (Figure 4 and Table S7). The most frequent "cellular component" term was "mitochondrion." Furthermore, these targets were subjected to KEGG pathway analysis. Among these enriched pathways, "purine metabolism, " "spliceosome, " and "plant hormone signal transduction" had high frequencies for targets of known miRNAs ( Figure S2). In addition, the "calcium signaling pathway" (ko04020) was enriched by KEGG pathway analysis (Table S6).

Target Identification of miRNAs by Degradome Sequencing
To confirm accuracy of the targets of miRNAs, two degradome libraries were generated from the spikelets of the PTGMS wheat line in different fertility condition (fertile and sterile). A total of 35,971,562 (98.16%) clean reads out of 36,647,546 raw reads and 26,622,338 (97.72%) clean reads out of 27,244,395 raw reads were produced in fertile line (FS) and sterile line (SS) libraries, respectively. Furthermore, 26,400,645 (73.89%) and 18,959,057 (71.21%) unique reads could be matched with the wheat genome in FS and SS libraries, respectively ( Table 6). The FS and SS libraries had 1,276,084 (20.39%) and 7,912,816 (12.64%) specific raw reads, respectively, of which 9,636,103 (46.21%) and 6,323,897 (30.33%) were specific unique reads, respectively ( Table 6). The reads that were mapped to wheat unigenes were subjected to further analysis.
A total of 5,357 targets in FS and 10,613 targets in SS were annotated. Among these targets, 1,062 (9.1%) and 6,316 (54.12%) were sample specific in FS and SS, respectively. A total of 5,357 targets in FS and 10,613 targets in SS were annotated ( Figure 5E). The signature abundance of each target was plotted using absolute read number along the target transcripts (t-plot) (German et al., 2008). Four representative t-plots of selected targets are shown in Figures 5A-D. Degradome sequencing data suggested that targets of a miRNA family belong to the same gene family. Moreover, different genes belonging to the same gene family can also be degraded by different miRNA families. For example, GAMYB (TRAES3BF027700010CFD_t1) was degraded by miR159m-3p, miR159a, and miR159m-5p, suggesting that miRNA regulation is complex. In addition, a total of 899 targets of differentially expressed miRNAs between FS and SS were obtained and analyzed by GO annotations and KEGG pathway analysis to more thoroughly investigate miRNA functions during male fertility transition in PTGMS wheat ( Table S7).
Among these targets, the targets of miR156, miR164, and tae-miR171a were predicted, and the others were confirmed using degradation sequencing. These miRNAs with their targets were examined by qRT-PCR. It was found that the expression of most miRNAs was negatively correlated with their targets (Figure 6) between fertile plants (FS) and sterile plants (SS) during pollen development. In addition, two novel miRNAs were selected to verify the correlation. Both novel-miR-964 and novel-miR-2186 target the pentatricopeptide repeat (PPR) gene family (Traes_7DS_16B6D3A80.2 and Traes_1BL_5BC43A438.1, respectively). qRT-PCR showed that the expression of these two novel miRNAs is also negatively correlated with their targets ( Figure 6B). These negative correlations suggest that miRNAmediated mRNA silencing during pollen fertility transition period probably caused male sterility.

DISCUSSION
Recently, many studies have indicated that some miRNAs regulate anther development and have some relationship to male sterility Yan et al., 2015;Zhang et al., 2016). The relationship between CMS and miRNAs have been studied in some plants such as cotton (Wei et al., 2013), maize (Shen et al., 2011), and rice (Yan et al., 2015). In addition, very recent research proposed a possible model for the control exerted by miRNAs on signaling pathways during fertility transition in a rice PTGMS line . In wheat,  analyzed the expression of miRNAs in TGMS wheat line when responsed to temperature (cold stress) to study the regulation of miRNAs . However, the investigation of miRNAs in PTGMS wheat lines is very limited. Therefore, characterizing the role of miRNAs in PTGMS wheat lines would be extremely useful and could contribute to an improved understanding of the molecular functions of miRNAs during male fertility transition in the male sterile wheat lines used in two-line breeding systems. In this study, we examined miRNA profiles determined by high-throughput sequencing and qRT-PCR to investigate the expression of miRNAs in the PTGMS wheat line BS366 during pollen development (male fertility transition). In order to construct miRNA-gene regulatory network to illustrate the role of miRNA in male fertility transition in the male sterile wheat line, we not only investigated the regulation of miRNAs which responsed temperature during pollen development, but also analyzed the contribution of photoperiod for pollen development. These analyses showed that miRNAs interacting with their targets may play an important role in pollen development. In the present study, a total of 58 wheat-specific known miRNAs were identified in the spike of BS366 as differentially expressed in sterile and fertile conditions ( Table 4).
Targets of these differentially expressed miRNAs include MYB family protein, kinases, RNA-binding proteins, PRR-containing proteins, and stress-resistance related proteins (Tables S5, S6), implying that these proteins may play a particular role in pollen development.
In plants, it is well known that miR164 targets the NAM gene family [CUC1 (CUP-SHAPED COTYLEDON1) and CUC2 (CUP-SHAPED COTYLEDON2)] which are a type of plant-specific transcription factor; they have been shown to participate in various physiological and biochemical processes during plant development and in response to biotic/abiotic stresses (Kikuchi et al., 2000;Fujita et al., 2004;Hu et al., 2006;Takasaki et al., 2010;Tran et al., 2010;Nuruzzaman et al., 2012;. To date, several NAM genes have been cloned and characterized in wheat, and most can be induced by biotic/abiotic stresses (Takasaki et al., 2010;Xia et al., 2010;. miR164 is also involved in age-dependent cell death in Arabidopsis leaves by cleaving ORE1, which is another NAM transcription factor (Kim et al., 2009). In addition, mir164-CUC1 and mir164-CUC2 plants both exhibit leaf shape and polarity defects, extra petals, missing sepals, and reduced fertility (Laufs et al., 2004;Mallory et al., 2004). In addition, miR164 was up-regulated under UV-B radiation treatment in maize, suggesting that the expression of miR164 might be regulated by light (Casati, 2013). In the present study, tae-miR164 was found to be down-regulated in the pollen at later developmental stages (SS2 and SS3) in sterile plants ( Table 4) but not in SS1/FS1, and this finding was verified by qRT-PCR (Figure 6A), indicating that miR164 may participate in fertility regulation during pollen development. The expression of target Traes_2BL_6AEE8AC28.1, which was annotated as the CUP-SHAPED COTYLEDON (CUC)/NO APICAL MERISTEM (NAM) gene family, was examined by qRT-PCR. This target was down-regulated in the sterile group, especially in later anther developmental stages (SS2 and SS3) ( Figure 6A) demonstrating that the target was silenced by miRNA.
Recently, it has been reported that miR156 (Schwab et al., 2005), miR159 (Achard et al., 2004), and miR172 (Aukerman and Sakai, 2003;Chen, 2004 are involved in control of flowering time. miR156, which targets SPL (SQUAMOSA PROMOTER BINDING PROTEIN-LIKE) transcription factors, promotes flowering and plays important regulatory roles throughout the growth and development stages (Cardon et al., 1997;Schmid et al., 2003). Indeed, recent studies have shown that vegetative phase changes in Arabidopsis could be regulated by miR156 which promotes the expression of the juvenile phase and represses the expression of the adult phase . miR172 promotes flowering primarily by repressing a set of AP2-like genes (TOE1, TOE2, and TOE3) (Aukerman and Sakai, 2003), whereas miR172-overexpressing plants exhibit early flowering under both long days (LDs) and short days (SDs). In addition, it has been reported that the expression of miR172 could be mediated by Cryptochrome 1 (CRY1) and Cryptochrome 2 (CRY2) after blue light stimulation involving GIGANTEA (GI)-mediated photoperiodic flowering (Jung et al.,     2007; Zhou et al., 2016). In the present study, the target of miR172 (Traes_1DS_2BB06F54C.1), which was confirmed by degradome sequencing, was negatively correlated with miR172 and down-regulated in SS, whereas miR172 was significantly up-regulated. miR156 has been known to play a crucial role in the development of sporogenic tissues with its targets in rice  and Arabidopsis (Xing et al., 2010). The targets of miR156 (Traes_2BL_186EA570A), annotated as the SBP family (for SQUAMOSA-PROMOTER BINDING PROTEIN) which is a sequence specific DNA-binding domain found in plant proteins, was nearly silenced in SS. miR159 targets MYB transcription factors, playing important roles in flower development (Achard et al., 2004;Jones-Rhoades et al., 2006). In rice, mutations of OsGAMYB resulted in defects in the anthers and pollen development (Kaneko et al., 2004). Moreover, studies have shown that the abundance of miR159 is directly related to fertility in Arabidopsis (Achard et al., 2004) and rice (Tsuji et al., 2006). Furthermore, overproduction of miR159 resulted in late flowering specifically under SDs. In this study, the target of tae-miR159/tae-miR319 (Traes_1AL_041EBB4A6.2) was down-regulated significantly suggesting that tae-miR159 might be involved in regulation of anther development. We also found that there was differential expression of the calmodulin-binding protein ( Table S7). As an important second messenger, calcium ions (Ca 2+ ) participate in a number of physiological responses in both animal and plant cells. The production and transduction of calcium signals are realized through the change in calcium distribution and concentration in cytoplasm and organelles (Dodd et al., 2010). Studies demonstrated that Ca 2+ signals contribute to cryptochrome (CRY) and phytochrome (PHY) signaling (Spalding, 2000;Guo et al., 2001). In addition, calcium concentrations that are inappropriately located may cause pollen development defects, resulting in male sterility in PGMS (photoperiod-sensitive genic male-sterile) rice (Tian et al., 1998). There are significant differences in calcium distribution between fertile and sterile anthers. In rice, some studies indicated that osa-miR1432 and osa-miR812d may be involved in Ca 2+ -mediated signaling pathways through interaction with their target genes (Yan et al., 2015). Osa-miR5976 showed up-regulated expression in rice sterile lines . In our study, tae-aly-miR825-5p was verified as targeting calmodulin binding protein by degradome sequencing (Traes_1BL_F0F995FE3.1), and expression of tae-aly-miR825-5p was up regulated in sterile lines ( Table 4 and Table S6). Transcripts of Traes_1BL_F0F995FE3.1 were positively correlated with the expression of tae-aly-miR825-5p ( Figure 6A). Therefore, the Ca 2+ signal pathway also may be regulated by environmental factors and miRNA, causing male sterility in wheat.
There are many studies demonstrating that miR167 with its targets plays a crucial role in floral organ development. Overexpression of miR167 leads to Arabidopsis male sterility (Siré et al., 2009). Mutation of auxin response factor genes ARF6 and ARF8 (targets of miR167) could cause flower organ developmental defects (Nagpal et al., 2005;Peng et al., 2006). The targets of tae-miR1120 and tae-miR1130 were annotated as the lipid phosphate phosphatase (LPP) gene, showing differential expression between FS and SS ( Table S5). The Arabidopsis miR171 family with their targets are involved in shoot branching (Wang et al., 2010;Manavella et al., 2013) and gibberellin (GA)-DELLA signaling-mediated chlorophyll biosynthesis (Ma et al., 2014). Moreover, miR171 target protein LiSCL (L. longiflorum Scarecrow-like), also called GRAS, is expressed specifically at the premeiotic phase within anthers (Morohashi et al., 2003).
Recently, studies reported that miR171 and miR167 showed diurnal oscillation expression during the day: up-regulated in the light and down-regulated in darkness (Siré et al., 2009). Another study showed that blue light alters miR167 expression and microRNA-targeted auxin response factor genes in Arabidopsis thaliana suggesting that miR171 and miR167 may be regulated by light or the circadian clock. In this study, tae-miR171a and tae-miR167a show differential expression between FS and SS at different developmental stages (Table 4), as well as their targets (Traes_6DL_26DDCA106.1 for tae-miR171a; Traes_2AL_A7941CB12.2 and Traes_2BL_B2D711406.2 for tae-miR167a) (Table S5). Furthermore, there are consistent negative relationships between expression of tae-miR171a and its target in different developmental stages, suggesting that tae-miR171a is a negative regulator in post-transcriptional gene silencing through base pairing with target mRNAs, which leads to mRNA cleavage or translational repression.
Lots of target genes of miRNAs have been reported in the environmental stress response and are directly involved in some hormone-related signaling pathways that regulate plant development (Shinozaki, 2007). As AP2 transcription factors, AREB/ABFs are involved in the ABA-dependent signal transduction pathway Shinozaki, 1994, 2005;Uno et al., 2000;Kang et al., 2002;Fujita and Yamaguchi-Shinozaki, 2005;Furihata et al., 2006). MYB2 function in ABA-and JA-inducible gene expressions (Abe et al., 2003). The RD26 NAC transcription factor is involved in ABA-and JA-responsive gene expression in stress responses (Abe et al., 2003;Tran et al., 2004;Wang et al., 2016). Recently, it has been reported that up-or downregulation of JAZ (JASMONATE-ZIM DOMAIN) genes might be associated with the abnormal anther dehiscence in TGMS wheat line (Wang et al., 2017). Based on the above analysis, we developed a model to try to describe miRNAs-mediated gene regulation networks and the signaling pathway involved in male sterility of PTGMS wheat line BS366 (Figure 7). As shown in Figure 7, external environment signals, including light and temperature, affect wheat through regulation of miRNAs and their targets. In this model, miR825, miR172, miR156, and miR171 are mainly regulated by light, whereas miR1130/miR1120, miR398, miR159, miR164, and two novel-miRNAs (novel-miR964 and novel-miR2186) may be regulated by light. Moreover, miR167 is regulated by both light and temperature. miR172, miR156, and miR171 interact with their respective target genes (AP2, SPL, and SCL), participating in the GA/abscisic acid (ABA) signaling pathway and GA/auxin signaling pathway to regulate flowering time; miR825, miR167, and miR1120/miR1130 interact with their respective target genes (CaBP, ARF, and LRR), participating in the CRY/PHY signaling pathway, auxin signaling pathway and JA signaling pathway, to modulate pollen development. These interactions result in morphological changes of the wheat pollen. In contrast, the interactions of miR398, miR159, miR164, and novel-miR964/novel-miR2186 with their corresponding target genes (EXPB, MYB, NAM, and PPR) participate in the auxin signaling pathway and the GA/ABA signaling pathway to modulate pollen germination, auxin/IAA responded gene expression and male fertility transition to result in defects of pollen fertility. Photo-thermosensitive genic male sterility in wheat might be the result of the combination of the regulatory interactions described above. It is therefore clear that miRNA-mediated photo-thermosensitive genic male sterility in wheat is a very complicated network regulation system involving the expression of genes, proteins, and regulation of hormone signaling. This study may offer important clues for understanding the different ways gene expression is regulated during male fertility transition by functional description and target analysis of miRNAs in wheat. These miRNA profiles also provide an important theoretical foundation for two-line hybrid wheat molecular breeding and may produce insights for further investigations of PTGMS wheat lines.

CONCLUSIONS
In this study, data from small RNA-seq and degradome-seq were analyzed to examine the regulation of miRNAs of the PTGMS wheat line BS366 during male fertility transition. The expression and correlation of 11 selected flowering and/or male sterility-related miRNAs/targets were detected by qRT-PCR in different fertility conditions during male fertility transition. Most miRNAs/targets had negatively correlated expression and may participate in the signaling pathway during male fertility transition. In addition, a miRNA-mediated regulation mechanism during PTGMS in the wheat line BS366 was developed via more extensive function analysis. A number of newly identified miRNAs were also detected. Our data could be used as a benchmark for future studies of the molecular mechanisms of PTGMS in other crops.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2017. 01370/full#supplementary-material Figure S1 | The first nucleotides analysis of the differentially expressed 21nt miRNAs showed that the bases rates of U, C, A and G.