Identification of Sex Differentiation-Related microRNAs in Spinach Female and Male Flower

Sex determination and differentiation is an important biological process for unisexual flower development. Spinach is a model plant to study the mechanism of sex determination and differentiation of dioecious plant. Till now, little is known about spinach sex determination and differentiation mechanism. MicroRNAs are key factors in flower development. Herein, small RNA sequencing was performed to explore the roles of microRNAs in spinach sex determination and differentiation. As a result, 92 known and 3402 novel microRNAs were identified in 18 spinach female and male flower samples. 74 differentially expressed microRNAs were identified between female and male flowers, including 20 female-biased and 48 male-biased expression microRNAs. Target prediction identified 22 sex-biased microRNA-target pairs, which may be involved in spinach sex determination or differentiation. Among the differentially expressed microRNAs between FNS and M03, 55 microRNAs were found to reside in sex chromosome; one of them, sol-miR2550n, was functionally studied via genetic transformation. Silencing of sol-miR2550n resulted in abnormal anther while overexpression of sol-miR2550n induced early flowering, indicating sol-miR2550n was a male-promoting factor and validating the reliability of our small RNA sequencing data. Conclusively, this work can supply valuable information for exploring spinach sex determination and differentiation and provide a new insight in studying unisexual flower development.


Introduction
Sex is a special character of biology. Most animals are unisexual and the mechanism of sex determination mammalian is already clearer. In plant kingdom, only 6% of angiosperms are dioecious plants and their sex determining mechanism are not well studied [1]. The gender origins of dioecious plant are not uniform, which enhances the difficulty to uncover plant sex determining mechanism. At present, there are two main ideas about the mechanism of sex determination in dioecious plants: one is the two-gene mutation model, and the other is the single-gene mutation model. For two-gene mutation model, the appearance of dioecious plants may be due to the mutation of two genes, M gene, which controls stamen fertility, and F gene, which regulates pistil development [2,3]. For single-gene mutation model, the emergence of dioecious plants may be caused by the two independent mutations of the same gene Q: one is loss-of-function mutation resulting in male sterility, and the other is gain-of-function mutation inducing female sterility [4][5][6]. From these two models, it can be seen that the mutation of flower development-related genes is the most important factor leading to the generation of plant sex. To date, the sex-determining genes or candidates have been identified in a few angiosperms such as Vitis vinifera [7,8], Diospyros lotus [9], Populus trichocarpa [10][11][12], Ficus carica [13], Date palm [14], Fragaria octoploids [15], Actinidia chinensis [16] and Asparagus officinalls [17]. The function of the sex-determining genes was studied only in few plants and the mechanism of sex determination remains largely elusive. Sex differentiation is an important process along with sex determination. Researches on sex differentiation may supply clues to uncover sex determination. Many genes related to sex differentiation have been identified in different plants such as in Silene latifolia [18,19], Carica papaya [20].
MicroRNAs (miRNAs) are a class of endogenous single-stranded non-coding small RNA molecules with a length of about 21-24 nt. The initial transcription product pri-miRNA is spliced into a single strand RNA precursor (pre-miRNA) with hairpin structure, and then formed the mature miRNA with regulatory function. MiRNAs negatively regulate their target messenger RNAs (mRNA) through mRNA degradation, translation inhibition or chromatin modification [21][22][23]. MiRNAs play important regulatory roles in plant morphogenesis, nutrient balance, biological and abiotic stress [24][25][26][27]. More and more studies have shown that miRNAs played important roles in plant flower development.
Flowering is an important process in plant growth cycle and miRNAs are involved in flower formation processes. In Arabidopsis, miR164c regulated two NAC transcription factors CUP-SHAPED COTYLEDON1 (CUC1) and CUC2 to control the boundary foundation of floral organs [28]. MiR172 can regulate the formation of sepal and petal primordium via inhibiting the translation of floral homeotic gene APETALA2 [29,30]. MiR159, miR167 and miR319 can form a regulation module to influence flower development by targeting transcription factors MYB33, TCP4 and Auxin Response Factor6/8 [31]. In tomato, sly-miR160a was involved in flower development via regulating auxin response factor SlARF10A [32]. The flowers of corn, a hermaphrodite plant, are bisexual in the early development stage but become unisexual in the late stage as the development of pistil or stamen is blocked. It has been reported that zma-miR172E was a member of miR172 family and an insertional mutation in the promoter region of its coding gene resulted in a corn recessive mutant tasseleseed4, in which carpel developed in tassel with no stamen [33]. The sex determining gene of a dioecious plant, persimmon (Diospyros lotus) was identified as a miRNA, OGI, targeting a homeodomain transcription factor MeGI to determine plant sex [9,34]. These studies fully demonstrate that miRNAs participate in flower development, even in the formation of unisexual flowers and moreover miRNA can serve as sex determining gene.
Spinacia oleracea (2n = 12), a dioecious plant, has a pair of homomorphic sex chromosomes (XY) and is at the early stage of sex chromosome evolution [35]. Many sex-linked, sex chromosome-specific and sex-determining-gene-linked molecular markers were developed in spinach, such as T11A, V20A, 45 s rDNA, SP_0018, SpoX and so on [36][37][38][39]. These molecular markers are helpful for locating the accurate sex-determining region (SDR) and further identifying sex-determining genes. Qian et al. (2017) [40], Yu et al. (2021) [41] and Ma et al. (2022) [42] found the SDR locating on the sex chromosome using different method. Moreover, 166 sex-differentiation-related genes and 12 Y-specific genes were successively reported in spinach [40,43]. However, the function of these genes was still unclear and it is still controversial that spinach sex determination is regulated by single-gene mutation model or two-gene mutation model [43]. Spinach genome was firstly published in 2017 and constantly updated with the development of high throughput technology, which supply valuable information for functional genome study [41,42,[44][45][46].
To verify whether miRNA can serve as sex determining gene in spinach, miRNAs resource should be firstly harvested. In 2017, spinach miRNAs were analyzed via in silico prediction in vegetative tissues of spinach [47]. Flower development-related genes are crucial for sex determination and differentiation, so we performed small RNA sequencing using 18 spinach flower samples at three early female and male flower development stages to identify some candidate miRNAs for spinach sex determination or differentiation. Virusinduced gene silencing (VIGS) technology is an efficient way to study gene function in plant [48]. In addition, this technology has been successfully applied in spinach [49,50]. Moreover, RNA interference technology and overexpression technology have been used to study miRNA function [51]. Hence, we performed VIGS and overexpression technology to analyze the function of sol-miR2550n, one of the candidate miRNAs, and to further validate the accuracy of our sequencing data.

Small RNA Sequencing
The roles of miRNAs in flower development are not clear in spinach till now. To characterize miRNAs in spinach, we performed small RNA sequencing using 18 spinach flower samples and harvested 13,911,586 high quality reads (mean value of 18 samples data) ( Table 1). After filtering, at least 91.37% of high quality reads were selected as small RNA clean tags, the primary length distribution of which is 24 nt (Figure 1). All small RNA clean tags were aligned in GeneBank database (Release 209.0), Rfam database (11.0), spinach genome (http://www.spinachbase.org/, accessed on 23 July 2021) and miRBase database (Release 21). Finally, 92 known miRNAs and 3402 novel miRNAs were identified in 18 spinach flower samples ( Table 2).   3402 Note: a , the ratio of known miRNA Tags number to Tags total; b , the ratio of novel miRNA Tags number to Tags total; c , the total number of known or novel miRNAs identified in all samples.

Sex-Biased miRNAs of Spinach
According to the expression of miRNA, differentially expressed miRNAs (DE miRNAs) were analyzed among female and male flower samples. There were 431 DE miRNAs between FNS and M03, 451 DE miRNAs between FNB and M05, and 781 DE miRNAs between FYS and M10 ( Figure 2a). Sex-biased genes, which exhibit significantly higher expression in one sex than in the other sex, always act downstream of  Note: a , the ratio of known miRNA Tags number to Tags total; b , the ratio of novel miRNA Tags number to Tags total; c , the total number of known or novel miRNAs identified in all samples.

Sex-Biased miRNAs of Spinach
According to the expression of miRNA, differentially expressed miRNAs (DE miR-NAs) were analyzed among female and male flower samples. There were 431 DE miRNAs between FNS and M03, 451 DE miRNAs between FNB and M05, and 781 DE miRNAs between FYS and M10 (Figure 2a). Sex-biased genes, which exhibit significantly higher expression in one sex than in the other sex, always act downstream of sex-determining gene. Hence, identification of the sex-biased genes is helpful to uncover the sex-determination mechanism. Herein, 74 DE miRNAs (7 were known and 67 were novel) were identified between female and male flowers at three early developmental stages (Figure 2a; Table S1); and of them, 20 miRNAs displayed female-biased expression and 48 miRNAs displayed male-biased expression (Figure 2b). Moreover, there were 9 female-specific expression miR-NAs and 17 male-specific expression miRNAs (Figure 2b). These sex-biased and sex-specific expression miRNAs may be involved in spinach sex determination or differentiation. sex-determining gene. Hence, identification of the sex-biased genes is helpful to uncover the sex-determination mechanism. Herein, 74 DE miRNAs (7 were known and 67 were novel) were identified between female and male flowers at three early developmental stages ( Figure 2a; Table S1); and of them, 20 miRNAs displayed female-biased expression and 48 miRNAs displayed male-biased expression ( Figure 2b). Moreover, there were 9 female-specific expression miRNAs and 17 male-specific expression miRNAs ( Figure 2b). These sex-biased and sex-specific expression miRNAs may be involved in spinach sex determination or differentiation.

Target Genes of miRNAs
Target genes of miRNAs were analyzed according to our previous transcriptome data [52]. Totally, 20,460 genes were predicted as targets of 2855 miRNAs and 49,591 target sites were found in these target genes for miRNAs binding (Table 3). It can be seen that one miRNA can regulate more than one gene and one gene can be controlled by more than one miRNA. MiRNA and its target always display antagonistic expression profile. As previously reported [52], 1946 male-biased genes and 961 female-biased genes were found between female and male flowers at three early developmental stages. Herein, we found ten male-biased genes were targeted by five female-biased miRNAs and six female-biased genes were targeted by ten male-biased miRNAs, composing 22 miRNA-target pairs ( Figure 3, Table S2). The 22 miRNA-target pairs may be involved in spinach sex determination or differentiation.

MiRNAs Residing in Sex Chromosome
In the published spinach genome data, chr4 is the sex chromosome where the sexdetermining genes reside [43]. Blasting with spinach genome [44], we found 307 miRNAs (16 known and 291 novel) residing in Chr1, 275 miRNAs (6 known and 269 novel) residing in Chr2, 415 miRNAs (12 known and 403 novel) residing in Chr3, 496 miRNAs (8 known and 488 novel) residing in Chr4 (sex chromosome), 282 miRNAs (17 known and 265 novel) residing in Chr5, 240 miRNAs (2 known and 238 novel) residing in Chr6, and 1479 miRNAs (31 known and 1448 novel) residing in scaffolds (Figure 4). Among 496 miRNAs residing in sex chromosome, there were 55 DE miRNAs between FNS and M03, 56 DE miRNAs between FNB and M05, and 114 DE miRNAs between FYS and M10. However, no intersection was existed among these three data sets, i.e., no DE miRNA was identified between female and male flowers at three developmental stages. Considering that sex determination and differentiation occurs at the earlier flower development stage in spinach, so we further analyzed DE miRNAs between FNS and M03, the earliest stage of the three flower developmental stages. As a result, 14 differential expression genes were found as targets of ten DE miRNAs between FNS and M03 (

Functional Analysis of Sol-miR172 and Sol-miR2550n
To uncover the potential function of these miRNAs residing in sex chromosome, we performed plant transformation experiments for the first time. sol-miR172 (miR172-y) is a DE miRNA between FNS and M03 and showed higher expression level in M03 stage; its homolog in Arabidopsis has been reported to regulate flower development [29], so we selected sol-miR172 as a positive control. sol-miR2550n (novel-m2550-5p) resided in sex chromosome and displayed higher expression level in M03 stage than in FNS stage (Table  S3), so we selected it as a sex-determining candidate. Virus-induced-gene-silencing and heterologous-overexpression were performed for functional analysis of sol-miR172 and sol-miR2550n in spinach. Silencing of sol-miR172 or sol-miR2550n in spinach resulted in abnormal male flower (anther abortion) (Figure 5a,c,f). Overexpression of sol-miR172 or sol-miR2550n in Arabidopsis induced early flowering (Figure 5b,d,e,g,h). The phenotype of overexpression sol-miR172 in Arabidopsis was consistent with Aukerman and Sakai's report [53]. These results verified the accuracy of the small RNA sequencing and analysis. Moreover, this work will supply abundant valuable information to explore the mechanism of spinach sex determination and differentiation.

Discussion
Sex determination and differentiation is an important question for monoecious and dioecious plant. In some species, the sex determining genes have been identified, such as asparagus, kiwifruit and persimmons. However, sex determining gene is still unclear in spinach; it may be transcription factor or noncoding RNA. Hence, we performed small RNA sequencing to search some clues.
It was the first time to sequencing small RNAs in spinach. Totally, 92 known miR-NAs and 3402 novel miRNAs were identified in 18 spinach flower samples (Table 1), which enriched plant miRNA resource. To analyze the miRNAs related to spinach sex determination or differentiation, 74 differentially expressed miRNAs were screened out among female and male flower at three early developmental stages (Figure 2a), includ- (h) statistics of flowering time of sol-miR172-overexpressed Arabidopsis plants, n ≥ 20; "*" represents p-value < 0.05; "**" represents p-value < 0.01; "***" represents p-value < 0.001; error bar represents standard error.

Discussion
Sex determination and differentiation is an important question for monoecious and dioecious plant. In some species, the sex determining genes have been identified, such as asparagus, kiwifruit and persimmons. However, sex determining gene is still unclear in spinach; it may be transcription factor or noncoding RNA. Hence, we performed small RNA sequencing to search some clues.
It was the first time to sequencing small RNAs in spinach. Totally, 92 known miRNAs and 3402 novel miRNAs were identified in 18 spinach flower samples (Table 1), which enriched plant miRNA resource. To analyze the miRNAs related to spinach sex determination or differentiation, 74 differentially expressed miRNAs were screened out among female and male flower at three early developmental stages (Figure 2a), including 20 femalebiased miRNAs, 9 female-specific miRNAs, 48 male-biased miRNAs and 17 male-specific expression miRNAs (Figure 2b). Moreover, 22 miRNA-target pairs were found through target prediction (Figure 3). These 74 sex-biased or sex-specific miRNAs can be served as candidates of spinach sex differentiation, but not sex determination as they were not reside in sex chromosome. Genes residing in sex chromosome, especially in sex-determining region, are important clues to explore spinach sex determination and differentiation. Hence, we analyzed the miRNAs residing in sex chromosome; and then we found 496 miRNAs in sex chromosome, but none was sex-biased miRNA (Figure 4). Considering that sex determinant plays a role in early spinach flower development stage, so we screened out ten DE miRNAs between FNS and M03, two earliest stages among samples (Table 4). 14 targets with opposite expression trend to these ten DE miRNAs were also identified ( Table 4). The function of these ten DE miRNAs residing in sex chromosome will be studied in our future work.
Herein, one of the DE miRNAs residing in sex chromosome, sol-miR2550n, were firstly studied for its function. Through VIGS and overexpression method, we found that silencing of sol-miR2550n in spinach induced abnormal male flower (anther abortion) and overexpression of sol-miR2550n in Arabidopsis resulted in early flowering ( Figure 5). sol-miR2550n showed higher expression level in male flower (M03) than in female flower (FNS) (Table S3), hence down-regulation of its expression influenced the male flower development. However, up-regulation of its expression didn't affect flower structure but induced early flowering. In spinach, the flowering time of male plant is earlier than female plant [54]. Hence, an interesting hypothesis is that up-regulation of a male factor may promote the male trait, such as early flowering. Hence, sol-miR2550n may be a malepromoting factor. The exact molecular mechanism needs to be explored in future work. Meanwhile, sol-miR172, a homolog of the well-known flower-related factor miR172 [53,55], was also studied in this work as a positive control. In spinach, sol-miR172 also showed higher expression level in M03 (male flower) than in FNS (female flower) but not resided in sex chromosome. Overexpression of sol-miR172 in Arabidopsis resulted in early flowering, in accordance with ath-miR172 ( Figure 5). However, down-regulation of sol-miR172 in spinach resulted in anther abortion of male flower, which is not similar with ath-miR172 (influencing flower whorls pattern) ( Figure 5). Such difference may be due to the architecture difference between unisexual flower and bisexual flower and also indicated a new potential regulation pathway of miR172 in unisexual flower development. The molecular mechanism will be studied in our future work.

Plant Materials and Growth Conditions
All plants used in this study were Spinacia oleracea L. cv DA JIAN YE BO CAI, a dioecious plant. Seeds were obtained from U.S. National Plant Germplasm System (https: //npgsweb.ars-grin.gov/gringlobal/search.aspx?, accessed on 23 July 2021, accession number is PI 527332) and grown in an experimental field at Henan Normal University, Xinxiang, China (113.90 • E, 35.32 • N). Female and male flowers collected at three stages were separately used for Small RNA Sequencing and qRT-PCR. Flower samples were same with our previous published paper [52]. Spinach

Library Construction and Sequencing
After total RNA was extracted by TRIzol (Thermo Fisher Scientific, Waltham, MA, USA), the RNA molecules in a size range of 18-30 nt were enriched by polyacrylamide gel electrophoresis (CWBIO, Taizhou, China). Then the 3 adapters were added and the 36-44 nt RNAs were enriched. The 5 adapters were then ligated to the RNAs as well. The ligation products were reverse transcribed by PCR amplification and the 140-160 bp size PCR products were enriched to generate a cDNA library and sequenced using Illumina HiSeq TM 2500 (Illumina, San Diego, CA, USA) by Gene Denovo Biotechnology Co. (Guangzhou, China).
Reads obtained from the sequencing machine included dirty reads containing adapters or low quality bases which would affect the following assembly and analysis. Thus, to acquire clean tags, raw reads were further filtered using fastp software (https://github. com/OpenGene/fastp, version 0.12.4, 1 March 2022) according to the following rules: (1) Removing low quality reads containing more than one low quality (Q-value ≤ 20) base or containing unknown nucleotides (N); (2) Removing reads without 3 adapters; (3) Removing reads containing 5 adapters; (4) Removing reads containing 3 and 5 adapters but no small RNA fragment between them; (5) Removing reads containing ployA (the content of A base in a reads is higher than 70%) in small RNA fragment; (6) Removing reads shorter than 18 nt (not include adapters). After removing rRNA, scRNA, snoRNA, snRNA, tRNA, mRNA degradation fragments and repeat sequences, the clean tags were then searched against miRBase database (http://www.mirbase.org/, 1 March 2022, Release 21) using bowtie (http://bowtie-bio. sourceforge.net/index.shtml, 1 March 2022, version 1.1.2, parameters: -v 0 -best -strata -a) to identify exist and known miRNAs. Then the unidentified clean tags were mapped to reference genome [44]. According to their genome positions and hairpin structures predicted by software Mireap_v0.2 (https://sourceforge.net/projects/mireap/, 1 March 2022), the novel miRNA candidates were identified. The default parameters of software Mireap_v0.2 were as follows:

MiRNA Expression Profiles
Total miRNAs consist of existing miRNAs, known miRNAs and novel miRNAs, based on their expression in each sample, the miRNA expression level was calculated and normalized to transcripts per million (TPM). The formula is as follows: TPM = Actual miRNA counts/Total counts of clean tags*10 6 .
To identify DE miRNAs across samples or groups the formula was shown as follows: We identified miRNAs with a fold change ≥ 2 and p value < 0.05 in a comparison as significant DE miRNAs.

Target Gene Prediction
Based on the sequences of the exist miRNAs, known miRNAs and novel miRNAs, the candidate target genes were predicted. The software patmatch (ftp://ftp.arabidopsis. org/home/tair/Software/Patmatch/, 1 March 2022, version 1.2) was used to predict target genes. The default parameters were as follows: (1) No more than four mismatches between sRNA & target (G-U bases count as 0.5 mismatches) (2) No more than two adjacent mismatches in the miRNA/target duplex Mir-X miRNA First-Strand Synthesis Kit (TaKaRa, Shiga, Japan) was used to perform reverse transcription of miRNA. Prime-Script™ RT reagent Kit with gDNA Eraser (TaKaRa, Shiga, Japan) was used to synthesis the cDNA of total RNA. Mir-X miRNA qRT-PCR TB Green Kit (TaKaRa, Shiga, Japan) and TB Green ® Premix Ex Taq™ II (Tli RNaseH Plus) (TaKaRa, Shiga, Japan) were used to perform the qRT-PCR of miRNA and mRNA, respectively. qRT-PCR was carried out on LightCycler ® 480 System (Roche, ) according to the manufacturer's instructions. SpoEF, SpoUBQ, At5sRNA, AtACTIN2 were used as an internal control for normalization. The relative expression level of transcript was calculated with 2 −∆∆CT method. All primers sequences are listed in Table S4.

Vector Construction and Plant Transformation
The pre-miRNA was cloned using Trans-Start ® FastPfu Fly DNA Polymerase (TRANS-GEN, Beijing, China) and inserted into the pTRV2 vector (VIGS) or pLP100-35s vector (overexpression) with ClonExpressIIOne Step Cloning Kit (Vazyme, Nanjing, China). These vectors were transformed into Agrobacterium competent cells (GV3101, WEIDI, Shanghai, China). Spinach seedlings with four leaves were injected with Agrobacterium infection buffer including pTRV2-pre-miRNA vector; flower phenotype was observed after about 30 days. Arabidopsis flower buds were dipped in Agrobacterium infection buffer including pLP100-35s-pre-miRNA vector; positive T1 plants were used to observe phenotype.

Conclusions
It was the first time to sequencing small RNAs in spinach. 92 known and 3402 novel miRNAs were identified, which enriched plant miRNA resource. 74 DE miRNAs including male-biased/specific and female-biased/specific expression miRNAs were identified between male and female flowers. Moreover, target prediction identified 22 sex-biased microRNA-target pairs, which may be involved in spinach sex determination or differentiation. Genes residing in sex chromosome, especially in sex-determining region, are very important for sex determination or differentiation. Hence, miRNAs residing in sex chromosome were analyzed and 55 DE microRNAs between FNS and M03 were identified; and one of them, sol-miR2550n, was analyzed via genetic transformation. Silencing of sol-miR2550n resulted in abnormal anther while overexpression of sol-miR2550n induced early flowering, indicating sol-miR2550n may be a male-promoting factor. Conclusively, our work can supply valuable information for exploring spinach sex determination and differentiation and provide a new insight in studying unisexual flower development.

Data Availability Statement:
Sequencing data that support the findings of this study has been deposited in the NCBI SRA database [https://www.ncbi.nlm.nih.gov/sra/PRJNA796451, accessed on 17 March 2022]. All data generated or analyzed during this study are included in this published article.