A lab-cultured pure line PH-40 of P. haitanensis was used in the experiments to eliminate the interference caused by genotypic differences. The sporophytes (conchocelis) of this pure line were directly developed from a single somatic cell that was isolated using 2% snail enzyme from a farmed gametophyte collected from Putian, Fujian Province, China. The gametophytes of P. haitanensis were cultured in a light incubator under the following conditions: 20 ± 1°C with 50–60 μmol photons⋅m^–2⋅s^–1 of illumination during a 12: 12 h light: dark cycle. Provasoli-Enriched Seawater (PES) medium was changed every 3 days. Material processing was performed when the blade was 10–15 cm long. For miRNA identification, gametophytes were exposed to dehydration and rehydration. The gametophytes were harvested in the following conditions: normal conditions as control (CON), 30 ± 5% water loss rate (moderate water loss, MWL), 50% ± 5% (high water loss, HML), and 80% ± 5% (severe water loss, SWL). The gametophytes from SWL were then put back into seawater to rehydrate for 30 min (REH) according to a previous study (Wang et al., 2015). The different water loss levels were determined according to a previous study (Kim et al., 2009). Three biological replicates were performed for each treatment. After treatment, all samples were stored in liquid nitrogen for subsequent RNA extraction.

Total RNA was separately isolated from gametophytes (control and osmotic-treated samples) using Trizol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s protocol. Concentrations and qualities of the isolated RNAs were measured using a NanoDrop spectrophotometer (Thermo Fisher Scientific). Approximately 2.5 μg of total RNA was used to prepare a small RNA library, according to the protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, United States). Then, the libraries were sequenced by Illumina Hiseq2500 50SE (single-end) at LC-Bio (Hangzhou, China).

The raw reads were processed using Illumina’s Genome Analyzer Pipeline software (Solexa 0.3) and the ACGT V3.1 program developed by LC Sciences (LC Sciences, Houston, TX, United States). After removing adapter dimers, junk, and low complexity sequences from the small RNA reads, total read counts for the small RNAs were established. Then, small RNAs that mapped to common RNA families [ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA)] were discarded. Subsequently, the remaining sequences were mapped to the assembled genome of P. haitanensis using the Bowtie software (Langmead et al., 2009; Cao et al., 2019). For alignment, one mismatch in the first 16 nt was allowed. Mapped reads were used for further analysis. The mapped sequences that were 18–25 nt long were mapped to specific species precursors in miRBase 21.0 by Bowtie search to identify known miRNAs and novel 3p- and 5p- derived miRNAs. Length variation at 3′ and 5′ ends and one mismatch inside of the sequence were allowed in the alignment. Meanwhile, the remaining sequences that did not match known miRNAs were used to identify novel miRNAs in P. haitanensis.

miRNA differential expression based on normalized deep sequencing counts was analyzed by selectively using DeSeq. The miRNA expression level fold-change was calculated using the formula mentioned in a previous report (Marsit et al., 2006): Fold-change = log₂ of RPKM value. The significance threshold (p-value) was set at 0.01 and 0.05 for each test. The fold-change (>1) and p-values (p < 0.05) were combined to identify differentially expressed genes, which were defined as osmotic stress-responsive miRNAs.

To further understand the function of the differentially expressed miRNAs, computational target prediction algorithms (TargetScan5.0 and miRanda 3.3a) were used to predict the miRNA binding sites (potential target genes) (Yang et al., 2012). The function of the target genes of these differentially expressed miRNAs were annotated using the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways² databases.

To validate the osmotic stress-responsive miRNAs identified by high-throughput sequencing and expression trends, the differentially expressed miRNA and its target gene was confirmed by qPCR (PC-3p-99769_194-ph08980.t1). The gametophytes were collected as described for miRNA library preparation and sequencing. Then, small RNAs (<200 nt) were extracted from the samples following the manufacturer’s instructions for the RNA isolation for Small RNA Kit (TaKaRa). Equal amounts of high-quality small RNAs were reverse transcribed to complementary DNA (cDNA). Reverse transcription (RT) was performed using the PrimeScript RT reagent Kit (TaKaRa, China) according to the manufacturer’s instructions but with specific stem-loop RT primers for miRNAs (Supplementary Table S1). The 5.8S and 18S subunits were used as internal controls for miRNA and mRNAs, respectively. The expression levels of these genes were analyzed by qRT-PCR with a Light-Cycle^® 480 Real-Time PCR System. The reactions were performed using 20 μL volumes containing 10 μL of SYBR^® Premix Ex Taq (TaKaRa Biotech Co., Dalian, China), 0.6 μL of each primer (forward and reverse primer), 2 μL of diluted cDNA, and 6.8 μL of RNA-free water using the following protocol: 95°C for 30 s, 40 cycles at 95°C for 5 s, 57°C for 30 s, and 72°C for 30 s. The relative expression levels were calculated using the comparative Ct (2^–ΔΔCt) method.

To profile small RNAs and identify the osmotic stress-responsive miRNAs in P. haitanensis, small RNA sequencing for Pyropia gametophytes subjected to varying degrees of dehydration and rehydration (CON, MWL, HWL, SWL, and REH) were performed. After filtering low-quality raw reads and mRNA, Rfam, and Repbase mappable reads, the average reads of these samples were 2,507,531, 2,705,002, 2,620,137, 2,741,703, and 2,842,799 bp, respectively (Table 1). These reads were used for further miRNA identification and prediction.

Small RNA sequences were mapped to miRBase and the P. haitanensis genome to identify conserved and novel miRNAs. Forty-seven miRNAs were predicted and 95.74% of these were 20–24 nt long, with the 21 nt as the dominant length (Figure 1). In the first nucleotide bias analysis of miRNAs in P. haitanensis, the first residue of the 19 nt miRNAs was either predominantly uridine (U) or cytosine (C). The first residue of the 20, 21, and 22 nt miRNAs was predominantly U, whereas the 24 nt miRNAs predominantly started with guanine (G).

Among the 47 miRNAs, 2 could be mapped to specific miRNA sequences in miRBase and the extended Pyropia genome sequences from their genomic loci could form hairpins. Thus, they were defined as conserved miRNAs for this study. The other 45 are novel miRNAs (Table 2).

To trace the evolutionary scenario of miRNAs in the red algae lineage, we downloaded the published miRNA data for the Pyropia species P. yezoensis and the sister red algal clade Floriceadens Chondrus crispus for phylogenetic comparative analysis. We discovered that the pre-miRNA sequences of 15 miRNAs could be aligned to both P. yezoensis and Porphyra umbilicalis, which indicated that these sequences were well conserved in Pyropia species. Only two were also observed in C. crispus (Supplementary Table S2). The target genes of these conserved miRNAs are mainly involved in DNA synthesis (DNA repair proteins and DNA mismatch repair proteins), signal pathways (mitogen-activated protein kinase and serine/threonine protein kinase), transporters (ABC transport system and zinc transporter proteins), and some transcriptional regulators.

Many genes undergo significant expressional variation in response to osmotic stress (Wang et al., 2015). Our identification of miRNAs in the P. haitanensis genome caused us to investigate whether miRNAs contribute to the post-transcriptional regulation of key genes that lead to physiological modulation to adapt to stressful conditions. Therefore, we further investigated the expression levels of miRNAs under different degrees of dehydration and subsequent rehydration. Twelve miRNAs exhibited significant variation in at least one dehydrated or rehydrated condition (Figure 2). According to their transcriptional patterns along the time course of osmotic stresses, they can be classified into three clusters (Figure 2). In cluster I, three miRNAs exhibited slight variations in response to most of the stress conditions, while the levels of PC-3P-45934_389 and PC-3P-75147_255 were transiently down-regulated in HWL and SWL, respectively. There were four miRNAs in cluster II. Although there was no significant change in their expression under dehydrated conditions, they were induced when Pyropia thalli were rehydrated after severe water loss. The final five miRNAs were grouped in cluster III. All of which showed a sharp increase in expression when Pyropia underwent water loss, except for PC-5P-180717_97, which was upregulated only until HWL. When rehydrated, the levels of most of the miRNAs in cluster III barely changed, and only PC-5p-210644_78 showed marked downregulation.

To understand the specific functions of these differentially expressed miRNAs under osmotic stress in P. haitanensis, we predicted the target genes of differentially expressed miRNAs and analyzed the enriched functional categories that they encoded in terms of GO and KEGG pathways. Overall, 60 target genes were predicted for 12 miRNAs. According to the GO classification system, these genes were classified into molecular function categories related to catalytic activity, binding, transporter activity, transcription factor activity, signal transducer activity, and antioxidant activity (Supplementary Table S3). The results of KEGG metabolic pathway annotation showed that cell cycle, meiosis, nucleotide excision repair, DNA replication, glycerolipid metabolism, and base excision repair pathways were significantly enriched (P < 0.01) (Supplementary Table S3). Nucleotide excision repair, base excision repair, cell cycle, and DNA replication were also significantly enriched in differentially transcribed genes under osmotic stress in P. haitanensis as previously described (Wang et al., 2015). The consistency between the two datasets supports the regulatory role of miRNAs in response to osmotic stress.

miRNA PC-3p-99769_194 from cluster III was predicted to bind to the 523–542 bp region of the CA gene and may regulate its expression by degrading mRNA. The miRNA PC-3p-99769_194 was upregulated during the whole dehydrated process. Accordingly, the target CA gene had declined in transcriptional activity as observed in the published transcriptomic data under dehydration in P. haitanensis. The CA gene exhibited 0.77- and 0.44-fold downregulation under 30 and 80% dehydration and 1.3-fold upregulation under rehydration conditions.

PC-5p-335121_37 and PC-5p-350045_35 were highly expressed in samples under dehydration stress. The target genes of PC-5p-335121_37 are ph07921.t1 and ph07931.t1 which encode the U5 snRNP complex subunit and superoxide dismutase, respectively. PC-5p-350045_35 targets cid-thiol ligase and glycosyltransferase (ph09862.t1 and ph00598.t1), respectively. PTC-MIR482D-P3_2SS14CG18CG exhibited upregulated expression in rehydration samples. The target gene of ptc-MIR482D-P3_2SS14CG18CG is ph04985.t1, which encodes epoxyqueuosine reductase (Supplementary Table S3).

To validate the expression of miRNAs obtained from the high throughput sequencing approach, qPCR was performed on three miRNAs among samples of CON, MWL, HWL, SWL, and REH. PC-3p-99769_194 exhibited 2. 21-, 3. 46-, 1. 71-, and 1.83-fold upregulation in P. haitanensis after MWL, HWL, SWL, and REH, respectively. The expression of PC-5p-350045_35 in dehydration and rehydration samples were higher than those in the control (Figure 3). Its target gene CA (ph0980.t1) was downregulated 0. 87-, 0. 84-, and 0. 51-fold during MWL, HWL, SWL, respectively. And underwent a slight 1.39-fold increase during REH. These inverse transcriptional profiles suggest that PC-3p-99769_194A might function as a negative regulator of CA expression.

To unravel the RNA silencing machinery for miRNA generation and function in Pyropia, we identified genes encoding the two proteins that mediate the central activities required for miRNA and small RNA biogenesis and functioning, Dicer and Argonaute (AGO), respectively. Dicer proteins are characterized in many organisms, have multiple functional domains from the N- to C-terminus, and are approximately 2000 amino acids (aa) long. In P. haitanensis, we performed a Basic Local Alignment Search Tool (BLAST) search against the genome using Dicer proteins in model organisms as the query and discovered two candidate genes that exhibited high sequence similarity (ph07342.t1 and ph05047.t1 at 562 aa and 1163 aa long, respectively). Ph07342.t1 contains a DEAD-like helicase domain (DEXDc) and a helicase C (HELIc) domain and was aligned to the N-terminal part of reference Dicer. Ph05047.t1 has two parallel RIBO (ribonuclease III) domains. The two genes are located in two separate loci and are not adjacent to each other, which is different from the gene arrangement in other species (Figure 4). This phenomenon was verified in Po. umbilicalis. However, red macroalgae C. crispus and Gracilariopsis chorda, as well as the unicellular red algae Porphyridium purpureum, have one full-length canonical Dicer gene with multiple functional domains. To trace the evolutionary history of the two segmented Dicer genes, we identified Dicer homologs in available species and performed phylogenetic analyses. In the evolutionary trees constructed based on model Dicer sequences for ph07342.t1 and ph05047.t1, both closely clustered with the two Po. umbilicalis proteins and formed a unique branch, whereas the Dicer proteins from other species united as a single group (Supplementary Figure S1). For AGO proteins, there were 1, 2, 1, 3, and 3 candidate genes found in P. haitanensis, Po. umbilicalis, G. chorda, C. crispus, and P. purpureum genomes, respectively. All had canonical structures with an AGO-specific domain, a PAZ domain, and a PIWI domain at the C-terminal. P. haitanensis AGO exhibited close phylogeny with homologs in other red algae (Supplementary Figure S2). Conserved amino acids in the PIWI domain that facilitate siRNA binding were observed in P. haitanensis AGO (Figure 5). Moreover, it harbored a classical glutamate as a key residue of the catalytic DDE motif in RNase H. The DDH motif is required for metal ion coordination in the catalytic center, albeit “H” residue was not found in the corresponding locus (Figure 5). The identification of the key components of miRNA biogenesis machinery in P. haitanensis, as well as their expressional variation under desiccation, indicated the potential involvement of an miRNA-based regulatory mechanism in response to osmotic stress and gene engineering via artificial miRNA and siRNA. Although AGO that exhibited contractive gene copies, albeit well-conserved sequences, the distinctive structure of Dicer proteins lets us envisage some unique features possibly existing in Pyropia and Porphyra. In addition, we searched the transcriptional levels of two Dicer genes and the AGO gene in the published transcriptomic data under dehydration in P. haitanensis and discovered that the two Dicer genes were increased under dehydration and rehydration conditions. The AGO gene was 1.14-fold upregulated under medium dehydration, whereas AGO was 0.69- and 0.52-fold downregulated under severe dehydration and rehydration conditions.

The expression of miRNAs, an important class of gene regulators, is altered by abiotic stress treatment such as light, temperature, moisture, salt, and heavy metal ions (Jia et al., 2009; Lv et al., 2010). However, as most of these studies were performed using model organisms, little is known about seaweeds. CA catalyzes the interconversion between carbon dioxide and dissolved bicarbonate. It is the key component of the carbon concentration mechanism in aquatic photosynthetic organisms by facilitating the accumulation of sufficient carbon dioxide in the vicinity of rubisco (Badger and Price, 1994). However, the exact molecular mechanism involved in the control of the expression of CA is unclear. In this study, we demonstrated a post-transcriptional regulatory mechanism of CA and identified an miRNA as the key element involved in the regulation of the expression of the CA gene. The raised level of expression of PC-3p-99769_194 under dehydration stress suggested its increased negative regulation of the target CA gene. This is further supported by the decreased transcriptional level of CA gene in our previous study (Wang et al., 2015). Three enzymes involved in photosynthesis (rubisco, fumarate and nitrate reductase, and CA are downregulated in desiccated individuals compared to naturally hydrated controls by 2DE and LC-MS/MS analyses (López-Cristoffanini et al., 2015). Under high salinity treatment, the algae Nitzschia closterium f. minutissima adjusts its photosynthetic pigment content to decrease photosynthesis and inhibits the activity of CA (Yu et al., 2011). Combined with our study, we speculate that Pyropia reduce their photosynthetic activities in response to osmotic stress via miRNA-based regulatory mechanisms.

Rehydration can induce various metabolic processes that responding to osmotic stress in seaweed. The target gene of upregulated ptc-MIR482D-P3_2SS14CG18CG is related to epoxyqueuosine reductase, which catalyzes the last step of the two-electron reduction of epoxyqueuosine to queuosine during the synthesis of tRNA (Nishimura, 1983) and plays an important role in translational efficiency, accuracy, and structural stabilization of tRNA (El Yacoubi et al., 2012). In our study, genes encoding epoxyqueuosine reductase showed high expression when P. haitanensis was subjected to osmotic stress and its expression was suppressed during rehydration. Salt treatment in Anabaena fertilissima leads to the hyperaccumulation of epoxyqueuosine reductase (Rai and Swapnil, 2019). Therefore, we speculate that suppressed miRNA leads to the increased expression of epoxyqueuosine reductase in response to osmotic stress in Pyropia.

Dicer and AGO are the key components in miRNA biogenesis and are characterized in almost all organisms that have miRNA regulatory mechanisms (Bernstein et al., 2003; Yang et al., 2013). Both have well-conserved homologs with similarly organized functional domains and motifs. However, in this study, we identified that in Pyropia, Dicer was divided into two genes, ph07342.t1 that harbors the N-terminal functional domains and an open reading frame (ORF) that contains the other domains. The segmentation of Dicer was also observed in other Bangiales including P. yezoensis and Po. umbilicalis, whereas in the Florideophyceae genomes, such as C. crispus, G. chorda, and P. purpureum, full-length canonical Dicer was present. This particular distribution of segmented Dicer protein suggested that it might be generated after the evolutionary diversification of Bangiales. The length of a small RNA is mainly related to the type of enzymes involved in its processing. For example, a 21 nt segment is generally generated by Dicer-like 1 (DCL1), a 22 nt segment by RNase endonuclease, and a 24 nt segment by Dicer-like 3 (DCL3) (Xie et al., 2004; Vazquez, 2006; Ou et al., 2012). The characteristic length distribution of small RNAs, with 21 nt as the dominant length in Pyropia, suggests that Pyropia Dicer protein might function in a DCL1-like way in miRNA biogenesis.

Partial and truncated Dicer with a C-terminal fragment was observed in Caenorhabditis elegans before and characterized to be a negative regulator of miRNA biogenesis via its competitive interaction with AGO proteins (Sawh and Duchaine, 2013). However, the truncated small Dicer was generated by in vivo proteolytic cleavage, rather than by alternative splicing or individual “truncated” genes as in Pyropia. Thus, we envisage that they might function in a different way. The action mode is “two-peptides,” which means that the two peptides encoded by the two ORFs serve as the “subunits” and form a heterodimer. In Pyropia AGO, although a close phylogenetic relationship with its red algal counterparts was observed, it was encoded by only one gene, whereas three were observed in other red algae. Therefore, the segmented structure of Dicer and the contracted copy numbers of AGO represented the distinctive feature of miRNA biogenesis in Pyropia and imply a diverse miRNA regulatory mechanism. The latter was further reinforced by the dominant length of miRNA and a relative lack of universally conserved miRNAs.