Differential Expression of MicroRNAs Between 21A Genetic Male Sterile Line and Its Maintainer Line in Cotton (Gossypium hirsutum L.)

MicroRNAs (miRNAs) are a family of 19to 25-nucleotides small non-coding RNAs that play important roles in regulation of many developmental processes in animals and plants, including flower development. In this study, the cotton (Gossypium hirsutum L.) genetic male sterile (GMS) line and its maintainer line were used to detect miRNAs related to male sterility. Using a bioinformatic approach and miRNA microarray analysis, we identified 33 potential miRNAs belonging to 27 families in cotton. miRNA microarray and quantitative real-time PCR assays showed that miRNAs expressed differently at the sporogenous cell, pollen mother cell and pollen grain stages between the cotton GMS and its maintainer line. These cotton miRNAs may regulate 35 potential target genes involved in cotton growth and development, signal transduction and metabolism pathways. The expressions of four targets were contrary to the expression of their corresponding miRNAs. These findings enhance our understanding of the roles of miRNAs during fertile and sterile anther development.


Introduction
Male sterility is not only a ubiquitous biological phenomenon, but also an important way to utilize heterosis in plants. The process of pollen development mainly includes the formation of the microspore and male gametogenesis, and is similar throughout higher plant species (Wilson & Zhang, 2009). Through meiosis, each pollen mother cell can form four microspores in a tetrad. After the microspores lose the thick wall of the tetrad, the free uninucleate microspore undergoes two mitotic divisions to form the male gametophyte also known as a three-celled organism (McCormick, 1993(McCormick, , 2004. Anther development involves a series of complex physiological and biochemical process and, therefore, abnormality at any point in this process can cause male sterility (McCormick, 2004).
MicroRNAs (miRNAs) are a newly discovered class of 19-to 25-nucleotides (nt) small non-coding RNAs that regulate gene expression by mRNA cleavage or translational repression in eukaryotes (Carrington & Ambros, 2003;Schwab et al., 2005;Axtell et al., 2007). They are processed from their double-stranded precursors by RNase III enzyme Dicer (Cobb et al., 2005). Approaches including cloning, bioinformatics and deep sequencing have been used in various organisms to identify miRNAs. Currently, there are 21,264 miRNA entries in the miRBase (Release 19, August 2012; http://microrna.sanger.ac.uk/), which includes 39 miRNAs identified and confirmed in cotton (Gossypium hirsutum L.). In plants, miRNAs have been shown to play important roles in organismal development, hormonal response, signal transduction and stress response (Pasquinelli et al., 2005;Sunkar, 2010). Many miRNAs have been identified and analyzed in the pollen development of Arabidopsis, rice, maize and hickory in recent years. Chambers and Shuai (2009) identified 37 known miRNAs in Arabidopsis pollen by using miRCURY LNA array technology. Similarly, Grant-Downton et al., (2009) not only detected 33 different miRNA families in mature pollen independently by using quantitative RT-PCR (qRT-PCR) but also identified seven putative novel miRNAs. Peng et al. (2012) identified 47 known and 57 novel miRNAs during male gametophyte development in rice from sRNA-Seq datasets, and also found that most of the new miRNAs were pollen-specific and not conserved among species. Seven novel miRNA families and one known miRNA were cloned by constructing a small RNA library from a mixture of anthers from a cytoplasmic male sterile line and its maintainer in maize, all of which showed different expression patterns according to qRT-PCR analysis between the two lines (Shen et al., 2011). In hickory, 39 conserved, two novel and ten potential novel miRNAs were identified by using Solexa technology to sequence two small RNA libraries from two floral differentiation stages . In addition, there is increasing evidence showing that the function of miRNAs such as miR156, miR159, miR164, miR167, miR172 and miR319 is significant during flower development, from the earliest to very late stages Palatnik et al., 2007;Glazinska et al., 2009;Nag et al., 2009). Through researching the mechanism of miRNA action, we can better understand the relationship between miRNA and plant male sterility.
Genetic male sterility (GMS), controlled by nuclear genes, is useful for hybrid seed production in crops. In the past few years, the mechanism of GMS has been intensively studied by using male sterility mutants in cotton. Most research has involved the differences in carbohydrate, free amino acids, enzymes, endogenous hormones and gene expression between the process of fertile and sterile anther development and molecular mapping of genic male-sterile genes (Song et al., 2004;Ma et al., 2007;Chen et al., 2009). With the development of bioinformatics and experimental methods, miRNAs were increasingly identified in cotton (Qiu et al., 2007;Zhang et al., 2007;Li et al., 2012;Yin et al., 2012;Wei et al., 2013). However, to our knowledge, the potential miRNAs involved in cotton anther development remain few.
A cotton 21A male sterility line was found to be under the control of a recessive gene. There was no pollen or a little in the mature anthers of male sterile plants, and no seed set by self-crossing. The microsporocyte development of 21A plant was abnormal at meiosis, and anther abortion occurred mainly at the pollen grain stage. Allelic tests indicated that the male sterility observed in 21A was non-allelic to all male sterile genes identified previously in upland cotton (Feng et al., 2010). In this study, 21A and its maintainer line were used to detect miRNAs and their target genes from cotton anthers using computational approaches. Expression patterns of potential miRNAs and their targets were further analyzed by miRNA microarray assay and qRT-PCR techniques.

Plant Materials
The new upland cotton (G. hirsutum L.) 21A GMS line (SL) and its maintainer line (ML) were grown respectively at different locations of Shandong Agricultural University, Taian, Shandong Province, China. At full bloom stage, developing anthers at the sporogenous cell, pollen mother cell and pollen grain stages (bud lengths < 3.5, 3.5-5.0 and > 5.0 mm, respectively) were separately harvested. The anther samples were quickly frozen in liquid nitrogen and stored at -80 °C.

RNA preparation
Total RNA was isolated from each sample using the TRIzol® reagent (Invitrogen). The quality of the RNA was assayed using the 2100 Expert_ Eukaryotic Total RNA Pico assay.

miRNA and Target Gene Prediction
Computational approaches were used to predict the potential miRNAs and their target genes, using methods described by Yin et al. (2008). The major steps for searching potential miRNAs and target genes in cotton are shown in Supplementary Figure 1. The cotton expressed sequenced tag (EST), genomic survey sequences (GSS) and mRNA databases were obtained from the NCBI databases. To identify new potential miRNAs, all previously known plant miRNA sequences obtained from the miRNA Registry Database (Release 18, http:/www.mirbase.org) were used in BLAST searches against the cotton EST and GSS databases. All BLAST results were saved. The sequences with ≤ 3 mismatches compared with the query miRNA sequences were chosen manually for further options. The hairpin structures of the selected sequences were predicted using Mfold 3.2. Redundant nucleotides were deleted in order to obtain the mature sequences of predicted miRNAs. To search potential target genes in cotton, the mature miRNA sequences identified were used for a BLAST search in the EST database and genome sequence database of G. Raimondii . Sequences with ≤ 4 mismatches, compared with the query miRNA sequences, were chosen. Because the Unigene database contains many transcript sequences that appear to come from the same transcription locus (Yin et al., 2012), we removed redundant sequences.

miRNA Microarray Assay
Microarray assay was performed using a service provider (LC Sciences). To determine differential expression of www.ccsenet.org/jps Journal of Plant Studies Vol. 3, No. 1; miRNAs between 21A GMS line and its maintainer line samples, 1014 previously known and 33 predicted miRNA probes in triplicate were contained in each miRNA microarray. Three chips were separately used for groups of the sporogenous cell, pollen mother cell and pollen grain stages. The assay started from 2-5 µg of total RNA, which was size fractionated using a YM-100 Microcon centrifugal filter (Millipore) and the small RNAs (< 300 nt) isolated were 3'-extended with a poly(A) tail using poly(A) polymerase. An oligonucleotide tag was then ligated to the poly(A) tail for later fluorescent dye staining. Hybridization was performed overnight on a µParaflo™ microfluidic chip using a micro-circulation pump (Atactic Technologies). The hybridization melting temperatures were balanced by chemical modifications of the detection probes. Hybridization used 100 µL of 6 × SSPE buffer (0.90 M NaCl, 60 mM Na 2 HPO 4 and 6 mM EDTA at pH 6.8) containing 25% formamide at 34 °C. After hybridization, detection was by fluorescence labeling using tag-specific Cy3 and Cy5 dyes. Hybridization images were collected using a laser scanner (GenePix 4000B, Molecular Device) and digitized using Array-Pro image analysis software (Media Cybernetics). Data were analyzed by first subtracting the background and then normalizing the signals using a LOWESS filter. For two color experiments, the ratio of the two sets of detected signals (log 2 transformed, balanced) and p-values of the t-test were calculated; differentially detected signals were those with p < 0.01.

qRT-PCR Analysis
qRT-PCR was used to analyze the relative expression of miRNAs and target genes. The specific miRNAs and target genes primers used are given in Supplementary Table 1. For qRT-PCR of miRNAs, however, Uni-miR qPCR Primer (TAKARA, Dalian, China) was used instead of the PCR forward primer. U6 small nuclear RNA served as an internal standard in miRNA amplification experiments, and Ubiquitin gene (UBI) in experiments of gene qRT-PCR. For the first-strand cDNA synthesis experiment, One Step PrimeScript® miRNA cDNA Synthesis Kit (TAKARA) was used following the manufacturer's instruction. From each sample, 4 μg of total RNA was converted to cDNA in a 20-μL reaction system which contained 10 μL of 2 × miRNA reaction buffer mix, 2 μL of 0.1% BSA and 2 μL of miRNA PrimeScript® RT Enzyme Mix. qRT-PCR for miRNAs and target genes was performed using SYBR ® Premix Ex Taq TM II (TAKARA) and carried out in Bio-RAD iCycler iQ5 Machine with following cycling profile: 95 °C at 30 s, followed by 40 cycles of 95 °C at 15 s and 62 °C at 30 s. The reaction mixture was prepared in triplicate. Each 25-μL reaction mixture contained 12.5 μL of SYBR®Premix Ex TaqTM II (2×), 1 μL of PCR forward primer (10 μM), 1 μL of PCR reverse primer (10 μM) and 2 μL of fivefold diluted cDNA template. After reactions were completed, the 2 -ΔΔCt method [ΔCt = Ct (miRNA or target gene) -Ct (U6 or UBI), ΔΔCt = ΔCt (SL) -ΔCt (ML), 2 -ΔΔCt = Relative Expression] was used to calculate relative expression of miRNAs or target genes (Livak et al., 2001).

Identification of potential cotton miRNAs
In our study, a homology search was used to identify conserved miRNAs in cotton. Following a set of strict filtering criteria, a total of 33 cotton miRNAs belonging to 27 miRNA families were identified in cotton EST and GSS databases (Table 1). The sequences of these newly identified miRNAs were 20-22 nt in length. A total of 24 miRNA mature sequences began with 5'-uridine a characteristic of plant miRNAs. In addition, the length of miRNA precursors was 48-513 nt, but a majority of them (72.7%) required 65-151 nt, which is similar to that observed in many other species (Nasaruddin et al., 2007;Yin et al., 2008). Furthermore, A + U contents of those pre-miRNAs sequences were calculated in the range of 36.76-74.77%. All of those miRNA precursors can be folded to the typical hairpin structure of the miRNA family (Supplementary Figure 2). The minimal folding free energy index (MFEI) is an available standard for distinguishing miRNAs from other types of coding or non-coding RNAs. In our results, the average MFEIs of the 33 pre-miRNAs was 0.81; a majority of the identified cotton pre-miRNAs had MFEI values in the range 0.71-1.29, which is remarkably higher than that for tRNAs (0.64), rRNAs (0.59) and mRNAs (0.62-0.66) (Zhang et al., 2006). For the 33 identified potential miRNAs, miR156, miR395, miR414 and miR426 had 2-3 family members each, while the remainder had only one member. The miR426 family and the other miRNAs, including miR403, miR418, miR477, miR529, miR834, miR846, miR857, miR1044 and miR1063 were identified in cotton for the first time. MFEIs, minimal folding free energy indexes; LP: length of miRNA precursor.

Microarray Analysis
Microarray technique was used to profile global miRNA expression between the 21A GMS line and its maintainer line during anther development. To confirm the predicted cotton miRNAs above, 33 probes complementary to the potential miRNAs were designed and assigned to each chip. Microarray data indicated that 504, 432 and 439 transcripts were detected perfectly at the sporogenous cell, pollen mother cell and pollen grain stages of anther development, respectively. Among them, all 33 potential cotton miRNAs were detected in both groups; 10 were differentially expressed at the three stages of anther development using miRNA microarray (Figure 1).
Analysis of microarray data showed that most miRNAs of miR156, miR164 and miR166 families were up-regulated in the SL compared with the ML group at the three stages of anther development. Compared with the ML group, miR159, miR172 and miR319 families were up-regulated at the sporogenous cell stage, but down-regulated at the pollen mother cell and pollen grain stages in the SL group ( Figure 1). These data suggested that anther development involved a series of expressions of conserved miRNAs. The regulation of these miRNAs from other species implied that they might play significant roles during anther development in cotton. Each row represents a miRNA and each column represents a stage of anther development. ARRY0X sporogenous cell, ARRY1X pollen mother cell and ARRY2X pollen grain stages. The clustering is performed on log 2 (ML/SL) ratio on variation across samples. The intensities of the color represent the relative magnitude of fold changes in log values: -3.0 indicates that the miRNA is highly suppressed; while +3.0 indicates that it is highly induced.

miRNA Expression Profiling With qRT-PCR
To validate the microarray data, 12 miRNAs were selected at random for qRT-PCR analysis for the relative expression of miRNAs at the three stages of the cotton anther development respectively. qRT-PCR analysis showed that the expression of 10 miRNAs fitted well with the pattern of microarray analysis (Figure 2). However, there was no obvious differential expression between the SL and ML samples for ghr-miR156b at the sporogenous cell stage. To further determine the relative expression of identified new cotton miRNAs at the different stages of anther development, qRT-PCR was also used for analysis -the expression patterns of these miRNAs were intriguing. Both in the SL and ML samples: (1) the expression levels of ghr-miR156a, ghr-miR156b, ghr-miR171, ghr-miR396 and ghr-miR398 had no significant differences between the sporogenous cell and pollen mother cell stages, but were up-regulated remarkably at the pollen grain stage; (2) the ghr-miR164, ghr-miR166 and ghr-miR482 showed decreased expression at the pollen mother cell compared to the other two stages; and (3) there was no dramatic differential expression for both ghr-miR162 and ghr-miR834 in these three stages. In addition, qRT-PCR analysis also showed no obvious differential expression for some miRNAs between the SL and ML samples, but a high expression for ghr-miR156a, ghr-miR156b, ghr-miR164, ghr-miR166, ghr-miR171 and ghr-miR482 at the pollen grain stage, and a low expression for ghr-miR398 at these three stages in the SL sample. The experimental data show that these miRNAs were differentially regulated at different stages of the cotton anthers between the GMS and the maintainer lines ( Figure 3). Therefore, these differentially expressed miRNAs might play important roles in the process of anther development.

Target Genes of miRNAs and Their Expression Analysis of qRT-PCR Assay
The newly identified miRNAs were used as query sequences for BLASTn searches to identify miRNA targets in the cotton mRNA database. In this study, we identified a total of 35 potential target genes for 20 of the new miRNAs from the cotton mRNA database (Supplementary Table 2). To examine the functional similarity of predicted targets, the cotton entries were selected from the Unigene database. Annotated mRNA targets had various biological functions in plant anther development, with the majority related to growth and development, transcriptional factors and metabolic pathways. Up to now, there have also been several predicted targets encoding unknown proteins. In this study, eight identified miRNA families in cotton targeted a variety of transcription factors, which play important roles in cotton anther development. miR156 targets the squamosa-promoter binding protein transcription factors and squamosa-promoter binding protein-like (SPL) genes. miR164 targets NAC transcription factor, and miR167 targets auxin response factors that play roles in plant signaling transduction and root development. It is interesting to note that those transcription factors are similar to those reported previously as miRNA targets for other plant species (Mallory et al., 2004;Arazi et al., 2005;Wu et al., 2006;Yang et al., 2006;Fu et al., 2012). In this study, we also found that scarecrow-like (SCL) protein 22 and mitochondrial transcription termination factor have near-perfect complementary sites with miR171 and miR426, respectively.
To confirm the prediction of the target genes, we performed qRT-PCR analysis of several target genes between the fertile and sterile lines. The mRNAs of four targets tested showed expression patterns of the potential target genes contrary to the expression of their corresponding miRNAs ( Figure 3). Therefore, our predictions of the partial target genes were preliminarily verified by qRT-PCR. Further research is needed to verify these target genes using 5' rapid amplification of cDNA ends.

Discussion
miRNAs, such as miR156, miR159, miR164, miR166, miR172 and miR319 play important roles during the development of flower organs and some participate in the process of male sterility. In this study, using bioinformatic and miRNA microarray approaches, 33 miRNAs were identified in cotton anther tissue. The miR426 family and the other miRNAs, including miR403, miR418, miR477, miR529, miR834, miR846, miR857, miR1044 and miR1063 were firstly identified in cotton. qRT-PCR showed that 10 miRNAs and their targets had a specific expression pattern differences between the GMS and maintainer lines during the three stages of cotton anther development.
During the process of fertile and sterile anther development, two miR156s, miR171, miR396 and miR398 were up-regulated remarkably at the pollen grain stage, but down-regulated at the sporogenous cell and pollen mother cell stages. Moreover, miR156s and miR171 were up-regulated at the pollen grain stage of sterile anther development, but miR396 down-regulated. Specifically for miR398, expression decreased at the three stages of sterile anther development. In some plant species, miR156s, miR171 and miR396 have been shown to target SPL transcription factors, SCL transcription factors and growth-regulating factors (GRF), respectively (Llave et al., 2002;Wang et al., 2009;Wang et al., 2011). The functions of these targets were found to be related to flower development, e.g. controlling flowering time (Llave et al., 2002;Xie et al., 2006). GRFs, a plant specific family of transcription factors known to control cell proliferation in leaves (Rodriguez et al., 2010), especially also play important regulatory roles in floral organ development. Reduced expression of NtGRF-like genes was reported to result in abnormalities such as many more petals, stamens and carpels, but lower fertility than in wild-type tobacco flowers (Yang et al., 2009). miR398 has been shown to target COX5b.1, two copper superoxide dismutase (CSD) enzymes and copper chaperone for superoxide dismutase (CCS1) in Arabidopsis (Jones-Rhoades and Bartel, 2004;Sunkar and Zhu, 2004;Bouché, 2010). Those targets were involved in different abiotic stresses. However, miR398 has not been previously reported to be related to male sterility. Significantly, we observed that expression of miR398 was down-regulated during the process of sterile anther development, suggesting that miR398 may be strongly related to male sterility.
miR164, miR166 and miR482 showed decreased expression from the pollen mother cell to the sporogenous cell stage, but increased from the pollen grain to the pollen mother cell stage in the GMS and maintainer lines; significantly, there was high expression for these miRNAs in sterile tissue at the pollen grain stage. Recent studies demonstrated that miR164, miR166 and miR482 target NAC transcription factor family genes, Class III HD-ZIP transcription factor family genes and nucleotide-binding site leucine-rich repeat (NB-LRR) genes, respectively (Emery et al., 2003;Mallory et al., 2004;Ko et al., 2006). Analysis of early extra petals1 Arabidopsis mutants confirmed that miR164c prevents extra petals in early-arising flowers by repressing CUC1 and CUC2, which were members of NAC transcription factor family (Baker et al., 2005). In addition, it was shown that miR166 regulated vascular development in Arabidopsis inflorescence stems through directed cleavage of ATHB15 mRNA (Kim et al., 2005), one of the predicted targets of miR166 in the present study. Through the phenotypes of a series of double mutants, the miR166/165 signaling pathways function in parallel to the WUS-CLV pathway in regulating the shoot apical meristem and floral development (Jung & Park, 2007). Eitas and Dangl (2010) found that NB-LRR proteins functioned in microbial recognition as activators of defense responses in plants. However, the function of protein encoded by targets of miR482 remained unknown in our prediction results. Both in the GMS and its maintainer line, there was no dramatic differential expression for miR162 and miR834 during anther development. These data suggest that miR162 and miR834 were not related to fertile and sterile anther development.
In conclusion, we identified 33 miRNAs belonging to 27 families from cotton EST and GSS databases based on bioinformatics and microarray analysis. Ten miRNA families were observed for the first time in cotton. These new miRNAs are highly conserved among plant species and some showed diverse expression patterns during the cotton anther development between the GMS and its maintainer line. We also obtained 35 potential miRNA targets which included transcription factors regulating plant growth and development, enzymes and proteins involved in cotton metabolism and signal transduction, implying that miRNA may play important roles in cotton anther development. Further research on function analysis of the target genes would give a better understanding of the anther abortion mechanisms of the GMS line in cotton.