Cataloguing the small RNA content of honey using next generation sequencing.

Honey adulteration is a problem that effects the global honey industry and specifically, has been discovered in the Australian market. Common methods of adulteration include dilution with sugar syrup substitutes and the mislabelling of the floral and geographic origin(s) of honey. Current authentication tools rely on the molecular variability between different honeys, identifying unique chemical profiles and/or DNA signatures characteristic of a particular honey. Honey is known to contain plant miRNAs derived from its floral source. To explore the composition and variability of honey RNA molecules, this is the first study to catalogue the small RNA content of Australian polyfloral table honey and New Zealand Leptospermum scoparium honey using next generation sequencing. The data shows that in addition to miRNAs, honey contains a variety of small non-coding RNAs including tRNA-derived fragments. Moreover, the honey small RNAs are derived from a range of phylogenetic sources, including from plant, invertebrate, and prokaryotic species. The data indicates that different honeys contain unique small RNA profiles, which suggests a novel avenue in developing molecular-based honey authentication tools.

Introduction

Honey is a sweetener produced by honeybees from nectar, which they forage from flowering plants. Rising consumer demand coupled with the limited production of this natural product may explain why commercially available honey is widely adulterated, with one study finding that 27% of honeys tested were diluted with sugar syrup substitutes (Zhou, Taylor, Salouros, & Prasad, 2018). Furthermore, to gain an economic advantage, some honey products are mislabelled with regards to its floral and/or geographical origin. Identifying fraudulent honey is therefore a prevailing challenge.

The molecular composition of honey is highly variable. For example, certain honeys derived from Leptospermum plants native to Australia and New Zealand have high antimicrobial activity due to the presence of specific chemicals methylglyoxal and dihydroxyacetone which are derived from the floral source (Adams et al., 2008, Cokcetin et al., 2016).

This molecular variability allows for the development of authentication tests that can identify unique chemical signatures (Ayton, Prenzler, Raman, Warren-Smith, & Meyer, 2019), to determine the floral source and geographical location of honey production (Chin & Sowndhararajan, 2020). However, these tests can be cheated (Ayton et al., 2019). Melissopalynology (pollen analysis), is an alternative method of honey classification (Bambara, 1991). However, both these tests are unsuitable for Australian honey analysis due to the unique nature of Australian native flora (Kale Sniderman, Matley, Haberle, & Cantrill, 2018). Pollen is however a rich source of plant DNA, and studies have suggested DNA markers as useful tools for honey authentication (McDonald et al., 2018, Prosser and Hebert, 2017).

Gismondi, Di Marco, and Canini (2017) was the first to show that honey also contains RNA molecules. This study identified mRNA and rRNA molecules within monofloral European honeys and determined that they originated from the floral source. Nine plant microRNAs (miRNAs), small (22 nucleotides), conserved, non-coding RNAs that function in post-transcriptional gene regulation (Reinhart, Weinstein, Rhoades, Bartel, & Bartel, 2002), were also isolated from the honeys. Interestingly, the miRNA concentrations varied between honey types, suggesting honeys may contain unique miRNA profiles. Further, as the identified RNAs did not show signs of degradation, Gismondi et al. (2017) concluded honey RNAs are stable and could retain functionality.

To further explore honey RNA, we used small-RNA next generation sequencing (NGS). Using this approach, we determined that the honey samples we tested contain processed tRNA fragments known as tRFs (Keam & Hutvagner, 2015). tRFs have been found in many eukaryotes, including a diverse range of plant species (Gupta, Singh, Zahra, & Kumar, 2018), and in a distinct mechanism from miRNAs, regulate the translation stage of gene expression (Keam & Hutvagner, 2015). In plants, expression of specific tRFs increases in times of stress (Park & Kim, 2018).

Additionally, analyses of the sequencing data revealed that the small RNA molecules within honey originate from a variety of species and is not limited to the floral source. These data suggest that this molecular constituent of honey is more diverse that previously reported.

Materials and methods

Honey samples

Four honey samples were examined, New Zealand Manuka honey (Leptospermum scoparium, Medihoney Batch #1815; n = 2) kindly provided by Comvita (Paengaroa, New Zealand) and a polyfloral Australian table honey (Batch: PKD Q 03/05/2017 23:54; n = 2) kindly provided by Capilano (Brisbane, Australia).

Phenol equivalence assay for Antibacterial activity

Antibacterial activity of honey was determined as described in Irish, Blair, and Carter (2011). Using Staphylococcus aureus ATCC 25923 (Oxoid, Hampshire, UK) as the test organism, honey solutions were prepared fresh to a 50% (w/v) solution in sterile water and further diluted to 25% (v/v) in either sterile water to test for total activity, or a 2 mg/ml catalase solution (3187 units/mg; Sigma-Aldrich, St. Louis, USA) for non-peroxide activity. Phenol (Sigma-Aldrich, St. Louis, USA) standards of 2%−7% v/v were prepared fresh in sterile water.

Aliquots (100ul) of each honey sample and phenol solutions were dispensed into the wells of assay plates in duplicate using quasi-Latin squares to randomise sample distribution. Plates were incubated (37 °C, 18–24 h) before the diameter of each zone of inhibition was measured using digital Vernier callipers (Kinechrome, Carole Park, Australia). The mean diameter of the zone of inhibition was calculated and squared, and a standard curve (concentration of phenol versus mean squared diameter) was generated using the phenol standards. The activity (total and non-peroxide) of each honey was calculated using the standard curve, multiplied by 4.69 to account for the dilution and density of honey, and expressed as the equivalent phenol concentration.

RNA isolation from honey

Honey (10 g) was diluted 1:1 in ultrapure water and shaken in an incubator (37 °C, 180 rpm) until homogenised. Solution was centrifuged (3900 rcf, 15 mins), supernatant was discarded, and pellet resuspended in 1 ml of Trizol (Invitrogen, Carlsbad, USA). Total RNA was purified using Direct-zol RNA MiniPrep kit (Zymo Research, Irvine, USA) following the manufacturer's instructions. Total RNA concentration was measured using Nanodrop Microvolume Spectrophotometer (Thermo Scientific, Waltham, USA) and miRNA levels quantified using the Qubit miRNA Assay (Invitrogen, Carlsbad, USA). RNA quality was assessed with the small RNA kit (Agilent, Santa Clara, USA) on a 2100 Bioanalyser (Agilent, Santa Clara, USA).

Quantitative PCR

Total RNA (50 ng) from each sample was used as a template to make cDNA using the HiCapacity cDNA Kit (Applied Biosystems, Foster City, USA) following the manufacturer's instructions. qPCR reaction was prepared with TaqManTM Universal PCR Master Mix (Applied Biosystems, Foster City, USA) and run in a QuantStudio Flex 6 (Applied Biosystems, Foster City, USA) with the following Taqman miRNA assays (Applied Biosystems, Foster City, USA): (ath-mir156a-5p; ath-mir162a-3p; ath-miR396a-5p). Individual biological samples were run in technical duplicates, averaged for analysis, and standard error calculated.

Sequencing of small RNA samples

Total RNA from each sample was extracted and sent to Novogene (Hong Kong) for small RNA sequencing. Small RNA libraries were prepared using the NEBNext® Multiplex Small RNA Library Preparation Kit (Illumina, San Diego, USA) according to the manufacturers protocol and paired-end 150 bp sequencing was performed on a HiSeq 4000 (Illumina, San Diego, USA). Reads were subsequently trimmed to 50 bp single reads.

Sequencing data analysis

Raw reads were processed so that all adapter sequences were removed and reads smaller than 10 bases were excluded from further analysis. For cross-species annotation, reads were first aligned to all known precursor miRNAs (miRbase v 22) using miraligner with one mismatch allowed. All remaining unmapped reads were then aligned with bowtie (using parameters -v 1 --all --best --strata) to a set of rRNA and tRNAs spanning all three domains of life. rRNAs were downloaded from SILVA (SILVA Ref NR collection v 132) and tRNAs from GtRNAdb (v 17). NCBI’s Taxonomy database (v 2019–07-01) was used for classifying species into divisions.

NGS-based target prediction was performed with the miR-PREFeR pipeline with default settings, using processed reads from all honey samples and the Eucalyptus Grandis genome (Egrandis1_0).

Results

Plant miRNAs are present in Australian polyfloral and New Zealand Leptospermum scoparium honey

We isolated total RNA from New Zealand monofloral manuka (Leptospermum scoparium) honey, which has broad spectrum antimicrobial activity (20% phenol equivalence) and is used in honey-based wound care products (Robson, Dodd, & Thomas, 2009). For comparison, we isolated honey from an Australian polyfloral table honey with no measurable antimicrobial activity or known clinical use (<5% phenol equivalence).

We first determined whether the samples yielded sufficient quality RNA for further analysis. Initial investigations (data not shown) indicated the honeys contain only small RNA molecules (<150 nucleotides (nt)). Further analysis revealed an RNA size distribution of approximately 20–100 nt in length (Fig. 1a).

Fig. 1

RNA isolated from honey. (a) Total RNA isolated from Australian polyfloral table honey (Sample 1 and 2) and New Zealand monofloral manuka (Leptospermum scoparium) honey (samples 3 and 4) was analysed with the small RNA kit (Agilent), showing the size distribution of RNA molecules in the samples. (b) Highly conserved plant miRNAs are present in both honeys. miRNA presence was ascertained using reverse transcription-PCR (Taqman miRNA assays) (n = 2; technical replicates). Nt – nucleotide. Ct – Cycle threshold value.

To verify if the purified small RNAs have maintained their integrity, we investigated if we could detect specific plant miRNAs within the sample. We chose highly conserved plant miRNAs with roles in flower development: miR156a, miR162a, both already detected in distinct honey samples (Gismondi et al., 2017), and miR396a (Shulga et al., 2017). All three miRNAs could be amplified from the samples (Fig. 1b), indicating these miRNAs are present in the honey.

The small RNA content of honey

While plant miRNAs have been identified within honey ((Gismondi et al., 2017), Fig. 1b), we catalogued the small RNA content of honey using next generation sequencing.

The bioinformatic pipeline used to analyse the small RNA sequencing data is illustrated in Fig. 2a. To identify small RNAs across different species, we aligned our reads to a comprehensive set of miRNAs, ribosomal RNAs (rRNAs) and transfer RNAs (tRNAs). We aimed to i) identify the different types of RNA present and ii) determine the phylogenetic classification of the RNA molecules within the honey.

Fig. 2

Small RNA sequencing of honey samples. (a) Bioinformatic approach. (b) length distribution (nucleotide) and RNA type (miRNA, ribosomal RNA (rRNA), transfer RNA (tRNA)) of the sequenced reads from polyfloral and monofloral honeys.

Our analyses revealed that the majority of RNA present in the polyfloral (68.9–76.1%) and monofloral (60.6–66.5%) honey samples did not align to either predicted or validated small RNA sequences and therefore their identity is unknown. One reason for this is the lack of fully annotated genomes available for Australian native species, including the plants commonly used as nectar sources for honey production (e.g. Corymbia species). The Leptospermum scoparium genome has been sequenced (Thrimawithana, Jones, Hilario, Grierson, & Ngo, 2019), however due to the cultural importance of this plant to the indigenous Māori of Aotearoa people of New Zealand, the genome sequence has not been published. Furthermore, honey RNA molecules may undergo editing and/or modifications during honey production and processing, which may influence their ability to be accurately identified. Despite these caveats, we confirmed that the plant miRNAs identified by RT-PCR (Table 1) are present, and successfully identified many others (Fig. 2b). We also identified other kinds of small RNA molecules. Interestingly, most were derived from tRNA molecules (22.9% and 27.6% for polyfloral and monofloral respectively), and a significant fraction from rRNAs (3.55% and 7.53% for polyfloral and monofloral respectively; Fig. 2b). miRNAs were the smallest identified fraction at 1.03% for polyfloral and 1.30% for monofloral honey. While it is yet to be discovered if these RNAs have any physiological function within honey, their roles in gene expression regulation and their highly conserved nature (Engels & Hutvagner, 2006), raises further questions about their possible role in the bioactive properties of honey.

Table 1

Highly conserved miRNAs identified in honey.

Polyfloral Australian table honey		Monofloral New Zealand manuka honey
miRNA	Total Read Count	miRNA	Total Read Count
miR156a	64.38	miR156a	130.34
miR162a	439.91	miR162a	241.40
miR396	535.59	miR396	650.88

Values represent counts per million (CPM) raw reads normalised to all reads which map to at least one plant sequence.

Finally, to identify plant-specific miRNAs that have not been previously annotated, we employed a miRNA prediction tool. As we are currently unable to access the Leptospermum scoparium genome, we used the phylogenetically related species, Eucalyptus grandis. On average, 58.01% of reads mapped to the E. grandis genome from the polyfloral and 48.17% from the monofloral honey, suggesting a significant portion of the unmapped reads from the cross-species analysis are either unannotated small RNAs or degraded RNAs from plants. This result was expected as eucalypt plants are a common nectar source for Australian honey. In total, this analysis identified 52 predicted plant miRNAs and was generally consistent with the miRNAs identified in our cross-species analysis (Supp Table 1).

Honeys have unique small RNA expression profiles

The small RNA composition between the two honey samples show clear differences (Fig. 2b). The polyfloral honey has a distinct tRF peak, 13 nucleotides in length. Conversely, the monofloral honey has a distinctive tRF peak, 35 nucleotides in length. A larger data set is required to determine if the RNA profile differences observed here are truly representative of the honeys tested. However, this disparity confirms that RNA composition of honey is not uniform, and therefore has the potential to distinguish between honey types.

tRFs have been classified into two broad groups according to their biogenesis. Half-tRNAs (30–35 nts) are generated by a single cleavage event in the anticodon loop of a mature tRNA. half-tRNA expression correlates with perceived stress, and function to repress translation in mammalian cells (Yamasaki, Ivanov, Hu, & Anderson, 2009). To determine if the 35-nucleotide tRFs observed in the honey samples are half-tRNAs, we mapped their read distributions across the full length tRNA gene (Fig. 3a). The data indicates that these tRFs are a product of a specific and controlled processing event, producing a stable half-tRNA.

Fig. 3

tRNA fragments in honey are formed from regulated processing of full length tRNA. (a) Example of a read distribution sequenced read across a full length tRNA gene which indicates a controlled endonucleolytic event forming a stable 35-nt tRNA fragment. (b) Example of a read distribution sequenced read across a full length tRNA gene which indicates controlled endonucleolytic event forming a stable 3ʹ tRF.

Precursor and mature tRNAs are processed via several mechanisms to form distinct tRF molecules. 5ʹ tRFs are formed by the endonucleolytic cleavage of the TψC-arm of a tRNA, while 3ʹ CCA tRFs are generated by cleaving the D-arm. These fragments typically range from 13 to 22 nucleotides long, and some function to regulate gene expression in human cells (Sobala & Hutvagner, 2013). From our read distribution data, the tRFs present in the polyfloral honey appear to be processed from the 3ʹ end of a full length tRNA (Fig. 3b).

RNA in honey comes from a variety of sources

The only documented origin of RNA found in honey is from the floral source (Gismondi et al., 2017). Our data confirms this as 21.1% and 17.0% of annotated polyfloral and monofloral small RNAs mapped exclusively to plants. However, when we performed species annotation of the sequenced reads, we determined that plants are not the only source of honey RNA (Fig. 4). We identified that many reads mapped predominately to species external to the plant kingdom, including invertebrate, bacterial and viral species. This includes miRNA sequences that map to invertebrate genomes such as mir-14, a miRNA conserved across invertebrates which has an important role in insect development (Mead & Tu, 2008) and mir-276, a critical regulator of the reproductive cycle (Lampe, Jentzsch, Kierszniowska, & Levashina, 2019). A miRNA specific to the genomes of eusocial Hymenoptera insects, which includes Apis mellifera (the European honeybee), is also present in the honey. European honeybees underpin the Australian and New Zealand honey industry, therefore identifying these reads gives credence to our method of species assignation of the RNA molecules within the honey.

Fig. 4

RNA in honey derives from a wide variety of species. Species annotation of the mapped and unmapped reads from the sequencing data shows that the RNA aligns with a wide variety of species from different phylogenetic domains. Counts were normalised and averaged across biological replicates.

Due to the conservation of sequences across species, this analysis cannot conclusively assign an RNA molecule to a particular species. However, this data does suggest the small RNA in honey originates from a range of species that span across different evolutionary kingdoms, which may provide further insight into honey composition and comprise a future avenue for honey authentication.

Discussion

The comprehensive approach of next generation sequencing revealed that honey contains several RNA types which derive from a range of different species. The RNA preserved within the New Zealand Leptospermum scoparium honey and the polyfloral Australian honey used in this study were small molecules of less than 150 nts in length. The absence of larger RNA molecules differs from Gismondi et al. (2017), who isolated an rRNA molecule that was 844 nts long from European honey. This discrepancy may be due to the distinct floral, geographic, and environmental conditions of honey production. Varied beekeeping practices, honey processing procedures and the conditions and duration of honey storage prior to RNA isolation may also affect the honey components. Furthermore, the differing chemical profiles of honey may influence its suitability in preserving RNA molecules of different lengths.

We identified small regulatory RNAs, including miRNAs and tRFs, within honey. Further research is required to determine if these molecules retain any physiological function or if they contribute to the nutritional properties of honey. It is possible that they have a role in attracting and manipulating mutualistic species and protecting the plant from species which would exploit the nectar resource. These qualities are an inherent aspect of nectar in plant-insect relationships (Nepi, Grasso, & Mancuso, 2018).

The successful amplification of plant miRNAs (Fig. 1b) indicates that the floral source(s) of the honey likely contributes a proportion of the RNA content. However due to the highly conserved nature of the RNAs identified, species-specific markers are required to confirm the plant source(s). Honey small RNAs derive from several phylogenetic groups including invertebrates and prokaryotes. The discovery of miRNAs which map to the honeybee genome (and other invertebrates) is a novel finding. These miRNAs may be suitable markers capable of validating unadulterated honey.

In plants, miRNAs are highly mobile molecules, capable of moving the length of the plant as they can enter and travel through the phloem (De Felippes, Ott, & Weigel, 2011). As plant miRNAs can survive and function in a variety of tissues it is perhaps not unsurprising that they are present in sugar-rich nectar and consequently in honey.

miRNAs are signalling molecules which form a highly intricate regulatory network essential to the development, growth and survival of the plant, and function to coordinate plant-wide responses (Skopelitis et al., 2018). Most of the plant-derived miRNAs that we identified from honey are highly conserved with understood functions. We investigated the known roles of the 10 most abundant plant-derived miRNAs identified in the honey samples. While the known functions of these miRNAs have not been confirmed in Leptospermum or Australian native species, reported regulatory functions of these miRNAs in other plant species can be broadly categorised into three main groups: tissue development; immunity, disease resistance and defence mechanisms; and stress responses, including drought (Table 2). Therefore, these miRNAs may have potential as markers which indicate the developmental stage, specific environmental stresses or health status of the plant at the time nectar was harvested for honey production. This insight into the health status of native Australian flora used as a nectar source may provide valuable information which could help to maximise Australian beekeeping practises and the maintenance of native Australian flora. Further research is required to determine if these factors influence the quality of honey produced.

Table 2

Functions of the 10 most abundant, identified, plant-derived miRNAs found in honey.

Function	miRNA(s)	References
Regulates plant tissue development	miR-166; miR-156; miR-482; miR-894; miR-396	Bai et al., 2017, Li et al., 2010, Liu et al., 2018, Shulga et al., 2017
Regulates plant immunity/defence mechanisms	miR-396; miR-482; miR-2089	Soto-Suárez et al., 2017, Zhai et al., 2011
Regulates plant stress response	miR-172c; miR-1511; miR-6300; miR-8175	Anna et al., 2019, Ferdous et al., 2015, Htwe et al., 2015;Liang, Liu, Hao, Li, & Luo, 2019

miRNAs are critical regulators of gene expression in both plants and animals. However, there are significant differences in the biogenesis and function of miRNAs between the two evolutionary kingdoms (Axtell, Westholm, & Lai, 2011). Plant miRNAs are methylated at their 3ʹ ends, whereas most animal miRNAs do not have this modification (Yu et al., 2005). It has been suggested that miRNAs have bioactivity potential i.e. being absorbed into blood from consumed food, entering the miRNA regulatory pathway and regulating gene expression. However, subsequent studies failed to confirm ingested miRNAs have regulatory function (Dickinson et al., 2013), therefore, there is no evidence that honey miRNAs contribute to its bioactive properties once consumed.

While many of the RNA molecules within the honey samples remain unmapped, our data indicates that honeys contain distinct RNA signatures, composed of a variety of RNA types and derive from a range of species. Expression of regulatory small RNAs is dynamic and changes in response to specific environments and stresses. Therefore, a comprehensive analysis of honey RNA signatures may provide further insight into the manufacture of honey, for example the geographical location or climate conditions of the floral source(s) when the nectar was collected.

Declaration of Competing Interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Comvita New Zealand and Capilano Honey Ltd. provided materials in the form of honey samples for the work described in the manuscript. 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.