Small RNAs, Degradome, and Transcriptome Sequencing Provide Insights into Papaya Fruit Ripening Regulated by 1-MCP

As an inhibitor of ethylene receptors, 1-methylcyclopropene (1-MCP) can delay the ripening of papaya. However, improper 1-MCP treatment will cause a rubbery texture in papaya. Understanding of the underlying mechanism is still lacking. In the present work, a comparative sRNA analysis was conducted after different 1-MCP treatments and identified a total of 213 miRNAs, of which 44 were known miRNAs and 169 were novel miRNAs in papaya. Comprehensive functional enrichment analysis indicated that plant hormone signal pathways play an important role in fruit ripening. Through the comparative analysis of sRNAs and transcriptome sequencing, a total of 11 miRNAs and 12 target genes were associated with the ethylene and auxin signaling pathways. A total of 1741 target genes of miRNAs were identified by degradome sequencing, and nine miRNAs and eight miRNAs were differentially expressed under the ethylene and auxin signaling pathways, respectively. The network regulation diagram of miRNAs and target genes during fruit ripening was drawn. The expression of 11 miRNAs and 12 target genes was verified by RT-qPCR. The target gene verification showed that cpa-miR390a and cpa-miR396 target CpARF19-like and CpERF RAP2-12-like, respectively, affecting the ethylene and auxin signaling pathways and, therefore, papaya ripening.


Introduction
Papaya (Carica papaya L.) is a popular and widely planted fruit in tropical and subtropical regions due to its unique taste, nutritional benefits, and medicinal benefits [1][2][3]. As one of the classic climacteric fruits, papaya fruits are highly perishable after harvest as a result of rapid ripening and softening and susceptibility to biotic or abiotic stresses [4,5], which usually result in a high percentage of product loss [6,7].
The ethylene antagonist 1-methylcyclopropene (1-MCP) has been widely used to maintain the quality and extend the shelf life of harvested products [8]. By binding to ethylene receptors, 1-MCP can inhibit the release of ethylene, inhibit the respiratory intensity of fruits and vegetables, and delay the ripening of various fruits and vegetables [9][10][11][12]. 1-MCP treatment on papaya at earlier development stages or for long-term/high-concentration treatments will lead to a "rubbery texture" phenomenon, where the fruits are unable to completely soften during later storage, leading to tasteless lack of flavor [13]. RNA-seq analysis has shown that improper 1-MCP treatment severely inhibits cell wall degradation and fruit softening by inhibiting ethylene signal transduction and cell wall metabolism pathways [10]. Integrated analysis of metabolomics and RNA-seq data has shown that various energy metabolites for the tricarboxylic acid cycle, glycolic acid cycle, flavonoids, and were disease-free, packaged into polystyrene boxes, and then transported immediately to the laboratory. The harvested fruit were dipped in 0.2% (w/v) hypochlorite solution for 10 min and then soaked in a 500 mg·mL −1 solution of mixed prochloraz (Huifeng, Jiangsu, China) and iprodione (Kuaida, Jiangsu, China) for 1 min to minimize the effect of microbes/microbe contamination. Treatments were then carried out with the same protocol as in our previous work in Zhu et al. [31]. A total of 600 papaya fruits were selected and divided into three sub-groups, and each group contained 200 papaya fruits. Two groups of the fruit were treated with 400 nL·L −1 1-MCP in a closed foam box for 2 and 16 h, respectively; air treatment for 16 h was used as a control. All fruits were then treated with 1000 µL·L −1 ethephon for 1 min. After air-drying for a few minutes, the papaya fruit was packed in an unsealed polyethylene bag (10 cm × 20 cm, 0.02 mm thick) and stored at 22 • C. All the treatments were conducted with three biological replicates.
Three sampling points at days 0, 1, and 6 were selected for measurement with three biological replicates. For each sampling point, three fruits representing the three biological replicates were collected. Samples were taken from the middle of the papaya fruit flesh, frozen immediately using liquid nitrogen, and then stored at −80 • C for further testing. Subsequent usage included RNA extraction, small RNA sequencing, and degradome sequencing.

RNA Extraction and cDNA Library Preparation
For RNA-seq, all RNA extraction and library preparation procedures were conducted as described in Zhu et al. [31].
The NEB Next Ultra-small RNA Library Preparation Kit for Illumina (NEB, #E7530, Ipswich, MA, USA) was used according to the manufacturer's protocol in order to generate a small RNA-seq library. Library quality was evaluated using an Agilent Bioanalyzer 2100 system. Samples from the control group, short-term 1-MCP treatment group, and longterm 1-MCP treatment group were selected for miRNA analysis after storage for 0, 1, and 6 days. Each sample time contained three biological replicates, and a total of 21 libraries were constructed. After removing the adapter sequences and low-quality sequence reads (including reads containing more than 10% N, reads without 3 linker sequences, and sequences shorter than 18 nucleotides or longer than 30 nucleotides), the clean reads were then mapped to the papaya reference genome (http://www.plantgdb.org/CpGDB/). Using a BLASTN search against miRbase (V21), known miRNAs were identified from mapped small RNA tags. For the new miRNA candidates, we used miRDeep2 for their identification. The total number of identified miRNAs in each constructed library was then normalized to TPM (number of transcripts per million of the clean tags). The DEGseq R package was used for differential expression analysis between the different groups. The q-value was used to adjust the p-value. Small RNAs with |log2(foldchange)| > 1 and p -value < 0.01 were assigned as differentially expressed. Target Finder software was used to predict the target genes of differentially expressed miRNAs (http://targetfinder.org/). The target gene identification was also based on the integration of small RNAs and transcriptome sequencing analysis. Normally, miRNA and mRNA pairs have a negative regulatory network relationship regarding their expression. GO software (http://www.geneontology.org/) was used to enrich and analyze the differential target genes between sample groups, and the KEGG (Kyoto Encyclopedia of Genes and Genomes) database was used to analyze the pathway annotation of differentially expressed miRNA target genes (http://www.genome.jp/kegg/).
The papaya samples, including fruits treated with 400 nL· L −1 1-MCP for 2 or 16 h and from the control group, were mixed and used to establish a degradation group library. After the samples were extracted, total RNA was extracted using the reagent Trizol (Invitrogen, CA, USA). The quantity and purity of the total RNA were analyzed using a Bioanalyzer 2100 and RNA 6000 Nano Lab Chip Kit (Agilent, CA, USA) with an RIN number > 7.0. Approximately 20 µg of total RNA was used for degradome library construction. Single-end sequencing (36 bp) was performed on an Illumina Hiseq2500 at LC-BIO (Hangzhou, China) according to the protocol recommended by the sequencing facility. Raw sequencing reads containing no adaptors nor low-quality reads were obtained using Illumina's software. The filtered sequencing reads were then used to identify potentially cleaved targets using the CleaveL pipeline. Degradome data and reads were mapped to the mRNA downloaded from JGI (http://www.plantgdb.org/CpGDB/). Only perfectly matched alignments to the given read were preserved for further degradation analysis.

RT-qPCR Verification
In order to verify the expression profiles of the identified miRNA and its target mRNA, 7 known miRNAs, 4 novel miRNAs, and 12 corresponding target genes were selected for real-time quantitative PCR (RT-qPCR) validation.
Total RNA and small RNAs were extracted using TRIzol reagent and Fruit-mate for RNA purification (Takara, Japan). Total RNA was reverse-transcribed using stem-loop qRT-PCR [32]. Stem-loop primers for reverse transcription and primers for RT-qPCR are listed in Tables S1 and S2. The relative gene expression value was normalized against the relative value of the CpTBP1 gene for mRNA expression [33] and 5 s RNA for miRNA expression [34]. SuperScript III reverse transcriptase (Invitrogen, Carlsbad, California, USA) was used to reverse transcribe the total RNA using the pulse reverse transcription program. The RT-qPCRs were performed using a total volume of 20 µL containing 10 µL SYBR Mixture (Promega, Madison, Wisconsin, USA), 3 µL cDNA template, 0.5 mM primers, and 6 µL ddH2O. PCR was performed on a Bio-Rad CFX96 real-time PCR system using the qPCR Master Mix Kit (Promega, Madison, WI, USA). Three biological replicates were used to determine the expression of each gene, and expression was calculated using the 2 −∆∆CT method [33].

Plant Expression Vector Construction and Transformation
Validated miRNAs were cloned into a pBI121 binary vector driven by a cauliflower mosaic virus 35S promoter (CaMV35S). The miRNAs' corresponding target genes were cloned into a pMS4 vector containing green fluorescent protein (GFP). Vector constructions were conducted in reference to Wang's method [35]. Confirmation of the target gene site was identified using the website http://plantgrn.noble.org/psRNATarget/). The primers are listed in Table S3. The constructed vectors were then transformed into tobacco mediated by Agrobacterium tumefaciens (GV3101). Young tobacco leaves were cultivated for five weeks and used for injection of Agrobacterium tumefaciens. Transformed leaves were photographed with ultraviolet light after 36 h of infection.

Statistical Analysis
Each treatment was conducted in three independent biological replicates, and the variance of the collected data was analyzed. Statistical differences between groups were assessed based on Duncan's Multiple Range Test (DMRT) in SPSS 19.0 (IBM, Armonk, NY, USA). The graphics were drawn using Prism 8 software (GraphPad Inc., La Jolla, CA, USA) and Adobe Illustrator software (CC 2017; Adobe Inc., San Jose, CA, USA).

Construction and Sequencing of Small RNA Libraries
Our previous work showed that 1-MCP treatment (400 nL·L −1 , 2 h) can delay the ripening of papaya fruit, but improper 1-MCP treatment (400 nL·L −1 , 16 h) causes an adverse "rubbery" texture [31]. Therefore, small RNA (sRNA) sequencing was performed. A total of 354.06 megabytes (MB) of clean reads was obtained from 21 different libraries. For each sample, low-quality sequences and reads with an unknown base (N) greater than 10% of the sequence were first removed and then 3 adaptor sequences, reads with low complexity, and reads homologous to t/rRNA sequence were removed. The data for each sample consisted of no less than 12.23 MB of clean reads (Table S4). In addition, 98.77 to 99.11% of the total reads ranging from 18 to 30 nt library lengths exactly matched the papaya genome sequence. The length distribution of the sRNAs was mainly concentrated between  Figure 1A). However, in the papaya injected with papaya ringspot virus (PRSV), 21 and 24 nt length sRNAs accounted for the majority of the total sRNAs present [36]. One possible reason for the observed difference in mature sRNA lengths across the different groups may be due to the differential activity across the various sRNA biogenesis pathways. Nucleotide sequence analysis of the identified miRNAs showed that uridine (U) was the most common nucleotide in sRNAs of 20-24 nt length, while cytosine (C) had a higher enrichment in the sRNAs of 19 and 25 nt lengths ( Figure 1B). C and U account for the majority of the total number of bases in papaya sRNA.
sample, low-quality sequences and reads with an unknown base (N) greater than 10% of the sequence were first removed and then 3′ adaptor sequences, reads with low complexity, and reads homologous to t/rRNA sequence were removed. The data for each sample consisted of no less than 12.23 MB of clean reads (Table S4). In addition, 98.77 to 99.11% of the total reads ranging from 18 to 30 nt library lengths exactly matched the papaya genome sequence. The length distribution of the sRNAs was mainly concentrated between 18 and 24 nt, where the papaya libraries were mainly composed of 20 and 21 nt length sRNAs ( Figure 1A). However, in the papaya injected with papaya ringspot virus (PRSV), 21 and 24 nt length sRNAs accounted for the majority of the total sRNAs present [36]. One possible reason for the observed difference in mature sRNA lengths across the different groups may be due to the differential activity across the various sRNA biogenesis pathways. Nucleotide sequence analysis of the identified miRNAs showed that uridine (U) was the most common nucleotide in sRNAs of 20-24 nt length, while cytosine (C) had a higher enrichment in the sRNAs of 19 and 25 nt lengths ( Figure 1B). C and U account for the majority of the total number of bases in papaya sRNA.

Identification of Known and Novel miRNAs
Known miRNAs in papaya fruit were identified through sequence matching with an miRNA database (miRBase v21). Read sequences that were exactly the same as the known miRNA were considered to be an identified known miRNA. A total of 213 miRNAs were identified in the present study, 44 of which are already known miRNAs, and 169 are newly predicted miRNAs (Table S5). As shown in Figure 1C, the distribution of known and new

Identification of Known and Novel miRNAs
Known miRNAs in papaya fruit were identified through sequence matching with an miRNA database (miRBase v21). Read sequences that were exactly the same as the known miRNA were considered to be an identified known miRNA. A total of 213 miRNAs were identified in the present study, 44 of which are already known miRNAs, and 169 are newly predicted miRNAs (Table S5). As shown in Figure 1C, the distribution of known and new miRNAs in sRNAs is different across different lengths. The number of new miRNAs discovered was significantly greater than the number of known miRNAs. Among the newly discovered miRNAs, the number of 21 and 24 nt length miRNAs was significantly greater than that of those of other size classes. The distribution of known miRNAs and new miRNAs from each sample is shown in Figure 1D. Fewer miRNAs were identified in the 1 day (s) after treatment (DAT) and 6 DAT samples in the control than in the other samples. Among all the known miRNAs, 14 miRNAs showed high abundance > 5000 TPM, including cpa-miR159a, cpa-miR159b, cpa-miR160d, cpa-miR162a, cpa-miR164a, cpa-miR166a, cpa-miR166d, cpa-miR319, cpa-miR390a, cpa-miR396, cpa-miR477, cpa-miR8137, cpa-miR8141, and cpa-miR8148. Five identified miRNAs from across all the miRNA libraries showed low abundance < 100 TPM, including cpa-miR160c-3p, cpa-miR164d, cpa-miR5211, cpa-miR8134, and cpa-miR8150. In addition, 169 novel miRNA candidates from 21 different libraries were predicted. It seems that 25.4% of the novel miRNA candidates showed lower abundance < 500 TPM. Among all the identified novel miRNAs, unconservative_supercontig_106_26263 and unconservative_supercontig_17_8265 were in greatest abundance (Table S6).

Identification of Differentially Expressed miRNAs
Compared to 0 DAT, 60 differentially expressed miRNAs (DEMs) (|log2(FC)| ≥ 1, pvalue ≤ 0.05) were identified (at 1 DAT vs. 0 DAT and 6 DAT vs. 0 DAT) (Figure 2A). There were 46 DEMs identified between 1 and 0 DAT in the control group, of which 40 DEMs showed increased expression and 6 DEMs showed decreased expression ( Figure 2A). There were 20 DEMs identified between 6 and 0 DAT in the control group, among which 5 DEMs had increased expression and 15 DEMs had decreased expression. Among the DEMs in both comparison groups, there were six DEMs that were shared across both comparisons during the papaya ripening stage, of which four DEMs had increased expression and two DEMs had decreased expression during fruit ripening. The number of DEMs decreased during fruit ripening ( Figure 2A). The enriched GO terms for the predicted target genes are presented in Figure 2B. "Biological process of metabolic process", "single-organism process and cellular process", "the cellular component of cell part and cell", and "the molecular function of catalytic activity and binding" had the most significant differences in gene expression changes during papaya fruit ripening. In the KEGG enrichment analysis, there were 36 top enrichments for metabolic/biological pathways ( Figure 2C). Among these enriched genes, "metabolic pathways of plant hormone signal transduction", "fatty acid metabolism", "biosynthesis of amino acids", and "starch and sucrose metabolism" are important for fruit ripening. Figure 3A shows a total of 60 miRNAs that were differentially expressed between 1-MCP treatments (long-term and short-term) compared to the control ( Figure 3A). The amount of DEMs decreased with the storage time and duration of the 1-MCP treatment. Compared to the control, a total of 45 miRNAs were differentially expressed during shortterm 1-MCP treatment. There were 35 DEMs specifically identified in the comparison of 1-MCP 2 h vs. CK at 1 DAT, 5 DEMs specifically identified in the comparison of 1-MCP 2 h vs. CK at 6 DAT, and 5 common DEMs shared across both DATs. Contrary to this, a total of 27 DEMs were identified in the comparison of long-term 1-MCP vs. the control. Of these, 11 and 12 DEMs were specifically identified after 1 and 6 DAT, respectively, and 4 DEMs were shared across both DATs. These results indicate that miRNAs are important in regulating fruit ripening.
In the abnormal fruit ripening case caused by 1-MCP treatment, the number of DEMs decreased dramatically compared to short-term 1-MCP treatment, indicating that the missing DEMs may be involved in the softening disorder. In order to fully understand the functions of DEM target genes, a GO enrichment analysis was performed ( Figure 3B). The enrichments for "metabolic process and cellular process", "the cellular component of cell part and cell", and "molecular function of catalytic activity and binding" biological processes showed the most significant differences in gene expression during 1-MCP-treated papaya fruit ripening. From the KEGG enrichment terms, there were 38 top enrichments for metabolic/biological pathways ( Figure 3C). Among these KEGG enrichment terms, "metabolic pathways of plant hormone signal transduction", "phenylpropanoid biosynthesis", "protein processing in endoplasmic reticulum", "carbon metabolism", and "starch and sucrose metabolism" genes were the most enriched, indicating the important roles of these terms in fruit ripening and 1-MCP-caused ripening disorder. The COG function for different treatments is shown in Figure S1. All DEMs were clustered into 26 functional categories based on four kinds of differential enrichment analyses. Among the 26 functional categories, "general function prediction only", "transcription", and "replication, recombination, and repair" were the most enriched items. The signal transduction mechanism was also enriched, which is an important functional protein for the ripening process of papaya ( Figure S1).  Figure 3A shows a total of 60 miRNAs that were differentially expressed between 1-MCP treatments (long-term and short-term) compared to the control ( Figure 3A). The amount of DEMs decreased with the storage time and duration of the 1-MCP treatment. Compared to the control, a total of 45 miRNAs were differentially expressed during shortterm 1-MCP treatment. There were 35 DEMs specifically identified in the comparison of 1-MCP 2 h vs. CK at 1 DAT, 5 DEMs specifically identified in the comparison of 1-MCP 2 h vs. CK at 6 DAT, and 5 common DEMs shared across both DATs. Contrary to this, a total of 27 DEMs were identified in the comparison of long-term 1-MCP vs. the control. Of these, 11 and 12 DEMs were specifically identified after 1 and 6 DAT, respectively, and   In the abnormal fruit ripening case caused by 1-MCP treatment, the number of DEMs decreased dramatically compared to short-term 1-MCP treatment, indicating that the missing DEMs may be involved in the softening disorder. In order to fully understand the functions of DEM target genes, a GO enrichment analysis was performed ( Figure 3B). The Ethylene and auxin signaling pathways showed an important role in the ripening of papaya treated with different durations of 1-MCP treatment [31]. According to the perfect or near-perfect complementarity of expression between miRNAs and their targets, DEMs were predicted to correspond to uni-genes as potential targets. Through sRNA sequencing analysis, Target Finder software was used for predicting targets. According to sRNA sequencing, miRNAs and their predicted target genes related to ethylene and auxin were found, and eight miRNAs and their predicted target genes with opposite expression

Combined Expression Analysis of miRNAs and Their Target mRNAs
Through differential expression analysis of the miRNAs and mRNAs in the sRNAs and transcriptome sequencing data, key miRNAs and genes were identified. The regulation of miRNAs and their target genes was uncovered. Compared to 0 DAT, 34 DEMs were identified, of which 20 DEMs were found at both 1 and 6 DAT. Six DEMs were shared across both DAT groups ( Figure 4A). The enriched target gene GO terms are presented in Figure 4B. The enrichments for "biological process of metabolic process and cellular process" and "the molecular function of binding" biological processes showed the most significant differences in gene expression during papaya fruit ripening. From the KEGG enriched terms, "metabolic pathways of plant hormone signal transduction", "carbon metabolism", and "starch and sucrose metabolism" were important for fruit ripening ( Figure 4C).
Combined with transcriptome analysis, a total of 34 miRNAs were expressed in the different 1-MCP treatments (long-term and short-term). The amount of DEMs increased with treatment storage time and duration of the 1-MCP treatment. In total, 22 DEMs were identified in the comparison of short-term 1-MCP vs. the control at 6 and 1 DAT. Among the 22 DEMs, 17 DEMs were identified at 1 DAT, 10 were at identified 6 DAT, and 5 were shared across both dates. The enriched GO terms for miRNA target genes are presented in Figure 5B. The enrichments for "metabolic process and cellular process" and "molecular function of binding" biological processes showed the most significant differences in gene expression during 1-MCP-delayed papaya ripening. From the KEGG enriched terms, "metabolic pathways of plant hormone signal transduction" and "starch and sucrose metabolism" were important for fruit ripening ( Figure 5C).
Combined with transcriptome analysis, the miRNA target genes were screened out if they had an opposite expression pattern to their respective miRNAs. In total, 12 miR-NAs were found to be involved in the signal transduction process of ethylene and auxin. The miRNAs and their corresponding target genes are shown in Table 2    in Figure 5B. The enrichments for "metabolic process and cellular process" and "molecular function of binding" biological processes showed the most significant differences in gene expression during 1-MCP-delayed papaya ripening. From the KEGG enriched terms, "metabolic pathways of plant hormone signal transduction" and "starch and sucrose metabolism" were important for fruit ripening ( Figure 5C).

Target Gene Identification of Papaya miRNAs by Degradome Analysis
Through degradome sequencing, the splice sites of miRNAs in the mRNA were found. The 1741 target genes corresponding to their miRNAs are summarized in Table S7. According to KEGG classification, the target genes of these key DEMs related to plant hormone signal transduction are summarized and shown in Figure 6A,B. The plant hormone signal pathways include ethylene, jasmonic acid, auxin, gibberella, MAPK, and abscisic acid. Among all these pathways, the miRNAs related specifically to the ethylene and auxin signaling pathways were the most enriched. The corresponding key target genes were CpCTR1, CpARFs, CpTIR1, and CpSAUR67. The expression of some miRNAs increased with fruit ripening, such as unconservative_supercontig_106_26263, cpa-miR156a, and cpa-miR160d and cpa-miR160a. Some miRNAs were closely related to papaya softening disorder, such as unconservative_supercontig_20_9672, unconservative_supercontig_13_6816, cpa-miR8148, unconservative_supercontig_152_30952, and unconservative_supercontig_119_27967. Four miRNA targets to CpCTR1 were identified: unconservative_supercontig_75_21810, unconservative_supercontig_2_1866, unconservative_supercontig_19_9148, and unconserva-tive_supercontig_3_2646. CpSAUR67 and CpTIR1 were involved in auxin biosynthesis and signal transduction and were putative targets for unconservative_supercontig_20_9672 and cpa-miR393, respectively. CpARFs were putative targets for cpa-miR160d, cpa-miR319, unconservative_supercontig_184_32956, cpa-miR156a, and cpa-miR160a. The expression of these target genes was neither increased with fruit ripening and repressed by 1-MCP treatment nor decreased with fruit ripening, but it was induced by 1-MCP treatments ( Figure 6B), indicating that these target genes may play an important role in the fruit ripening process.

Target Gene Identification of Papaya miRNAs by Degradome Analysis
According to the analysis of small RNAs, degradome, and transcriptome sequencing, the miRNAs related to ethylene and auxin signaling were selected and visualized. The cpa-miR160a and cpa-miR160d miRNAs co-targeted CpARF10, CpARF16-like, and CpARF17, thereby affecting auxin signaling and fruit ripening. The cpa-miR172a miRNA targeted CpERF RAP2-7 and Carica papaya auxin transport protein BIG, and unconserva-tive_supercontig_2_1866 targeted CpARF3 and CpCTR1, thereby affecting the signaling pathways of ethylene and auxin simultaneously (Figure 7).

Target Gene Identification of Papaya miRNAs by Degradome Analysis
According to the analysis of small RNAs, degradome, and transcriptome sequencing, the miRNAs related to ethylene and auxin signaling were selected and visualized. The cpa-miR160a and cpa-miR160d miRNAs co-targeted CpARF10, CpARF16-like, and CpARF17, thereby affecting auxin signaling and fruit ripening. The cpa-miR172a miRNA targeted CpERF RAP2-7 and Carica papaya auxin transport protein BIG, and unconserva-tive_supercontig_2_1866 targeted CpARF3 and CpCTR1, thereby affecting the signaling pathways of ethylene and auxin simultaneously (Figure 7).

Target Gene Identification of Papaya miRNAs by Degradome Analysis
Through the analysis of sRNAs and mRNA sequencing, the miRNAs and their target genes that are related to ethylene and auxin were selected and verified by RT-qPCR. There were four miRNAs related to ethylene signaling, namely cpa-miR172a, cpa-miR396, cpa-miR8140, and unconservative_supercontig_46_16464. Among these miRNAs, the expression of miR172a in the control group increased sharply at 1 DAT and then decreased. 1-MCP treatments severely repressed the expression of miR172a ( Figure 8A). CpERF RAP2-7-like was predicted to be the target gene of miR172a. The expression of CpERF RAP2-7-like decreased at 1 DAT and then increased at 6 DAT, with an opposite expression trend compared to miR172a ( Figure 8B). The expression of cpa-miR396, cpa-miR8140, and unconserva-tive_supercontig_46_16464 decreased with fruit ripening, while the expression of their corresponding target genes, namely CpERF RAP2-12-like, CpACS1-like, and CpACO4-like, increased first and then decreased, respectively ( Figure 8C-H). Among the four ethylenerelated genes, the expression level of CpACS1 significantly increased at 1 DAT in the control sample, signaling that these genes may play an important role in the early development stage of fruit ripening. The expression levels of CpACO4, CpERF4, and CpERF RAP2-

Target Gene Identification of Papaya miRNAs by Degradome Analysis
Through the analysis of sRNAs and mRNA sequencing, the miRNAs and their target genes that are related to ethylene and auxin were selected and verified by RT-qPCR. There were four miRNAs related to ethylene signaling, namely cpa-miR172a, cpa-miR396, cpa-miR8140, and unconservative_supercontig_46_16464. Among these miRNAs, the expression of miR172a in the control group increased sharply at 1 DAT and then decreased. 1-MCP treatments severely repressed the expression of miR172a ( Figure 8A). CpERF RAP2-7-like was predicted to be the target gene of miR172a. The expression of CpERF RAP2-7-like decreased at 1 DAT and then increased at 6 DAT, with an opposite expression trend compared to miR172a ( Figure 8B). The expression of cpa-miR396, cpa-miR8140, and uncon-servative_supercontig_46_16464 decreased with fruit ripening, while the expression of their corresponding target genes, namely CpERF RAP2-12-like, CpACS1-like, and CpACO4-like, increased first and then decreased, respectively ( Figure 8C-H). Among the four ethylenerelated genes, the expression level of CpACS1 significantly increased at 1 DAT in the control sample, signaling that these genes may play an important role in the early development stage of fruit ripening. The expression levels of CpACO4, CpERF4, and CpERF RAP2-7-like slightly decreased and were negatively correlated with fruit ripening. The expression of CpERF RAP2-12-like increased with fruit ripening and was positively correlated with fruit ripening (Figure 8). The RT-qPCR results were consistent with the RNA-seq data, and miRNA and the corresponding target genes showed opposite expression patterns.
Foods 2021, 10, x FOR PEER REVIEW 18 of 25 7-like slightly decreased and were negatively correlated with fruit ripening. The expression of CpERF RAP2-12-like increased with fruit ripening and was positively correlated with fruit ripening (Figure 8). The RT-qPCR results were consistent with the RNA-seq data, and miRNA and the corresponding target genes showed opposite expression patterns.

Target Gene Identification of Papaya miRNAs by Degradome Analysis
There were six miRNAs and seven target genes identified related to auxin signaling. Among the miRNAs and target genes, miR160a and miR160d worked together to regulate the expression of ARF10/16-like/17. Similar expression profiles in miR160d, miR167c, miR319, miR390a, and unconservative_supercontig_2_1866 were observed, which first increased and then decreased (Figure 9A,E,G,K). On the contrary, expression of the target gene ARF first decreased and then increased ( Figure 9B-D,F,H,L). The expression of un-conservative_supercontig_9_5033 first decreased and then increased; the expression of its target gene, auxin-binding protein T85, increased first and then decreased ( Figure 9I,J). The

Target Gene Identification of Papaya miRNAs by Degradome Analysis
There were six miRNAs and seven target genes identified related to auxin signaling. Among the miRNAs and target genes, miR160a and miR160d worked together to regulate the expression of ARF10/16-like/17. Similar expression profiles in miR160d, miR167c, miR319, miR390a, and unconservative_supercontig_2_1866 were observed, which first increased and then decreased (Figure 9A,E,G,K). On the contrary, expression of the target gene ARF first decreased and then increased ( Figure 9B-D,F,H,L). The expression of un-conservative_supercontig_9_5033 first decreased and then increased; the expression of its target gene, auxin-binding protein T85, increased first and then decreased ( Figure 9I,J). The expression of all six of these miRNAs was positively correlated with fruit ripening, but their targets showed a negative correlation with fruit ripening. 1-MCP treatment repressed the expression of these miRNAs but induced the expression of their target genes (Figure 9). expression of all six of these miRNAs was positively correlated with fruit ripening, but their targets showed a negative correlation with fruit ripening. 1-MCP treatment repressed the expression of these miRNAs but induced the expression of their target genes ( Figure 9).

Target Gene Identification of Papaya miRNAs by Degradome Analysis
According to the previous analysis, CpARF19-like and CpERF RAP2-12-like were potential targets of cpa-miR390a and cpa-miR396, respectively. To further test the relationship between these miRNAs and their target genes, transient co-expression assays in Nicotiana benthamiana leaves were conducted. As shown in Figure 10, the green fluorescent protein (GFP) fluorescence from tobacco leaf areas injected with 35s::Pre-cpa-miR390a + 35s::CpARF19-like-GFP and 35s::Pre-cpa-miR396 + 35s::CpERF RAP12-2-like-GFP was starkly lower than that of the negative controls ( Figure 10). These results provide evidence that the CpARF19-like and CpERF RAP2-12-like genes may be directly regulated by cpa-miR390 and cpa-miR396 miRNA, respectively.

Target Gene Identification of Papaya miRNAs by Degradome Analysis
According to the previous analysis, CpARF19-like and CpERF RAP2-12-like were potential targets of cpa-miR390a and cpa-miR396, respectively. To further test the relationship between these miRNAs and their target genes, transient co-expression assays in Nicotiana benthamiana leaves were conducted. As shown in Figure 10, the green fluorescent protein (GFP) fluorescence from tobacco leaf areas injected with 35s::Pre-cpa-miR390a + 35s::CpARF19-like-GFP and 35s::Pre-cpa-miR396 + 35s::CpERF RAP12-2-like-GFP was starkly lower than that of the negative controls ( Figure 10). These results provide evidence that the CpARF19-like and CpERF RAP2-12-like genes may be directly regulated by cpa-miR390 and cpa-miR396 miRNA, respectively.

Papaya miRNAs with Conserved and New Gene Targets
As an inhibitor of ethylene receptors, 1-MCP can effectively delay fruit ripening. This has wide applicability to post-harvest fruit preservation. Improper 1-MCP treatment can negatively affect the softening ability of fruit during ripening, resulting in papaya ripening disorder [37][38][39]. 1-MCP function can reduce the production of endogenous ethylene. In banana, 12 novel and 128 known miRNAs were identified. Among these miRNAs, 22 were differentially expressed after ethylene regulation treatment and were involved in the ethylene response [30]. In Wang's (2020) study between wild-type and LeERF1 transgenic tomato fruits, 9 miRNAs and 12 nat-siRNAs were found to be differentially expressed [40]. Degradome sequencing analysis further validated the nine miRNA targets, and six new targets were also identified. In our present study, 46 known miRNA families and 169 papaya-specific miRNAs were identified in 1-MCP-treated papaya samples using deep sequencing and computational analyses (Figure 1, Table S5). This difference in these two studies may be due to the type of fruit and the difference in post-harvest treatments. In papaya, a total of 1741 target genes of 178 known miRNAs were identified through degradome analysis (Table S7). Across the identified miRNAs, there were 40 known miRNAs, 138 novel miRNAs, and 29 miRNAs related to plant hormone signaling.

miRNAs Participate in Hormone Pathways during Papaya Fruit Ripening
Using miRNA sequencing analysis, 34 miRNAs related to papaya ripening and 60 miRNAs related to 1-MCP treatment were found (Figures 2 and 3). After target prediction (Table 1) and combined analysis of miRNAs and miRNA sequences, 12 miRNAs were found to directly act on ethylene-and auxin-related pathways. Of these 12 miRNAs, eight were known miRNAs and one was a novel miRNA. miRNAs are negative regulators of expression of their target genes. The miRNA mdm-miR160 targets the auxin response factors MdARF16 and MdARF17 and participates in the formation of adventitious root in apple rootstocks [27]. During papaya's ripening stage (6 DAT vs. 0 DAT in the control), the expression of cpa-miR160d decreased, while in papaya fruit with abnormal ripening (1-MCP 16 h 6 DAT vs. control 6 DAT), the expression of cpa-miR160d increased. These results indicate that cpa-miR160d is involved in the ripening process of papaya. The miRNA miR396 regulates plant growth and development by inhibiting the expression of growth-regulating factor (GRF) [41], participates in the development of wheat plants [42], and responds to environmental stressors such as cold injury, high temperature, and drought [43][44][45]. In our study, during normal fruit ripening (1 DAT vs. 0 DAT in the control), the expression of cpa-miR396 was found to be upregulated, while at the early development stage of 1-MCP-delayed ripening (1-MCP 2 h 1 DAT vs. control 1 DAT) in papaya, the expression of cpa-miR396 was downregulated. In Luan's (2018) study, by overexpressing miR172a/b in tomato plants, the expression of AP2/ERF transcription factor was inhibited, and the chlorophyll content and photosynthetic rate increased, as well as the development of higher resistance to phytophthora infection [46]. In the present work, the expression of cpa-miR172a showed no significant changes during normal ripening (1 DAT vs. 0 DAT in the control), while during the delayed papaya ripening process (1-MCP 2 h 6 DAT vs. 1 DAT) and the early development stage with abnormal ripening (1-MCP 16 h 1 DAT to control 1 DAT), the expression of cpa-miR172a was found to be upregulated. These results indicate that the identified miRNAs may play key roles in papaya fruit ripening.

Degradation Analysis Showed That miRNAs Are Involved in Regulation of Ethylene and Auxin Signaling Pathways
Through miRNA-targeted degradome analysis, miRNAs related to ethylene and auxin signaling were identified, as shown in Figure 6A,B. Their corresponding key target genes are included in Figure 6B: CpCTR1, CpARFs, CpTIR1, and CpSAUR67. The miRNAs and their corresponding target genes identified by degradome analysis are shown in Table S8. These miRNAs were found to regulate hormone signal transduction pathways and hormone homeostasis related to developmental processes. CTR1 is located downstream of the ethylene receptor and mediates the signal from the ethylene receptor by negatively regulating the ethylene response. Both miR1917 and miR171b target CTR1 and are important regulators for ethylene signal transduction in the ripening of tomato fruit [47,48]. In the present study, CTR1 was targeted by unconservative_supercontig_75_21810, unconservative_supercontig_2_1866, unconservative_supercontig_19_9148, and unconservative_supercontig_3_2646. The expres-sion patterns of unconservative_supercontig_75_21810, unconservative_supercontig_19_9148, and unconservative_supercontig_3_2646 decreased during fruit ripening. The expression of unconservative_supercontig_2_1866 decreased in 1-MCP-delayed fruit ripening.
At the early stages of auxin signal transduction, it was found that several gene families, including Aux/IAA (auxin/indole-3-acetic acid), SAUR (small auxin up RNA), and GH3 (Gretchen Hagen3), could respond to auxin treatments. The Aux/IAA protein acts as an important component of the auxin signaling pathway by participating in the regulation of expression for a large number of genes downstream of the auxin signaling pathway by releasing auxin response factor (ARF). Most of the Aux/IAA proteins contain four conserved domains, namely the I, II, III, and IV domains. Domain II is a key component that leads to instability of the Aux/IAA protein, which is later degraded by the ubiquitinproteasome protein (TIR1) pathway [49,50]. Therefore, TIR1, SAUR, and ARF are key factors in the auxin signaling pathway. In Arabidopsis thaliana, miR393 participates in auxin signal transduction and plant development by regulating TIR1 [50]. The miR165/166 miRNA determines the fate of A. thaliana root cells, is involved in plant hormone cross-talk, and regulates root growth via negative regulation of its target genes, such as the auxin response factors ARF10, ARF16, and ARF17 [51]. In the present study, ARFs were targeted by miR156a, miR160a, miR160d, miR319, and unconservative_supercontig_182_32956. The expression of these five miRNAs decreased during fruit ripening. TIR1 and SAUR67 were targeted by miR393 and unconservative_supercontig_65_20146, respectively. The expression of miR393 decreased while that of unconservative_supercontig_65_20146 increased during papaya fruit ripening.

Network Regulation Diagram of miRNAs and Target Gene Verification
According to the regulatory network diagram, RT-qPCR verification of key miRNAs and target genes was performed, showing that the expression of cpa-miR172a, cpa-miR319, cpa-miR390a, and cpa-miR396 greatly changed. Previous studies found that miR172amediated upregulation of ERF RAP2-7 during water submergence stress in maize roots restricted plant growth during flood-stress conditions [52]. In the present study, cpa-miR172a had an effect on ethylene and auxin signaling pathways. Here, cpa-miR396 acted on the ethylene signaling pathway, and cpa-miR319 and cpa-miR390a acted on the auxin signaling pathway. Similar results also showed that miR390 targeted ARFs and repressed their expression [53]. It has previously been reported that miR390 and ARFs form an auxin-responsive regulatory network to control lateral root growth in A. thaliana. The expression of miR390 is confined to the mesenchymal cells of the xylem prior to lateral root initiation. It was found that miR390 stimulates the production of tasiARFs, which repress the expression of their targets ARF3 and ARF4. There was also positive and negative feedback between ARF2/ARF3/ARF4 and miR390 to regulate lateral root growth [54]. Two highly expressed miRNAs were selected from the ethylene and auxin pathways, and their target genes were verified. We found that the expression signal of GFP protein with target gene binding in tobacco was weakened. These results indicate that CpARF19-like and CpERF RAP2-12-like are potential targets of cpa-miR390a and cpa-miR396, respectively, and they all may be important in the process of papaya fruit ripening.

Conclusions
Suitable 1-MCP treatment effectively delays the ripening of papaya fruit, and longterm 1-MCP treatment causes papaya ripening disorder. In this study, after different 1-MCP treatments on papaya fruits at different development stages, miRNA, transcriptome, and degradome analyses were performed. A total of 213 miRNAs and 1741 target genes of these miRNAs were identified. Among these, 11 different miRNAs related to ethylene and auxin and 12 corresponding target genes were found. The analysis and verification of cpa-miR390a and cpa-miR396, targeting CpARF19-like and CpERF RAP2-12-like, respectively, highlighted that these miRNAs and their target genes are likely partners in the ethylene and auxin signaling pathways. Our results indicate that these miRNAs may play an important role in regulating fruit ripening by targeting ethylene and auxin signaling pathways. Unsuitable 1-MCP treatment may disruptively repress miRNA function and cause fruit ripening disorder.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/foods10071643/s1, Figure S1: Histogram of COG classification; Table S1: RT-qPCR primers of miRNAs; Table S2: RT-qPCR primers of target genes; Table S3: Primers for plant expression vectors; Table S4: Sequencing data statistics table; Table S5: miRNA statistical results of each sample; Table S6:  miRNA expression level statistics table; Table S7: miRNAs and their targets in degradome analysis; Table S8: miRNAs and their corresponding target genes of the degradome analysis.
Author Contributions: X.Z. and X.L. designed the experiments; J.C. and Z.W. conducted the experiments; X.Z., Y.H., Y.L. and Z.S. interpreted the data; J.C. and X.Z. wrote the manuscript; Y.H., W.C., Y.L. and X.L. reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.