{'fraction_non_alphanumeric': 0.046044098573281456, 'fraction_numerical': 0.04423753937372615, 'mean_word_length': 4.91965999451604, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Heat‐responsive microRNAs and phased small interfering RNAs in reproductive development of flax

Abstract Plants will face increased heat stress due to rising global temperatures. Heat stress affects plant reproductive development and decreases productivity; however, the underlying molecular mechanisms of these processes are poorly characterized. Plant small RNAs (sRNAs) have important regulatory roles in plant reproductive development following abiotic stress responses. We generated sRNA transcriptomes of reproductive bud stages at three different time points to identify sRNA‐mediated pathways responsive to heat stress in flax ( Linum usitatissimum ). With added sRNA transcriptomes of vegetative tissues, we comprehensively annotated miRNA and phasiRNA‐encoding genes (PHAS) in flax. We identified 173 miRNA genes, of which 42 are newly annotated. Our analysis revealed that 141 miRNA genes were differentially accumulated between tissue types, while 18 miRNA genes were differentially accumulated in reproductive tissues following heat stress, including members of miR482/2118 and miR2275 families, known triggers of reproductive phasiRNAs. Furthermore, we identified 68 21‐PHAS flax loci from protein‐coding and noncoding regions, four 24‐PHAS loci triggered by miR2275, and 658 24‐PHAS‐like loci with unknown triggers, derived mostly from noncoding regions. The reproductive phasiRNAs are mostly downregulated in response to heat stress. Overall, we found that several previously unreported miRNAs and phasiRNAs are responsive to heat stress in flax reproductive tissues.

through impacts on pollen development, pollination, and embryo formation (Hedhly et al., 2009). Although there are studies that assess the impacts of higher temperature on growth, physiological processes, and yields on crops such as wheat, maize, rice, canola, and flax (Cross et al., 2003;Pokharel et al., 2020;Wang et al., 2016;Zhao et al., 2017), limited knowledge is available on the underlying molecular mechanisms that play roles in heat stress; such knowledge would enable the development of thermo-tolerant crops.
Heat shock proteins (HSPs) are major players in heat tolerance (Hasanuzzaman et al., 2013). HSPs are produced in response to heat stress after activation by heat stress transcription factors (HSFs) such as MULTIPROTEIN BRIDGING FACTOR 1C, WRKYGQK motifcontaining proteins, and basic leucine zipper proteins (He et al., 2019;Zhang et al., 2017). Multiple thermosensors initiate HSF activation after sensing heat stress. The conserved heat tolerance process in plants involves regulatory cascades of heat-responsive genes as well as other regulators such as noncoding RNAs that play a role in modulation of gene expression, and these noncoding molecules are yet to be identified.
Plant sRNAs are implicated in modulating expression of genes related to heat stress (He et al., 2019). These sRNAs are regulatory elements that modulate gene expression, genome-wide, by transcriptional and post-transcriptional gene silencing. Typically, sRNAs are 20-24 nucleotides (nt) in length and are diced by DICER-LIKE proteins (DCLs) from hairpin or double stranded RNA precursors formed by RNA-DEPENDENT RNA POLYMERASEs (RDRs). The sRNAs form a silencing complex with Argonaute (AGO) proteins to modulate the expression of genes and transposable elements (TEs). The major classes of plant sRNAs are as follows: microRNAs (miRNAs), phased secondary small interfering RNAs (phasiRNAs), and secondary small interfering RNA (siRNAs) (Axtell, 2013;Axtell et al., 2006;Fei et al., 2013;Reinhart et al., 2002;Zhai et al., 2011).
The miRNAs are the key regulators in the attenuation of gene expression, by transcript cleavage and degradation or translation inhibition in a homology-dependent manner (Zhai et al., 2011). There have been reports of heat-responsive miRNAs such as miR398, miR156, miR159, miR160, miR167, miR396, and miR408 during vegetative and reproductive stages of plants like Arabidopsis, sunflower, maize, and tomato (Giacomelli et al., 2012;Guan et al., 2013;He et al., 2019;Keller et al., 2020;Lin et al., 2018;Stief et al., 2014). PhasiRNAs are generated from mRNA transcripts of coding and noncoding regions of the genome mainly after the cleavage of trigger miRNAs (Fei et al., 2013;Johnson et al., 2009;Zhai et al., 2011). A special class of phasiRNAs that are produced from trans-acting siRNA loci (TAS) is called tasiRNAs (Fei et al., 2013;Howell et al., 2007). In Arabidopsis, tasiRNAs from the TAS1 locus are found to target transcripts of HEAT-INDUCED TAS1 TARGET 1 (HTT1) and HTT2 proteins that accumulate in response to heat stress (Li et al., 2014). Complete loss of reproductive 24-nt phasiRNAs after knock out of DCL5, a reproductive phasiRNA pathway gene, in maize, produces male sterility at high temperatures (Teng et al., 2020), while mutation in reproductive 21-nt phasiRNAs-generating loci produces photoperiod-and thermosensitive male sterility in rice (Fan et al., 2016;Zhou et al., 2012).
Flax is a widely cultivated crop species from the family Linaceae known for oil, fiber, and medicinally important compounds. Heat stress in flax reduces boll formation and seed set, and it may be due to the failure of pollen tubes to reach ovules (Cross et al., 2003). In this study, flax plants were subjected to 7-and 14-day heat treatments after flower initiation in a growth chamber to simulate heat stress conditions in the field. We harvested three reproductive bud stage tissues at three time points and generated sRNA transcriptomes to identify sRNA-mediated pathways responsive to heat stress. With added publicly available sRNA transcriptomes of leaves, stem, and root, we comprehensively annotated miRNA-and phasiRNA-spawning genes (PHAS) in flax. We found that several miRNAs and both types of reproductive phasiRNAs (21 and 24 nt) are mostly downregulated in response to heat stress. Overall, we found that several previously unreported flax miRNAs and phasiRNAs are responsive to heat stress in reproductive tissues.

| Comprehensive annotation of miRNA in flax
In the latest release of miRBase version 22 (http://www.mirbase.org), there are 124 miRNA precursors deposited for flax, mainly identified from computational approaches (Barvkar et al., 2013). There are studies that annotated flax miRNAs based on sRNA sequencing in control, nutrient, saline, and alkaline stress conditions (Melnikova et al., 2016;Yu et al., 2016), but their miRNA lists are not available from miRBase; these studies predicted a high number of novel candidate miRNAs , in addition to observing known miRNAs, suggesting that most of them may be false positives. Because the highly studied model plant Arabidopsis only has a few hundred miRNAs, caution may be warranted for reports of more than 100 novel miRNAs in a species (Axtell & Meyers, 2018). Further, all these flax miRNA identifications were carried out using an older scaffold-based genome (Wang et al., 2012). Recently, a chromosome-scale genome assembly has been released for flax (You et al., 2018).
To comprehensively annotate miRNA precursors in the latest flax genome assembly version 2 (ASM22429v2) (You et al., 2018), sRNA libraries were sequenced from leaves, stem, and root tissues. In addition, to identify heat-responsive miRNAs in flax, sRNA transcriptomes from reproductive bud stages (Figure S1a) at three time points were sequenced (see Section 4 for details). In total, we generated 24 sRNA libraries, and clean reads (after removing rRNA/tRNA sequences) were used for miRNA prediction (Table S1). To accurately identify miRNA precursors, we used ShortStack and miR-PREFeR, some of the bestperforming miRNA annotation tools in terms of precision and accuracy (Johnson et al., 2016;Lei & Sun, 2014). We identified 173 candidate miRNA precursors yielding 106 nonredundant mature miRNAs ( Figure 1, outer circle; Table S2). Among them, 131 precursors produced known miRNAs, while 42 precursors generated novel miRNAs which contained at least five variants (mismatches plus unaligned nucleotides) with Viridiplantae miRNAs derived from miRBase version 22. We identified miRNAs from 31 known miRNA families; among them, the families having the greatest numbers of precursors encoded into the genome were miR167 (14), miR482/2118 (12), miR156 (11), miR166 (9), and miR160 (9) (Figure S1b, Table S2).
2.2 | Differentially accumulating miRNA precursors in flax tissues and heat treatments miRNA species (miRNA-miRNA* pairs) are produced from precise cleavage of stem-loop RNA precursors, called miRNA precursors. The miRNA species accumulate differentially to perform their roles in developmental processes in vegetative and reproductive tissues. To assess miRNA precursor abundance, we summed the abundance of sRNAs that mapped to the entire region of the hairpin, and this was used also for subsequent differential expression (DE) analyses. We used this approach for DE of miRNA loci because it accounts for the abundance of the miRNA mature sequence, its variants, and the star sequence, instead of using the abundance value of the single mature miRNA. Here, we characterized the abundance of miRNA loci in different tissue stages and summarized the results (Table 1). In total, we found 140 miRNA loci varying in abundance between vegetative and reproductive tissues; among them, 89 were enriched in vegetative tissues ( Figure 2a, Table S3), while 51 were enriched in reproductive tissues ( Figure 2b, Table S3). Vegetative tissue-enriched miRNA loci F I G U R E 1 Genome-wide distribution of MIR genes and PHAS loci in flax. Chromosomal position of 173 MIR genes (outer circle, pink dots), 68 21-PHAS loci (middle circle, green dots), and 662 24-PHAS/24-PHAS-like loci (inner cycle, red dots). The MIR genes and 21-PHAS loci were annotated in each chromosome. The outermost cycle represents chromosomes from 1 to 15 as numbered were divided into four sets based on their enrichment into different vegetative tissues. Set I represents miRNA loci more abundant in both leaves and root and includes members of miRNA families miR160, miR156, miR166, miR393, miR395, miR845, miR482, and several novel miRNAs (22, 26, and 46). Set II represents miRNA loci enriched mostly in root tissues only and composed of members of miR164, miR156, miR395, miR159, miR399, miR394, and several novel miRNAs (20, 38, and 13). Set III, mostly enriched in both stem and root tissues, is composed of members of miR8175, miR396, miR395, miR169, miR172, miR171, miR159, and novel miRNAs (15, 16, 25,17, 10, 14, 11, 36, and 3). Set IV-related loci were enriched mostly in leaves only and composed of members related to miR167, miR171, miR166, miR162, miR172, miR396, miR530, miR396, miR6445, miR168, miR397, miR408, miR398, and several novel miRNAs (9, 31, 39, and 12). We divided the group of miRNA loci that were enriched more in reproductive tissues into three sets (Set V, VI, and VII). Set V, mostly enriched in all three reproductive tissues along with stem and root, includes members of miR156, miR319, and miR482/2118. Set VI, mostly enriched in reproductive tissues, includes members of miR171, miR390, miR1446, miR482/2118, miR2275, miR169, and novel miRNAs (30, 24, and 7). The last, Set VII, includes miRNA loci that were enriched in reproductive tissues and leaves, and consists of members of miR160, miR393, miR164, miR403, miR171, miR167, and novel miRNAs (5, 18, 33, 23, and 43).
In this study, we subjected flax to two heat treatments during flowering as reported in (Cross et al., 2003) (more in Section 4 for details); and these heat treatments resulted in differential abundances of a total of 18 miRNA loci (Figure 2c, Figure S2a, Table S4). These were mostly downregulated due to heat stress, with one exception: miR8175 ( Figure S2b). The downregulated miRNA loci under heat stress are members of the families of miR482/2118, miR390, miR168, and miR393a, miR399d, miR397b, miR408, and novel_7. Compared with other bud stages, bud stage 2 is the most sensitive to heat stress, and 11 miRNA loci were differentially accumulated after 7 days of heat treatment. There were no differentially accumulated miRNA loci identified between the two different treatments in all stages, but with an increased number of heat stress days, from 7 to 14, more miRNA loci were differentially expressed (Table 2). Overall, members of the miR482/2118 family were enriched in reproductive tissues and they initiate the production of 21-nt reproductive phasiRNAs (Pokhrel, Huang, Bélanger, et al., 2021).

| Genome-wide identification of phasiRNAs in flax
Plant genomes encode loci that produce secondary siRNAs that are generated from a miRNA cleavage site within a polyadenylated RNA; these are generated in a phased pattern and are called PHAS loci.
These loci can be derived from both protein-coding (mRNA) and noncoding (long, noncoding RNA or lncRNA) transcripts. We used small RNA libraries from both vegetative and reproductive tissues to identify PHAS loci in flax. We identified 68 21-PHAS loci (i.e., locigenerating 21-nt phasiRNAs) from both coding and noncoding regions ( Figure 1, middle circle; Table 3; Table S5). Out of 68 loci, 51 were enriched in reproductive tissues (Figure 3a), while 17 loci were enriched in vegetative tissues ( Figure S3a). We identified three loci of miR390-TAS3, the most conserved phasiRNA pathway in land plants (Xia et al., 2015) (example locus in Figure 3b). We found two noncoding TAS-like (TASL) loci enriched in leaves and reproductive tissues triggered by miR167h (example locus in Figure 3c, miR167 variants in Figure 3d) and another 12 noncoding loci enriched in reproductive tissues with unknown triggers (Figure 3a). We found triggers of at least 30 21-PHAS loci, including miR167 as a trigger for a locus encoding a metal tolerance protein, miR168 for AGO1, miR393 for AFB2, miR6445 for NAC, miR408 for a laccase-encoding gene, and miR482/2118 for noncoding reproductive 21-PHAS loci that were recently described (Pokhrel, Huang, Bélanger, et al., 2021). The other gene families for which triggers were not identified included genes encoding BHLH, SCD, RING-type, Zf-C3H4, and BZIP domaincontaining proteins, a eukaryotic translation initiation factor, and an E3 ubiquitin protein ligase (Table 3).
at specific stages of reproductive development, misidentified as 24-nt phasiRNAs as previously described (Polydore et al., 2018). However, a subset of these loci may be derived from a 12-nt phasing mechanism similar to the 24-nt phasiRNAs of Solanaceous species (Xia et al., 2019).

| Differentially expressed phasiRNAs
Because the two classes of phasiRNAs enriched in anthers of angiosperms demonstrate short periods of peak abundance (Pokhrel, (Table 4).
In the case of 24-PHAS loci, for bud stages 1 and 2, 14 days of heat treatment resulted in more downregulated 24-PHAS loci compared with 7 days of heat treatment, although the abundance of 24-PHAS loci in bud stage 3 did not change due to either of the heat treatments ( Table 5). Each of the four reproductive 24-PHAS loci triggered by miR2275 was downregulated in heat-treated samples, while 450 24-PHAS-like loci were also downregulated by heat, out of 658 total loci (Figure 4b, Table S8).

| DISCUSSION
We Furthermore, we identified 68 21-PHAS loci that include three TAS3 loci, two TASL loci, and additional loci from coding and noncoding regions. We found no PHAS loci related to PPR genes or NB-LRR defense genes, which is surprising as these are conserved and numerous PHAS loci in many eudicots (Zhai et al., 2011). There are cases of the near absence or absence of PHAS loci at PPR genes in all Solanaceous species, although some of the Solanaceous species (petunia and Nicotiana) harbor few PHAS loci related to NB-LRR defense genes (Baldrich et al., 2021). Instead of protein-coding PHAS loci triggered by miR482/2118 members, we found that these miRNAs initiate the production of reproductive tissue-enriched 21-nt phasiRNAs from noncoding loci in flax, a eudicot (Pokhrel, Huang, Bélanger, et al., 2021;and this study). We identified several miRNAs and phasiRNA pathways in flax that are responsive to heat stress. The two members of miR390, the conserved trigger of TAS3, were downregulated under heat stress, although the abundance of TAS3 derived phasiRNAs did not change significantly. We found no miR167 members to be differentially accumulated under heat stress, but we observed a reduction in phasiRNAs from miR167-triggered metal T A B L E 2 Differentially accumulated microRNA (miRNA) precursors in flax due to heat treatments Note: miRNAs that were significantly accumulated in different tissues were denoted as downregulated (D) and upregulated (U) (p < .05, FDR < .05, Table S4). B_2_4 indicates bud stage 1 with buds 2-4 mm in length, B_5_7 indicates bud stage 2 with buds 5-7 mm in length, and B_8_10 indicates bud stage 3 with buds 8-10 mm in length. "_S" and "_F" in the codes indicate 7 and 14 days of heat treatment. Note: Detailed information for these loci can be found in Table S5. Abbreviations: AFB2, AUXIN SIGNALING F-BOX 2; AGO, Argonaute protein; microRNA, miRNA; NAC, NAC domain-containing protein; N.D., not determined.

T A B L E 3 Summary of 21-PHAS loci and their triggers identified in flax
tolerance proteins phasing pathways that were enriched in reproductive tissues.
In recent studies, conservation of the miR2275-24-PHAS pathway in angiosperms was reported, including several eudicots Xia et al., 2019); the pathway is also conserved in flax, but with few (only four) 24-PHAS loci. Here, we identified a larger number (658) of 24-PHAS-like loci that are reproductive tissue-enriched but that lack obvious trigger miRNAs, possibly similar to the Solanaceous 24-nt phasiRNAs (Xia et al., 2019) or possibly some sort of loci of reproductive-tissue-enriched hc-siRNAs. We found that all 24-PHAS and 450 24-PHAS-like loci were heat responsive and all are downregulated due to heat stress.
Notably, we observed downregulation of the reproductive phasiRNA triggers, miR482/2118 and miR2275 in response to heat treatments, as well as the reproductive 21-and 24-nt phasiRNAs triggered by these miRNAs, respectively. In rice, mutations in the precursors of two 21-PHAS loci resulted in temperature-sensitive genic male sterility (Fan et al., 2016;Zhou et al., 2012), although in maize, knockdown of DCL5, a gene responsible for the production of reproductive 24-nt phasiRNAs, greatly reduced the phasiRNAs and resulted in temperature-sensitive male sterile phenotype (Teng et al., 2020).
Cross et al. reported that two heat treatments reduce boll formation and seed set in flax (Cross et al., 2003), the same conditions that we used; however, the molecular mechanisms underlying these results

| Data analysis and visualization
Preprocessing and filtering of sRNA libraries were conducted as previously described (Pokhrel, Huang, Bélanger, et al., 2021). Two miRNA identification tools: ShortStack and miR-PREFeR were used (Johnson et al., 2016;Lei & Sun, 2014) for genome-wide MIR gene identification. To increase the sensitivity of miRNA detection (a result of increased supporting data), in ShortStack, two libraries were used at a time using default parameters except Àmismatch 0 and Àmincov 10.
T A B L E 5 Differentially accumulated 24-PHAS loci in flax due to heat treatments 0D, 0 U Note: 24-PHAS loci that were significantly accumulated in different tissues were denoted as downregulated (D) and upregulated (U) (p < .05, FDR < .05, Table S8). B_2_4 indicates bud stage 1 with buds 2-4 mm in length, B_5_7 indicates bud stage 2 with buds 5-7 mm in length, and B_8_10 indicates bud stage 3 with buds 8-10 mm in length. "_S" and "_F" in the codes indicate 7 and 14 days of heat treatment.
24-PHAS loci with unknown triggers were annotated as 24-PHAS-like loci. Target prediction was carried out by sPARTA (Kakrana et al., 2014) with a target cut off score of ≤5. The sRNA abundance and phasing score were viewed at a customized browser (Nakano et al., 2020). The 250 base pairs were padded in both upstream and downstream of predicted PHAS loci and subsequently used as the de novo annotated set for analysis by BLASTX (Altschul et al., 1997) against UniRef90. The summed abundance of sRNAs from given locus was calculated by ShortStack using -locifile parameter. Differential analysis (p < .05, false discovery rate [FDR] < .05) of sRNA loci was conducted in R using the edgeR package (Robinson et al., 2010). Circular plots were constructed using OmicCircos (Hu et al., 2014) for the chromosomal distributions of miRNAs and phasiRNAs generating loci. pheatmap (https://rdrr.io/cran/pheatmap/) was used to draw heatmaps in R.

ACKNOWLEDGMENTS
We thank members of the Meyers lab for helpful discussions and Joanna Friesner for assistance with editing. We thank Mayumi Nakano for assistance with data handling. This work was supported by resources from the Donald Danforth Plant Science Center and the

University of Missouri-Columbia plus the William H. Danforth Plant
Science Fellow Award (to S.P.).

CONFLICT OF INTEREST
The authors declare no conflict of interest associated with the work described in this manuscript.

AUTHOR CONTRIBUTIONS
B.C.M. and S.P. designed the research; S.P. performed the experiments and analyzed the data; S.P. and B.C.M wrote and revised the manuscript.

DATA AVAILABILITY STATEMENT
All the sRNA data generated during this study were deposited in NCBI's SRA (Sequence Read Archive) under the accession numbers PRJNA669702and PRJNA768110.