extfieldsofstudy: Medicine Chemistry

text:
Significance of Squalene in Rice Bran Oil and Perspectives on Squalene Oxidation.

As an intermediate metabolite during the biosynthesis of sterols, squalene is found ubiquitously in plants and animals. In rice, squalene is contained in rice bran, and consequently, squalene in rice bran oil has gained attention. Studies have shown that the intake of squalene from food sources demonstrate various physiological benefits such as the prevention of cancer and cardiovascular disease. Squalene is also known as an effective antioxidant in edible oils. However, due to its chemical structure, squalene is susceptible to oxidation, which may cause a decline in the nutraceutical and antioxidative effects of squalene in edible oils. Oxidative degradation of squalene also results in the formation of scission products (i.e., aldehydes and ketones) which may lead to off-flavor. Since the rate of squalene oxidation depends on the factors that induce its oxidation (i.e., light or heat), emphasis on oxidation mechanisms is necessary. It has been demonstrated in previous studies that the oxidation products formed by the singlet oxygen oxidation and free radical oxidation of squalene are different, and more recently, we demonstrated that different squalene monohydroperoxide isomers are formed by each oxidation mechanism. We herein discuss the significance of squalene in rice bran oil as well as the oxidative degradation of squalene in edible oils with focus on oxidation mechanisms.

Nutraceutical and Antioxidative Properties of Squalene in Rice Bran Oil
Squalene is a C30H50 hydrocarbon that possesses a characteristic chemical structure consisting of six isoprene units (Fig. 1A). Squalene is known to be an intermediate metabolite during the biosynthesis of sterols and is ubiquitously found in plants and animals. Commercially, squalene is used for various purposes including vaccine adjuvants, dietary supplements, and cosmetic products, and is obtained mainly from shark liver oil. However, due to recent concerns in marine preservation, squalene from alternative sources, especially from plant origin, has attracted interest (1). In rice, squalene is contained in rice bran as an unsaponifiable constituent. Thus, studies have investigated the presence of squalene in rice bran oil, a product derived from rice bran. The squalene content in rice bran oil and its derivatives have been quantified with the use of different methods such as thin-layer chromatography (TLC) (2)(3)(4), gas chromatography with flame ionization detec-tion (GC-FID) (5-8), high-performance liquid chromatography with ultraviolet detection (HPLC-UV) (9), and high-performance liquid chromatography with diodearray detection (HPLC-DAD) (10) ( Table 1).
The significance of squalene in foods lies in the fact that it demonstrates various physiological effects including prevention of diseases such as cancer. For example, administration of a diet containing 1% squalene to rats was found to suppress the formation of aberrant crypt foci, an early precursor of colon cancer, suggesting the preventive effect of squalene on colon cancer (11). Dietary squalene was also shown to inhibit 4-(methylnitrosamino)-1-(3-pyridyl)-1-butanoneinduced lung tumorigenesis in rats (12). Other studies have also demonstrated the putative preventive effects of squalene on skin cancer (13) and certain types of breast cancer (14). Additionally, a diet supplemented with 2% squalene was found to attenuate myocardial infraction in rats (15) and intake of 1 g/kg squalene per day demonstrated a decrease in atherosclerotic lesion size in male ApoE knockout mice (16), which both illustrate the potential cardioprotective effect of squalene. Importantly, although squalene is an intermediate in the biosynthesis of cholesterol, animal studies as well as clinical trials have demonstrated that dietary intake of squalene does not affect serum cholesterol levels (17,18).
In addition to such properties as a food factor, the antioxidative property of squalene has also been thoroughly investigated. In a study that evaluated the effect of squalene on the auto and photooxidation of polyunsaturated fatty acids, the role of squalene as a peroxyl radical scavenger was suggested (19). Other studies have also demonstrated that squalene acts as an effective quencher of singlet oxygen (20) and that squalene demonstrates synergistic antioxidative properties in combination with a-tocopherol (21). A similar effect was observed during the photooxidation of olive oil, where it was suggested that the presence of squalene prevents the degradation of a-tocopherol via regeneration of a-tocopherol radicals (22). Additionally, a concentration-dependent antioxidative activity of squalene was observed during the heat-induced oxidation of rapeseed oil (23) and olive oil (24). Thus, the above studies suggesting the antioxidative effects of squalene imply, at least in part, that squalene may contribute to the prevention of oxidation in edible oils.

Oxidative Degradation of Squalene
Although squalene possesses the benefits described above, because of its characteristic chemical structure containing six double bonds, squalene is known for its susceptibility to oxidation. Since oxidation of squalene leads to the loss of its nutraceutical and antioxidative potential, it is essential to know the mechanisms by which squalene is oxidized and degraded, as well as the consequences of such oxidative degradation of squalene in foods.
Oxidative degradation of squalene in edible oils has been investigated in several studies. In a study that evaluated the long-term storage of olive oil under room temperature and diffused lighting, squalene content was found to significantly decrease, especially under conditions where the oil was exposed to air (i.e., nearly half of the squalene was lost after 6 mo of storage) (25). However, a contrasting result was obtained when the long-term storage of olive oil in dark was evaluated; loss of squalene was limited to 10% even after 24 mo of storage (26). Likewise, when the auto and photooxidation of pumpkin seed oil was investigated, while both conditions resulted in a considerable decrease in squalene content, loss of squalene was greater during photooxidation than that during autooxidation (27). Such results suggest that the rate of squalene degradation depends on the cause of oxidation (i.e., light or heat) which will be discussed further in the following section.
Another potential consequence of squalene oxidation besides the loss of its positive properties is the negative effect of secondary oxidation products. Oxidation of squalene at 60˚C for 24 h resulted in the formation of volatile organic hydrocarbons including methane, propene, isobutene, and isoprene (28). In a more recent study, heating of squalene at 120˚C for a shorter time period identified abundant scission products of squalene, namely ketones and aldehydes (29). Although the structures of these ketones and aldehydes were not identified, we assume that these may include compounds that share a similar structure with hexanal and nonanal, which are well-known volatile components of the off-flavor of edible oils (Fig. 2). Interestingly, in the same study, such secondary oxidation products were found to demonstrate a prooxidative effect on olive oil, suggesting the further complications caused by squalene oxidation. Additional studies have shown that the photooxidation of squalene also forms formaldehyde and malonaldehyde via scission products including 6-methyl-5-hepten-2one (30,31). Although such studies regarding the secondary oxidation products of squalene are rather limited, the above studies demonstrate the further need for investigation into ways to prevent squalene oxidation.

Evaluation of Squalene Oxidation Mechanisms: LC-MS/MS Analysis of Squalene Monohydroperoxide (SQOOH) Isomers
As described above, evaluation of squalene oxidation with focus on the factors that induce its oxidation (i.e., oxidation mechanisms) is crucial. Oxidation mechanisms, in a more general sense, can be considered as the type of reactive oxygen species that induce lipid oxidation. For example, light-induced oxidation (i.e., photooxidation) is typically caused by singlet oxygen, which directly reacts with double bonds via the ene reaction to form lipid hydroperoxides as the primary oxidation product (32,33). Such reaction of singlet oxygen with squalene is known to be rather reactive (20), and previous studies have demonstrated that squalene monohydroperoxide (SQOOH) was mainly formed (34,35). On the other hand, heat-induced oxidation (i.e., autooxidation) of lipids is generally caused by free radicals such as the hydroxyl radical and alkoxyl radical. Free radical oxidation is known to initiate from the abstraction of a hydrogen atom to form a lipid radical, and the subsequent addition of molecular oxygen results in the formation of lipid oxidation products (32,36). With regard to squalene, studies on the free radical oxidation of squalene has led to contrasting views. Whereas early studies identified cyclic dihydroperoxides as the main oxidation product (37), more recent studies have reported that alcohols and epoxides are mainly formed (29,38).
Meanwhile, we have been conducting studies on SQOOH, the primary oxidation product of squalene oxidation, with special focus on its isomers. Squalene, when oxidized, is known to form six SQOOH isomers (Fig. 1B), but the analysis of SQOOH at the isomeric level had not been achieved previously. We were first able to achieve SQOOH isomer analysis by using LC-MS/MS with atmospheric pressure chemical ionization (APCI), but due to the nature of APCI fragmentation, fragment ions with similar patterns were produced, making it difficult to selectively analyze SQOOH isomers (35). Hence, we recently developed an improved method that enabled for the highly selective detection of six SQOOH isomers (39), based on evidence that the use of sodium ions during LC-MS/MS analysis enables the selective analysis of other lipid hydroperoxide isomers (40)(41)(42)(43). With the use of this method, it was identified that squalene oxidation by singlet oxygen oxidation and free radical oxidation yield different SQOOH isomers; singlet oxygen oxidation formed six SQOOH isomers in nearly equal amounts whereas free radical oxidation mainly formed 2-OOH-SQ and 3-OOH-SQ (44). Moreover, when two antioxidants, a-tocopherol and b-carotene, were evaluated  for their antioxidative effect to prevent squalene oxidation, each antioxidant was shown to be effective mainly towards only one of either singlet oxygen oxidation or free radical oxidation. Such recent results suggest that identifying oxidation mechanisms by analyzing SQOOH isomers, and then applying an antioxidant that is effective towards the identified oxidation mechanism is an effective way to prevent squalene oxidation (Fig. 3).

Conclusions and Future Perspectives
Squalene contained in foods including rice bran oil, possesses positive characteristics including its nutraceutical benefits and antioxidative properties. However, squalene is known to be susceptible to oxidation, and oxidative degradation of squalene may cause a decline in such properties. Studies have suggested that the rate of squalene degradation depends on the cause of oxidation, and additionally, secondary oxidation products of squalene have been suggested to negatively affect the properties of various edible oils. As such, it is necessary to evaluate the oxidation of squalene in food products with focus on oxidation mechanisms. With regard to oxidation mechanisms, previous studies have demonstrated that the photooxidation (i.e., singlet oxygen oxidation) of squalene results in the formation of mainly SQOOH, whereas heat-induced oxidation (i.e., free radical oxidation) of squalene has led to contrasting views. Such information regarding the oxidation products derived from the free radical oxidation of squalene is yet to be clarified. Moreover, studies on the detection of squalene oxidation products from actual edible oil samples including rice bran oil is highly anticipated. Finally, our recent results concerning the analysis of SQOOH isomers for determination of squalene oxidation mechanisms may be valuable for prevention of squalene oxidation, and its application to edible oils are targets for future studies.

Disclosure of State of COI
The authors declare no conflicts of interest.


extfieldsofstudy: Medicine Biology

provenance: AgriculturalAndFoodSciences,2019-2020-Agricultural And Food Sciences-2019.jsonl.gz.gz

text:
Single-molecule real-time transcript sequencing identified flowering regulatory genes in Crocus sativus

Background Saffron crocus (Crocus sativus) is a valuable spice with medicinal uses in gynaecopathia and nervous system diseases. Identify flowering regulatory genes plays a vital role in increasing flower numbers, thereby resulting in high saffron yield. Results Two full length transcriptome gene sets of flowering and non-flowering saffron crocus were established separately using the single-molecule real-time (SMRT) sequencing method. A total of sixteen SMRT cells generated 22.85 GB data and 75,351 full-length saffron crocus unigenes on the PacBio RS II panel and further obtained 79,028 SSRs, 72,603 lncRNAs and 25,400 alternative splicing (AS) events. Using an Illumina RNA-seq platform, an additional fifteen corms with different flower numbers were sequenced. Many differential expression unigenes (DEGs) were screened separately between flowering and matched non-flowering top buds with cold treatment (1677), flowering top buds of 20 g corms and non-flowering top buds of 6 g corms (1086), and flowering and matched non-flowering lateral buds (267). A total of 62 putative flower-related genes that played important roles in vernalization (VRNs), gibberellins (G3OX, G2OX), photoperiod (PHYB, TEM1, PIF4), autonomous (FCA) and age (SPLs) pathways were identified and a schematic representation of the flowering gene regulatory network in saffron crocus was reported for the first time. After validation by real-time qPCR in 30 samples, two novel genes, PB.20221.2 (p = 0.004, r = 0.52) and PB.38952.1 (p = 0.023, r = 0.41), showed significantly higher expression levels in flowering plants. Tissue distribution showed specifically high expression in flower organs and time course expression analysis suggested that the transcripts increasingly accumulated during the flower development period. Conclusions Full-length transcriptomes of flowering and non-flowering saffron crocus were obtained using a combined NGS short-read and SMRT long-read sequencing approach. This report is the first to describe the flowering gene regulatory network of saffron crocus and establishes a reference full-length transcriptome for future studies on saffron crocus and other Iridaceae plants.

Background
Crocus sativus L, commonly known as saffron crocus, is prized for purple flowers that are well known for producing spice saffron from the filaments. Spice saffron is the most valuable spice used as a fabric dye and in traditional medicine with special medicinal effects of promoting blood circulation, cooling blood and detoxifying, thereby relieving depression and soothing nerves [1]. As a valuable traditional Chinese medicine, saffron is widely used in China and Europe. Saffron crocus blooms only once a year and unlike most spring-blooming plants, saffron crocus does not blossom until autumn. In China, the daughter corms began to grow at the end of January and matured at the end of May and subsequently, entered a dormant period until mid-August. During the period, the corms were dug out from the soil when the leaves turned yellow and wilted and moved into the door to store. Experiencing the high temperature treatment in summer (ranged from 23 to 27°C), the buds were broken up from dormancy in the middle of August and the floral primordia began to initiate. When the average room temperature fell to 15-17°C in mid-autumn, most apical buds were in blossom [2]. Basically, the corms had 1-3 apical buds and 6-10 lateral buds depending on their weight. Each apical bud germinated 1-3 flower primordia while lateral bud usually did not blossom. Occasionally, one or more lateral buds of corms weighing more than 30 g also blossomed. The corms weighing less than 6 g cannot blossom. As soon as all the flowers were picked up indoors, the corms were planted in the soils until the new daughter corms matured in the next May. Planting and harvesting corms as well as collecting red stigmas from flowers, is performed manually. To produce 1 kg of dry saffron, 110,000-170,000 flowers are harvested and 40 h of labour are needed to pick 150,000 flowers. Such labour-intensive cultivation practices make saffron a high expensive crop with prices ranging from $500 to $5000 per pound at wholesale and retail rates [3]. Due to limited natural resources for saffron crocus plants, inefficient cultivation, and low yield, saffron is becoming even moreexpensive and is well known as "red gold" [3]. It is highly important to explore comprehensive genetic information for breeding and improving its biological traits.
Increasing the flower number of saffron crocus is a viable way to produce more saffron to meet the everincreasing demand in the market [4,5]. Research has been conducted to investigate the factors that affect floral development including temperature, photoperiod, corm size, and bud position [2,6]. We can obtain samples of different flowering quantities by controlling these factors artificially. Therefore, C. sativa is a good material for studying the development of flowering. Many genes related to plant floral development have been discovered along with the rapid development of technology in molecular biology. For example, long-day conditions can promote Arabidopsis flowering through the function of FLOWERING LOCUS T (FT) protein, which is considered to be the main component of "florigen" [7,8]. The transcription factor Flowering Locus C (FLC) is a key regulator of the vernalization process of Arabidopsis thaliana. The transcription factor PIF4 is a major regulator of high temperature-induced flowering [9]. Using the FT gene in Arabidopsis thaliana as a reference, Tsaftaris et al. cloned a CENTRORADIALIS/TERMINAL FLOWER1 (CsatCEN/TFL1) like gene [10] and three FT-like genes [11] from the flowers, flower buds, leaves, and corms of saffron crocus, respectively, and further proved that their expression patterns were tissue-specific and depended on the flower developmental stage. Other studies found a serial potential flower-related gene in saffron crocus, for instance, B-class paleo AP3-like genes (CsatAP3-like) [12], AP1-like MADS-box genes [13], Bclass floral homeotic genes PISTILLATA/GLOBOSA [14], E-class SEPALLATA3-like MADS-box genes [15], and CsMYB1, a transcription factor belonging to the R2R3 family [16]. Later, NGS-based RNA-seq technology was widely used for gene discovery, which led to the identification and functional characterization of flowering genes in various species. For example, trehalose 6phosphate synthase and squamosa promoter-binding protein-like genes promoted the floral induction of apple trees [17]. A series of genes related to the circadian clock are important key regulators for the flower development of hibiscus [18]. Using NGS-based RNA-seq technology, Baba et al. [19] and Jain et al. [20] discovered the gene expression of saffron crocus involved in apocarotenoid biosynthesis and further explored the expression profiling of zinc-finger transcription factors [21]. However, the underlying molecular mechanism controlling and/or affecting the number of flowers of saffron crocus has not been determined. The genome has not been fully elucidated to date, even in the whole Iridaceae family, only de novo assembly based short-fragment transcriptome of saffron crocus was provided by Illumina RNA-seq sequencing [19][20][21].
Compared with the NGS platform, PacBio Iso-Seq can obtain a collection of high quality full-length transcripts without assembly, which is especially important for species without reference genome sequences. Some transcripts might contain repeat regions, whereas transcripts of different gene isoforms show high sequence similarity. The assemblies of short sequencing reads often encounter complications without reference genome sequences. The problem seems more severe for saffron crocus, because of its relatively larger genome size [29] (greater than 10 Gb) and polyploid characteristics [30] (2n = 3x = 24). Saffron crocus consists mainly of repetitive DNA sequences, such as retrotransposon and satellite DNA [31], resulting in particular challenges for the accuracy of short-read assembly. The PacBio Iso-Seq technology can overcome these difficulties by generating sequence information for the full length sequence as a single sequence read without further assembly.
In this paper, NGS and SMRT sequencing were combined to generate two sets of full-length transcriptomes of flowering and non-flowering saffron crocus. Moreover, differentially expressed full-length transcripts of flowering and non-flowering saffron crocus were identified and characterized.

Plant materials
Saffron crocus plants were cultivated at a research farm at South Tai Lake Agricultural Park, Huzhou (longitude 120.6°E, latitude: 30.52°N, elevation 0 m), using a twostage cultivation method: corms planted in soil to allow them to grow outdoors and be cultivated indoors without soil [32]. In May 2016, dormant corms were excavated from the field and stored indoors for approximately half a year until flowering.
Two sample pools were set up to establish the PacBio Iso-seq libraries of flowering saffron crocus and nonflowering saffron crocus separately. One sample pool was constructed for the full-length transcript set of flowering saffron crocus, which included 1) top bud tissues, 2) tuber tissues of flowering corms (5-7 mm, ≈20 g) (recently differentiated flower primordia and leaf primordia), 3) pistils, 4) stamens of flowering corms (≈20 g) when colours turned from yellow to red, and 5) leaves of flowering corms (≈20 g) when colours turned from white to green), and 6) purple petals of flowering corms (≈20 g). The other sample pool was constructed for the fulllength transcript set of non-flowering saffron crocus, which included 1) top bud tissues, 2) lateral bud tissues, 3) tuber tissues of non-flowering corms (5-7 mm, ≈20 g), 4) leaves of non-flowering corms (≈20 g) when turned from white to green, and 5) top bud tissues of nonflowering corms (5-7 mm, ≈6 g) (Additional file 1: Figure S1).
Meanwhile, an additional five groups of saffron crocus corms were prepared to construct higher-accuracy short-read libraries using an Illumina RNA-seq method. The sample pools included 1) top buds of flowering saffron crocus corms, 2) paired top buds of non-flowering saffron crocus corms (≈20 g) that were split into two parts and cultivated at room temperature (20-25°C, flowering phenotype) and 10°C (non-flowering phenotype) for 15 days, 3) lateral buds of flowering saffron crocus corms, 4) paired lateral buds of non-flowering saffron crocus corms (≈30 g), and 5) top buds of nonflowering saffron crocus corms (≈6 g) (Additional file 1: Figure S1). All five bud samples were collected when they were 5-7 mm long. A total of 15 plants, (three plants per group) were harvested to construct 15 Illumina RNA-seq libraries.
All the samples prepared for both PacBio Iso-seq and Illumina RNA-seq sequencing were immediately frozen in liquid nitrogen until RNA was isolated.

RNA preparation
All tissues were ground in liquid nitrogen and total RNA was extracted using an RNeasy@Plant Mini Kit (Qiagen Corporation, Hilden, Germany) according to the manufacturer's protocol. The isolated RNA samples were detected using 1% agarose electrophoresis to avoid degradation and genomic DNA contamination.

PacBio Iso-Seq library preparation and sequencing
PacBio Iso-Seq libraries of flowering and non-flowering saffron crocus were constructed separately. After RNA samples were tested, total RNAs from each set of sample pools (flowering/non-flowering saffron crocus) were mixed and isolated for Poly (A) RNA using a Poly (A) Purist™ MAG Kit (Invitrogen, Carlsbad, CA, USA). Poly (A) RNA was reverse transcribed into cDNA using a SMARTer® PCR cDNA Synthesis Kit (Clontech, Mountain View, CA, USA) with SMARTScribe® MMLV RT enzyme (Takara, Dalian, China). The cDNA products were further amplified with the optimal number of cycles using KAPA HiFi PCR Kits. The PCR products were screened using a BluePippin® Size Selection System (Sage Science, Beverly, MA, USA), and three fractions containing fragments of 1-2, 2-3, and > 3 kb in length were obtained. The sorted fragments of PCR products were amplified again using KAPA HiFi PCR Kits to produce enough DNAs for constructing sequencing libraries. The PCR products were subjected to construct SMRTBell libraries using a SMRTBell Template Prep Kit (Pacific Biosciences, Menlo Park, CA, USA) after fragment ends were repaired and the blunt hairpin adapters at both ends of the DNA fragments were connected. A total of 16 SMRT cells, that is, eight SMRT cells (3 cells for the 1-2 kb library, 3 cells for the 2-3 kb library and 2 cells for the > 3 kb library) run for each sample pool, were analysed using a PacBio RS II platform (Pacific Biosciences, Menlo Park, CA, USA). Figure 1a lists the workflow for the whole Pac-Bio Iso-seq data processing.
Illumina RNA-seq library preparation, sequencing, and Contigs assembly Fifteen RNA samples from saffron crocus buds were used for Illumina RNA-seq library construction and sequencing. Total RNA was enriched using Oligo (dT) magnetic beads and randomly broken into short fragments that were further used as a template to synthesize cDNA with random hexamer-primers. The cDNA products were endrepaired, A-tailed, and added with Illumina paired-end adapters. The fragments were selected using AMPure XP beads and PCR amplified to obtain sequencing libraries that were qualified and paired-end sequenced with an Illumina Hiseq 2000 (Illumina, San Diego, CA, USA).
The raw reads of the sequences were obtained by removing adapter reads, reads with length of < 100 bp, and reads with content of ambiguous bases 'N' > 5%. De novo assembly of transcriptome sequencing without reference genome, including steps of Inchworm, Chrysalis, and Butterfly with default parameters was conducted using Trinity software.
Quality control, error correction of PacBio reads and Contigs mapping between corrected PacBio reads and Contigs from RNA-seq The raw data from the PacBio RS II platform were filtered using SMRTLink software (version 4.0) to obtain Post-Filter Polymerase reads, namely, CCS, when the adaptors, subreads < 50 bp, polymerase reads < 50 bp and accuracy of polymerase reads < 0.75 were deleted. CCS were further self-corrected and filtered with the criterion of full passes > 1 and the predicted consensus accuracy > 0.8 toobtain high-quality reads of inserts (ROIs). ROIs were classified into non-full-length reads and full-length reads (including full-length non-chimeric reads and full-length chimeric reads) based on the presence and location of 3′ primer, 5′ primer and polyA. Full-length non-chimeric reads were corrected by the CEC algorithm and produced Unpolished Consensus Sequences (UCS). The UCS and the remaining ROIs were further corrected using Quiver software to obtain polished high-quality isoforms (accuracy > 0.99) and polished low-quality isoforms.
Subsequently, all Quiver-polished isoforms were mapped to Trinity-assembled contigs from RNA-seq to produce Trinity-corrected Pacbio Isoforms using LoR-DEC software [33]. By aligning the Trinity-corrected Pacbio Isoforms to contigs assembled by Trinity with a high level of similarity (> 99% threshold), the longest contigs were assigned to the duplication-removed and corrected long reads (DRCLR). The DRCLR was corrected to remove redundant information using CD-HIT software (version 4.6) and regarded as Unigene.

Unigene annotation
To predict unigene function, unigenes were searched against five databases, including Cluster of Orthologous Groups of proteins (COG), SwissProt, NCBI nonredundant (NR), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Functional annotation of unigenes was obtained from sequence similarity alignment using the BLAST algorithm with a criterion of E-value <1e-10.

Prediction of coding DNA sequence and protein
All the isoforms were used to predict the coding sequences (CDS) and protein sequences using ANGEL software with Arabidopsis thaliana and Phalaenopsis equestris (orchid) genomes as the reference genomes. The genome of the Iridaceae family has not been fully elucidated to date [19]. Among all species with known genomes released recently, Phalaenopsis equestris has the most homology with saffron crocus [34].

SSR annotation and long non-coding RNA identification
SSRs (simple sequence repeats) were searched using MISA software (version 1.0) [35]. Long non-coding RNAs (lncRNA) were predicted according to the guiding principles of lncRNAs pipeline (https://bitbucket.org/ arrigonialberto/lncrnas-pipeline) with PLEK (an improved k-mer scheme tool) as the core algorithm [36]. PLEK is widely used to discriminate protein-coding mRNAs and non-coding RNAs and has the ability of predict all possible open reading frames (ORFs) and translate the sequences into peptides.

Alternative splicing analysis and validation
The alternative splicing (AS) events were predicted based on the BLAST alignment of DRCLR to the Trinity-assembled contigs from RNA-seq sequencing using default parameters. AS events were defined when the alignment gaps were longer than 50 bp and were at least 100 bp from the 3′ and 5′ ends [33]. The specific AS presented in only the PacBio Iso-seq library of flowering or non-flowering saffron crocus were screened separately. To validate the accuracy of the AS detected with PacBio Sequencing, RT-PCR of three randomly selected unigenes, PB.174, PB.313 and PB.988,was performed. Total RNA of saffron crocus buds was extracted as described above. The PrimeScript II 1st Strand cDNA Synthesis Kit (TaKaRa, Japan) and SYBR Premix Ex Taq II (TaKaRa, Japan) were used for reverse transcription reaction and PCR assay. Specific primers (Additional file 2: Table S1) of the chosen genes were designed using Primer Premier 5.0 software (Premier, Vancouver, British Columbia, Canada) according to the homologous sequences at the upstream and downstream ends of all the different alternative splicing fragments. The PCR amplification procedure included 98°C 10 s, 56°C 30 s, 72°C 3 min for 30 cycles and then 72°C extended for 5 mins. PCR products were monitored by 1% agarose gel electrophoresis. Sequencing of the PCR products further confirmed the correctness of the amplification.

Screening differentially expressed Unigenes and GO and KEGG enrichment analyses
The expression levels of all the unigenes in fifteen samples were assayed based on the Illumina short reads dataset, and reference sequences were the unigene libraries. Relative gene expression levels of each unigene were determined by FPKM (fragments per kilobase of transcript per million mapped reads) and differentially expressed unigenes were screened using DESeq2 R with parameter cutoff p-value < 0.05, FDR < 0.01 and fold change ratio > 2.
Differentially expressed unigenes were also employed for the enrichment analyses of GO and KEGG pathway with adjusted p-value (q-value) < 0.05 serving as the standard for significantly enriched pathway.
Validation of differentially expressed Unigenes using realtime qRT-PCR Twenty (8 flowering and 12 non-flowering) top buds and ten (4 flowering and 6 non-flowering) lateral buds of saffron crocus with various corm weighst and bud lengths were used to validate differentially expressed unigenes using real-time quantitative reverse transcription PCR (qRT-PCR). Eleven differentially expressed unigenes between flowering and non-flowering samples were selected for validating key flower unigenes. All buds were ground in liquid nitrogen, and total RNA was prepared using an RNeasy@Plant Mini Kit. The Prime-Script II 1st Strand cDNA Synthesis Kit (TaKaRa, Japan) and SYBR Premix Ex Taq II (TaKaRa, Japan) were used for reverse transcription reaction and qRT-PCR assay. Specific primers of the chosen genes were designed using Primer Premier 5.0 software (Additional file 2: Table S2). PCR products were verified by dissociation curves, and data were normalized with endogenous reference tubulin gene to obtain ΔCt values. Water was used as a negative and quality control, and each sample was measured in triplicate.

Expression analysis of the flower-related genes in tissues and organs
The expression analysis of the flower-related genes in different tissues and organs was performed with qRT-PCR. Total RNA from the top and lateral buds (0.5-1 cm in length), the inner immature flowers (obtained from top bud when it grew to 1.5-3 cm in length), the corms, leaves, petals, stigmas, stamens and the remaining protective sheath of the full-bloom flowers, were extracted using an RNeasy@Plant Mini Kit and the following reverse transcription reaction and qRT-PCR assays were conducted according to the above description. The expression levels of flower-related genes in each sample were normalized to the tubulin gene to obtain ΔCt values. The top bud was used as a control sample, and the relative expression levels of target genes in the other samples were analysed using the 2 -ΔΔCt method: ΔΔCt = ΔCt other sample (Ct target gene-Ct tublin)-ΔCt control sample (Ct target gene-Ct tubulin).
Time course expression analysis of flower-related genes during the flower development Total RNA from four different stages of top buds from 20 g corms, including resting bud (1-2 mm in length), early stage of shoot growth (2-5 mm in length), late stage of shoot growth (5-10 mm in length), and stage of visually distinguishable flower organ formation (10-15 mm), were extracted using an RNeasy@Plant Mini Kit and the following reverse transcription reaction and qRT-PCR assays were conducted according to the above description.

Data availability
The raw data were uploaded to Sequence Read Archive (SRA) (http://www.ncbi.nlm.-nih.gov/) with a reference of PRJNA528829.

Long-length Transcriptome of saffron Crocus from PacBio Iso-seq
High-quality RNAs from top buds, tubers, pistils, stamens, petals and leaves of flowering saffron crocus were combined to acquire the PacBio Iso-seq libraries. Meanwhile, PacBio Iso-seq libraries of non-flowering saffron crocus were constructed using leaves, lateral buds, tubers, and top buds of non-flowering corms (20 g and 6 g). Multiple size-fractionated cDNA and cells (3 cells for 1-2 kb, 3 cells for 2-3 kb, 2 cells for > 3 kb) were prepared to construct flowering/non-flowering Iso-seq libraries separately. This approch avoids loading bias and obtaining more RNA sequences representing the gene expression profiles in flowering and non-flowering saffron crocus.
A total of 22.85 Gb of clean data were obtained from all sixteen cells with 1,325,207 raw polymerase reads and 23.9 billion nucleotides. After the adaptor and lowquality sequences were filtered, a total of 12,433,006 subreads were obtained, among which 7,178,336 and 5, 254,670 subreads were in the libraries of flowering and non-flowering saffron crocus, respectively (Additional file 2: Table S3). High quality ROIs were further generated from CCS after filtering with full passes and accuracy. The numbers of ROIs from the flowering saffron crocus libraries were 224,710 for 1-2 kb, 199,782 for 2-3 kb, and 106,171 for 3-6 kb, respectively, which were more than those of the corresponding nonflowering saffron crocus libraries (179,712 for 1-2 kb, 73,160 for 2-3 kb, 52,904 for 3-6 kb) (Additional file 2: Table S4). In total, 394,653 (74.4%) and 252,850 (82.7%) full-length non-chimaera reads (FL non-chimaera, fulllength reads with 3′ primer, 5′ primer and polyA reads after chimaera was filtered) were produced from ROIs of flowering and non-flowering saffron crocus libraries, respectively, with average lengths of 1223 bp, 2333 bp and 3512 bp in corresponding flowering saffron crocus libraries and 1188 bp, 2236 bp and 3322 bp in that of non-flowering saffron crocus libraries (Fig. 1b, Additional file 2: Table S4)).
After classification and correction by Clustering for Error Correction (CEC) and Quiver programs, 79,841 high-quality (Accuracy > 0.99) and 219,720 low-quality polished CCS were generated from ROIs. CCS were further corrected using the de novo assembly reads derived from Illumina RNA-seq. Ultimately, a total of 216,419 isoform level transcripts and 75,351 unigene transcripts were obtained after two-step CD-HIT classification of both flowering and non-flowering PacBio libraries. The length distribution of polished CCS, isoform and unigene is shown in Fig. 1c, with a majority of sequences ranging from 1 kb to 4 kb. The libraries of flowering and non-flowering saffron crocus were constructed separately, and the specific isoforms in each library and the differential expression profiles between flowering and non-flowering saffron crocus plants were obtained. The number of isoforms that expressed in both flowering and non-flowering saffron crocus was 174,369, while the number of isoforms that only expressed in flowering saffron crocus (30,188) were considerably more than those in non-flowering saffron crocus (11,862). These isoforms may provide a novel avenue to clarify the underlying molecular mechanism of floral development of saffron crocus.
Total 125 mRNAs derived from saffron crocus were reported on NCBI database at present. All the 75,351 full-length unigene transcripts were homologously aligned with them using BLAST. The results showed total 108 previously reported mRNAs were identified and matched with their highly homologous sequences in our data, with 86.4% coverage rate (Additional file 2: Table S5). Among them, 44 unigenes have a sequence identity of 99% or more and the identity of 88 unigenes were more than 95%, which suggested a full-length unigene database of saffron crocus with satisfactory coverage and accuracy was obtained in this study.

Saffron Transcriptome of short-reads from Illumina RNAseq
Fifteen Illumina RNA-seq libraries constructed from saffron crocus with different numbers of flowers (0-3) were sequenced to correct the polished CCS of PacBio Iso-seq and to quantify full length transcripts obtained from PacBio Iso-seq. After trimming process and screening with a high quality score, a total of 745 million clean reads were produced from all samples. Over 575 million short reads were successfully mapped back to the fulllength of PacBio Iso-seq with an average mapping ratio of 77.2% (Additional file 2: Table S6), which suggested that the full-length transcripts derived from PacBio Isoseq data method represented the majority of the genetic information of both flowering and non-flowering saffron crocus.
A total of 14,159 (21.9% of annotated unigenes) unigenes were associated with 34 pathways in KEGG pathway analysis. A high percentage of unigenes were assigned to "Translation" (10.3%) and "folding, sorting and degradation" (9.3%) of the genetic information process as well as "signal transduction" of the environmental information process (9.7%) (Additional file 3: Figure S2a).
A total of 64,562 unigenes (85.7%) were successfully matched to known sequences in at least one database. There were 99.5% matched unigenes in the NR database, 82.0% in SwissProt, and 72.0% in COG (Additional file 3: Figure S2d).
A total of 1193 GO terms were assigned to 33,117 unigenes (51.3% of annotated unigenes) with 454 biological processes, 159 cellular components and 580 molecular functions. In the class of biological processes, the top three GO terms were "metabolic process", "cellular process", and "single-organism process". In the cellular component, "cell" was dominant and then "cell part" and "organelle". In the class of molecular functions, a high percentage of the unigenes were enriched in "binding", "catalytic activity" and "molecular function regulator" (Additional file 3: Figure S2c).

CDS, SSR, and LncRNA prediction
The candidate coding sequence (CDS) in the PacBio transcript isoforms was analysed by retaining only open reading frames (ORFs ≥100 aa) using the ANGEL software. Both Arabidopsis thaliana and Phalaenopsis equestris genomes were used as the training sets. As shown in Fig. 2a, 50,197 CDS were obtained from the Arabidopsis thaliana genome with lengths ranging from 300 bp to 5400 bp and an average length of 1005 bp, while training with Phalaenopsis equestris genomes, ANGEL obtained a total of 289,377 predicted CDS with lengths ranging from 300 bp to 5400 bp and an average length of 1081 bp. Because saffron crocus is more closely related to orchids, more comprehensive information on encoded proteins would be obtained using orchid as the training set.
The PLEK workflow of lncRNA-pipeline was used to discriminate between coding and non-coding transcripts and then identify lncRNAs using PacBio data from species with no reference genome. To obtain more putative lncRNA candidates for saffron crocus, 216,419 isoform transcripts were used to predict lncRNAs in this study. A total of 72,603 (33.5%) PacBio non-coding transcripts were obtained and the length ranged from 194 bp to 6860 bp with an average length of 1367 bp. Similar to other species, the length abundance is concentrated at 500-1500 bp (54,296, 74.8%) (Fig. 2c).

Alternative splicing analysis and validation
Most mRNA precursors of eukaryotic genes produce only one mature mRNA that is thus translated to only one molecular protein. However, some mRNA precursors can produce different mRNA splice isoforms by different splicing sites, which is known as alternative splicing (AS). AS is an important mechanism of regulating gene expression and producing proteome diversity. At present, it is still challenging to reconstruct full-length splice isoforms using Illumina-based transcriptome assembly [37,38]. Splice isoforms with multiple introns make it difficult to identify alternative splicing using short read lengths, which were constrained by cufflink-based assemblies. One of the most important features of PacBio Sequencing is the ability to identify alternative splicing by directly comparing isoforms of the same gene without de novo assembly and thus avoiding artificial mistakes. Among the 75,351 unigenes identified in saffron crocus, 33.7% (25,400) have two or more isoforms. The number of AS events ranged from 2 to 217, and the distribution of AS events is shown in Fig. 3a. GO enrichment analysis showed that these AS genes were enriched in 120 pathways, with the top three being "Binding", "Heterocyclic compound binding" and "Organic cyclic compound binding" (Fig. 3b). It was interesting that the top 20 pathways were involved in various binding activities of molecular function ontology and a few AS genes were annotated with functions of catalytic activity and transport activity (Additional file 2: Table S7). Only two KEGG pathways, "Citrate cycle (TCA cycle)" and "Carbon metabolism", were significantly enriched by AS genes (Fig. 3c).
Furthermore, different splicing events between flowering and non-flowering saffron crocus were screened in PacBio-Seq libraries because AS is highly tissue-specific. A total of 42 AS events were identified in the PacBio-Seq libraries of flowering saffron crocus only, while 28 AS events were found in those libraries of non-flowering one (Fig. 3d, Additional file 2: Table S8). The AS events of three randomly selected unigenes, PB.174, PB.313 and PB.988, were validated using RT-PCR. Four, three, and two AS events were shown in the gel electropherogram (Fig. 3e), and the correctness of the amplification was proved by sequencing of the PCR products.

Differentially expressed Unigenes between flowering and non-flowering saffron Crocus corms
The expression levels of all unigenes from PacBio-Seq libraries in five groups of saffron crocus corms with different flowering phenotypes were quantified using an FPKM method based on the results of Illumina RNA-seq. In general, when FPKM < 0.1, the gene does not express; 0.1 ≤ FPKM < 3.75, the gene is considered to be expressed at a low level; 3.75 ≤ FPKM < 15, the gene is at a median level; FPKM ≥15, the gene is expressed at a high level. RNA-Seq2 software was used to screen differentially expressed unigenes between flowering and non-flowering saffron crocus corms, and the threshold of significance was determined as p-value < 0.05, false discovery rate (FDR) < 0.01 and log2 fold-change < − 1 or > 1.
The three top bud samples of flowering saffron crocus were compared to the three matched samples of nonflowering saffron crocus with cold treatment. A total 1677 DEGs (Additional file 2: Table S9), 70.8% (1187) of which were up-regulated, were screened with log 2 foldchanges ranging from − 9.4 to 12.6. All the DEGs were significantly enriched in 109 GO terms (q-value < 0.05), in which the function of the top five pathways was on binding and transcription factor activity with nucleic acid. These proteins play critical roles as transcription factors or at the transcriptional level in low temperature stress and flowering signalling pathways (Fig. 4a).. In addition, 11 unigenes were specifically highly expressed in three top buds of non-flowering saffron crocus with cold treatment while no expression was observed in any of the matched flowering samples. In contrast, 238 unigenes were highly expressed in three flowering samples exclusively (Additional file 2: Table S10). The DEGs with high fold-changes and strict consistency in biological repeat samples were ideal candidates for further largerscale validation. A total of 1086 significant DEGs (Additional file 2: Table S9), 81.7% (887) of which were up-regulated, were identified in top buds between flowering saffron crocus and non-flowering saffron crocus of 6 g corms, with log 2 fold-changes ranging from − 14.2 to 10.6. A total of 141 GO terms were clustered, and the top 5 pathways in each main ontology are shown in Fig. 4b. Within the cellular component ontology, the DEG-coding proteins were mainly located in the "macromolecular complex", "intracellular part" and "cytoplasm". For the biological process ontology, DEG-coding proteins were mainly involved in the transcription processes, such as "regulation of transcription, DNA-templated", "regulation of RNA metabolic process" and "regulation of nucleic acidtemplated transcription". In terms of molecular function ontology, DEG-coding proteins were binded with nucleic acid and functioned as transcriptional regulation factors  (Fig. 4b). In addition, 70 and 32 unigenes were exclusively highly expressed in all samples of large corms of flowering saffron crocus and small corms of nonflowering saffron crocus, respectively (Additional file 2: Table S11).
hree paired samples of lateral buds from the same corms of flowering and non-flowering saffron crocus were collected to screen key genes that regulate the flowering of lateral buds. A total of 267 DEGs (Additional file 2: Table S11), 83.5% (223) of which were down-regulated, were identified with log 2 fold-changes ranging from − 8.5 to 9.4. All 267 DEGs were assigned to 75 GO terms. The only GO terms in the cell component ontology was "nucleus". The most highly enriched pathway in biological progress ontology was "microtubulebased movement" (GO level 5), which was in the sub-grade GO level of both top 2 ("movement of cell or subcellular component", GO level 4) and top3 ("microtubule-based process", GO level 4) pathways, followed by the pathway of "DNA replication". In the molecular function ontology, most unigenes were highly enriched in a series of pathways from GO level 4 (hydrolase activity, acting on acid anhydrides) to GO level 7 (Nucleoside-triphosphatase (NTPase) activity). NTPases activation may regulate the nuclear export of mRNA and are known to be involved in defense signaling pathways and apoptosis regulation [39] (Fig. 4c).. In addition, 15 unigenes (Additional file 2: Table S11) were highly expressed in three lateral buds of non-flowering saffron crocus specifically but were not expressed in lateral buds of flowering saffron crocus. Only 7 unigenes were exclusively highly expressed in lateral buds of flowering saffron crocus (Additional file 2: Table S12). These 22 unigenes could be preferentially used to explore the molecular mechanism of the flowering process of lateral buds.
As for KEGG pathway analysis, although the DEGs involved in each experiment are different, "plant hormone signal transduction" was assigned by all the three comparable experiments, which even is the only pathway enriched in flowering vs. non-flower lateral buds (Fig.  4c). The changes of flowering phenotype caused by nutritional status (normal vs. small top buds), low temperature stress (normal vs. cold treatment) or natural variation (flowering vs. non-flower lateral buds) may be involved in different pathways, but they shared the same regulation network of plant hormones, including the response to hormone auxin, abscisic acid and cytokinin (for example, auxin-responsive protein IAA10-like, probable indole-3-acetic acid-amido synthetase GH3.5, auxin Both DEGs of normal vs. small top buds and normal vs. cold treatment were assigned to the KEGG pathways of "starch and sucrose metabolism" and "linoleic acid metabolism". It has been recognized for long time that temperature could modify carbohydrate metabolism of some plants by regulating starch and sugar levels [40]. Furthermore, altered carbohydrate metabolism impacted plant biomass production, as well as flower development. Actually, altering the expression of multiple sucrose metabolism-related genes in tobacco, such as UGPase, SuSy and SPS which were also found in DEGs of saffron crocus, not only enhanced primary growth but also altered flower morphology [41]. Based on the fact that flower numbers of saffron crocus were heavily depended on temperature and corm weight, the results of KEGG enrichment suggested that carbohydrate metabolism may play important roles in flower development of saffron crocus. There are many endogenous or exogenous factors affecting the flowering process of saffron crocus and driving different gene expressions. Genes shared in different pathways were initially speculated to regulate the flowering process of saffron crocus. All up and down regulated DEGs from the three comparable experiments (normal vs. cold-treated corm; big vs. small corm; top vs. lateral bud) were analysed using Venn software. Among all the DEGs, only six unigenes were shared in three groups. Unexpectedly, neither shared gene was identified in the upregulated nor downregulated DEGs (Fig. 5a). All the upregulated DEGs in the lateral buds of flowering saffron crocus were downregulated in the top buds of flowering saffron crocus and vice versa (Fig. 5b). Similar results were found in the DEGs shared in the two groups. For example, 78 shared DEGs (72 in two groups and 6 in three groups) were found in the lateral bud group and cold-treated top bud group, while only one shared DEG had consistent expression in both groups (Fig. 5a). Almost all shared DEGs had consistent expression in cold-treated top buds and normal/small top buds (475/476). The results suggested that there were different networks in regulating the flowering process of lateral buds of saffron crocus, and these networks may have been shared in the top buds.

Identification and expression analysis of flower-related DEGs in saffron Crocus
A total of 62 DEGs were functionally annotated as putative flower-related genes (Additional file 2: Table S13) and we analysed the expression patterns of all the transcripts using RNA-seq data (Fig. 5c). We identified 39 flower-related genes (assigned with predicted gene names) in the normal vs cold treatment group, 24 of which were significantly highly expressed in flowering corms. CsVRN (B3 domain-containing transcription factor) is a key gene involved in the vernalization pathway, which may affect flower development when the saffron crocus was treated with cold.
In the normal vs. small top buds group, 33 flowerrelated genes were identified and interesting, putative flowering time control protein (CsFCA) and squamosa promoter-binding-like protein 12 (CsSPL12), which are involved in autonomous and age pathways, were significantly highly expressed in the larger corms. The decreasing expression levels of CsFCA and CsSPLs may partly explain why the small corms (weight less than 6 g) cannot bloom. Furthermore, SPL12, strongly regulated together with CONSTANS (CO), were negatively regulated by miR156b and miR156h, suggesting micro-RNAs biogenesis or microRNA alteration maybe involved in the process of floral initiation [42].
Only five flower-related DEGs were identified in the flowering vs. non-flowering lateral bud groups and all of them were significantly lower in the flowering samples. There may be some novel genes, and even new regulatory pathways trigger the lateral bud bloom.
Because few flower-related gene sequences were submitted to the public database, among the 62 flowerrelated DEGs that we identified in this study, only PB.59337.3 (putative Apetala1-like MADS-box transcription factor, assigned as CsAP1-like) was matched to Crocus sativus. To clarify the flower-related gene regulatory network of saffron crocus, we plotted the schematic representation based on the known flowering gene network of Arabidopsis thaliana (Fig. 6): the flowering process of saffron crocus, similar to most other species, was involved in vernalization, gibberellins, photoperiod, autonomous and age pathways.
Validation of novel flower-related genes in larger-scale samples using real-time RT-PCR Except for the 62 genes that are homologous to known flowering gene sequences in public databases, we attempted to select and validate more novel flower-related genes in DEGs. Therefore, the fifteen samples were reclassified into six flowering samples (three flowering top buds and three flowering lateral buds) and nine non-flowering samples (three cold-treated top buds, three small top buds and three lateral buds) to identify DEGs using RNA-seq2 software with the same significance threshold as above.
Thirty-three DEGs were screened out, and only 3 unigenes were downregulated in the flowering group (Fig. 7a). Though the differences in DEG expression levels in top buds were more significant than those in lateral buds, the DEGs related to flowering processes were likely involved in downstream pathways or common nodes in flowering regulatory networks.
Ten DEGs related to the flowering process were selected to validate the key flowering genes using an additional 30 saffron crocus buds. The biological parameters of saffron crocus buds were collected, including the type of bud (top or lateral) number of flowers (0-2), length of bud (2 mm -5 cm), and fresh corm weight. Real-time RT-PCR was used to The expression pattern is shown in Additional file 4: Figure S3 and the results of the correlation analysis are shown in Table 1 Both genes were relatively highly expressed in immature flowers, which was similar to the expression pattern of previously isolated flower-related genes in saffron-  sequences were downloaded and screened for Phylogenetic analysis. The deduced amino acid sequences of PB20221.2 and PB38952.1 were aligned with the selected listing protein sequences using the multiple sequence alignment program ClustalX and phylogenetic tree were constructed with MEGA-X software using neighborjoining-method with 1000 bootstrap replicates. Aligning with homologous sequences from 38 different species, Phylogenetic studies revealed that PB20221.2 and PB38952.1 belonged to heat stress transcription factor family (HSF) and clustered in HSF-A2b sub-family (Additional file 6: Figure S5), indicating they were probably associated with the regulation of heat stress in flowering process of saffron crocus and further studies were needed to investigate their specific functions.

Discussion
Comprehensive, high-quality, full-length transcriptome sequences were presented for saffron crocus using SMRT technology in this study. A total of 75,351 fulllength saffron crocus unigenes were identified with an average length of 2049 bp in the flowering saffron crocus library and 1802 bp in the non-flowering library. Among these unigenes, 64,562 (85.7%) were functionally annotated. Based on RNA-seq data, Baba et al. identified 64, 438 transcripts with an average length of 609 bp and a functional annotation ratio of 58.5% (37,696) [19]. Subsequently, Jain et al. found 105,269 transcripts with an average length of 1047 bp [20], and 54% of them were functionally annotated. Basically, illumina RNA-seq can obtain higher coverage of genome and more transcripts than SMRT-Seq due to relatively higher sequencing depth. Our results further confirmed that SMRT-seq was more effective in recovering full-length transcripts [22]. The higher functional annotation ratio of unigenes in our study suggested that SMRT-seq data could provide higher accuracy and more effective information on saffron crocus transcriptomes.
Based on the PacBio-seq platform, 79,028 SSRs were identified in 34,895 unigenes (46.3% of total unigenes) in this study, which was considerably more than that predicted in the Illumina platform (16,721 SSRs identified in 13,407 transcripts). As popular molecular markers, SSRs are widely used for determining genetic variations, and the most important advantage is locus specificity, which is highly suitable for allopolyploid species, such as saffron crocus [43]. It should be noted that the status of genetic variation in saffron crocus is still controversial since several studies provided contrasting results [44]. Thus, the novel candidate SSRs of saffron crocus may provide pools to explore new alternative SSRs. In addition, the comprehensive CDS, lncRNAs and AS prediction of saffron crocus were reported for the first time in this study.
PacBio Iso-Seq libraries of flowering and non-flowering saffron crocus were constructed separately to identify as many candidate genes as possible without missing flowering suppressor genes that are expressed in non-flowering corms only. More unigenes were identified exclusively in flowering saffron crocus (30,188) than in non-flowering saffron crocus (11,862). This result suggested that saffron crocus plants might have additional genes to regulate floral development. The result might also indicate that suppressor genes for floral development may exist in some tissues of non-flowering saffron crocus only. The individual PacBio Iso-Seq library of flowering and non-flowering saffron crocus facilitated the clarification of specific AS in each library. AS functioned to regulate gene expression by introducing a premature termination codon andled to special RNA splicing isoform degradation. Increased AS was found in various plants with the application of highthroughput sequencing while the biological significance of  [45]. A total of 25,400 AS events in the full-length transcriptome were identified, and 42 and 28 AS events were identified in the flowering and non-flowering libraries, respectively. The underlying molecular mechanism for lateral buds of saffron crocus to flower was possibly different from top buds. The GO and KEGG pathway analyses of DEGs in lateral buds suggested potential genetic variation of flowering lateral buds, based on the enriched pathways of "cell cycle", "meiosis", "mismatch repair", "DNA replication" and "nucleoside-triphosphatase activity". The DEGs in top buds were found to function at transcriptional levels as well as in the signal pathways responding to hormones, nutrients, and other environments. Furthermore, the DEGs in lateral buds were completely different from top buds. The six shared DEGs were expressed in opposite trends between lateral and top buds. Saffron crocus multiplies by means of corms, but corm multiplication does not generate genome variation, except for occasional somatic mutations. However, phenotypic variations of flowers have been frequently observed in saffron crocus, as demonstrated by different numbers of stigmas, different aspects of tepals or blooming in lateral buds [30,46]. As the genetic variability in saffron crocus is still under debate, whether the phenotype variations are influenced by genetic variability or the environment warrants further investigation [47].
In this study, PB.20221.2 and PB.38952.1 were significantly correlated with flowering quantities and were suggested as flower-related genes. Phylogenetic analysis showed both genes encoded for HSF-A2b. Heat stress transcription factor family are the terminal components of a signal transduction chain mediating the activation of genes responsive to heat stress, including more than twenty members and classified as A, B and C in plants. HSF-A2 is a strong transcription activator through interacting with heat shock proteins (HSPs) during long-term heat stress [48]. Heat stress-regulated flowering was recently recognized as a floral transition pathway [49], together with the fact that HSP family was required for floral meristem formation [50]. Thus, as the central regulator of HSP expression, HSF-A2 may play important roles in the induction of inflorescence meristem formation. In fact, all the corms of saffron crocus can blossom only after enough heat accumulation of high temperature higher than 23°C, which also consistent with our results.
In conclusion, full-length transcriptomes of flowering and non-flowering saffron crocus were obtained using a combined NGS short-read and SMRT long-read sequencing approach. This method enables the generation of predicted comprehensive databases of AS, lncRNA, SSR, CDS, and DEGs in samples with phenotypic differences in flowering traits. These results were used to further identify genes related to flowering, including 62 genes homologous to known flowering gene sequences in public databases and novel flower-related genes,such as PB.20221.2 and PB.38952.1. Our study represents a first step to establish a reference full-length transcriptome for future studies of the gene atlas of saffron crocus and other species in the Iridaceae family. In the future, these genes, which are closely related to the number of flowers selected out in this study, maybe used as biomarkers for screening multi-flowering varieties and monitoring the optimal environmental conditions in field production.


text:
Direct and correlated responses from selection for early sexual maturity in Japanese quail

Base population after third generation of selection conclude total number of 788 birds at 6 weeks of age were used in this study and randomly assigned to two mating groups. First groupe was selected for early age at sexual maturity (EASM) and the second was as a control line (RBC) to study the direct and correlated responses for this direct of selection. The results can be summarized as follow: 1. EASM was highly significantly (P≤0.01) decreased from 39.54 days to 35.23 days. 2. The actual and the expected selection differentials were 2.52 and 2.29g / generation. 3. The actual response to selection in the (EASM) line was1.46g /generation. aHeritability of EASM estimated ranged from 0.29 to 0.35 bThe realized heritability ranged from 0.44 to 0.79. cCorrelated responses for EASM were as follow: dBW0, BW2, BW4, and BW6 increased significantly (P≤0.05) from 9.95, 40.58, 98.84 and 170.47g to 10.47, 47.01, 101.87 and 181.42g. eEgg weight (EW) was decreased (P≤0.05) from 10.74g to 10.04g. But Body weight at sexual maturity (BWSM) was increased from 136.32 g to 139.32 g. Total egg weight at first ten week (TEN10) decreased from 55.25eggs to 54.25egg. fFertility (FR %) was increased (P≤0.05) from 80.58% to 80.90%. and Hatchability (HA %) were increased (P≤0.05) from 70.42% to 70.92%. gMeat, dressing and giblets % were increased significantly (P≤0.05) from 42.33, 62.04 and14.27% to 43. 76%, 63.99 and 14.99%.

INTRODUCTION
The primary purpose of egg production selection trials is to improve production performance by increasing the number of birds reproduced from each hen (Bahie El-Deen et al. 1998).
The comparatively poor reproductive performance of birds selected for faster development rate is a frequent concern in all avian species studied. This poor performance is typically manifested as a reduction in overall egg production as well as an increase in the reproduction of eggs that are unsuited for certain environments, such as softshelled membrane and double-yolked eggs (Nestor andBacon, 1982 andStrong et al, 1978).
Selection experiments provide a framework for studying complicated trait inheritance and allow for the testing of theoretical predictions by comparing observations to expectations.
Short-term experiments can be used to estimate genetic variances and covariance, verify their consistency from various sources of information, and assess the magnitude of initial rates of response to selection, depending on the time scale (Varkoohi and Kaviani, 2014).From generation zero to generation four of selection, the current study aimed to quantify direct responses to selection for early sexual maturity and to connect those responses with specific productive and reproductive qualities. In addition, the genetic and phenotypic relationships between early sexual maturity and each of the attributes of growth, egg production, reproduction, and carcass are being investigated. The poultry sector is particularly interested in growth, and the Japanese quail could be valuable in body weight research. Because growth is age-dependent, genetic control of body weight is expected to change with age. Because the physiological and chronological ages of quail and hens are not the same, Different breeding strategies, like as selection, inbreeding, and bridization, can be used to optimize egg production and reproductive features.

MATERIALS AND METHODS
From May 2017 to September 2018, data for this study were collected on a flock of Japanese quail (Coturnix coturnix japonica) kept by the Department of Animal Production, Faculty of Agriculture, Al-Azhar University, Cairo, Egypt.
Breeding plan and management: A total of 788 birds (380 males and 408 females) were selected at random from the flock to serve as the study's parents at the age of six weeks.
In the basic population, the mating pattern was one male to two females, preventing full and half-sib mattings. The birds were randomly assigned to one of three mating groups: the first was selected for late age at sexual maturity (LASM), the second was kept as a randomly bred control line, and the third was kept as a random bred control line (RBC).
Family selection was done within the chosen line, however lower sexual maturity families (one male and two females) were chosen based on their divergence from the line mean.
Feed and water were freely available. According to N R C (1994), the minerals and vitamins were adequately given to meet the requirements (1994).
For the first week after hatching, the temperature was 37.5°C, then dropped 2-3°C weekly to 26-28°C in the fourth week of life till the end of the brooding period.
The breeding stocks were not subjected to any immunization or beak cutting operations.

Statistical analysis:
Harvey's least-squares and maximum likelihood computer programme (Harvey, 1987, model type 5) was used to evaluate data on age at sexual maturity, with each line, generation, and sex as fixed effects and sire and dam within sire as random variables. It was decided to use the following model: Yijmkln = µ + Gi + Lj + Sxm+Sk + Dlk + eijmkln. Where: Yijmkln = the observation of ijmkln th bird. µ = the overall mean, common element to all observations. Gi = the fixed effect of ith generation. Lj = the fixed effect of the jth line. Sxm = the fixed effect of mth sex. Sk = the random effect of kth sire. Dlk = the random effect of lth dam nested within a random effect of kth sire. eijmkln = the random error term.
Henderson Method 3 (Henderson, 1953) was utilized to estimate the genetic variance components for the different traits studied

1-Direct response:
In this paper, we describe the real methods of LASM for four generations of selection in the selected and control lines (Table, 1).Sexual maturity in the selected line (LASM) has recently grown, from 50.24 days at zero generation to 63.23 days in the fourth generation of selection. With a duration of 57.33 days on average. This outcome was in line with expectations.  Anthony et al., (1996). The control line, on the other hand, fluctuated arbitrarily from generation to generation, however there was a fluctuation from 50.24day at the zero generation to 50.84day at the fourth generation. With a daily average of 53.02 hours.

Selection differential:
For the LASM line, the expected and actual selection differentials, as well as the ratio between the actual and expected selection differentials (Actual/Expected), are provided in (Table, 2).From generation zero through generation four, however, had no significant differences were detected between the actual and expected selection differentials. This indicates that the LASM distribution closely resembles the features of a normal distribution, a Natural selection may have a favorable effect on EASM selection, as the ratio of actual to expected selection differentials ( Table 2) was greater than unity. Flowing for somewhat accurate prediction. Aboul-Hassan (2001b) and El-Fiky (2005) noted the similar pattern when selecting Japanese quail for high BW6, as did Abu-Mosallam (2014) when selecting Japanese quail for increased feed and water consumption during three generations of selection. Furthermore, for six generations, Marks and Lepore (1968) chose two lines of Japanese quail for body weight at 4 weeks of age.

Selection response:
In this table, the actual and cumulative selection responses are shown (Table 3).During the four generations of selection, the real response to selection in the line (LASM) declined from generation to generation.
Following the first generation of selection, the estimated response to selection was 1.82 days, then rapidly decreased as selection progressed to 1.62, 1.38, and 1.10 days after the second, third, and fourth generations of selection among the LASM line. With a 1.48-second average. This regularity of the selection response has been found in a number of published selection experiments. After four generations of selection, the cumulative response to selection in the LASM line was 5.92d. Heritability estimates for lately age at sexual maturity: Components of variance: The LASM line has heritability estimates for recently age at sexual maturity determined from components of variation, i.e. sire (paternal half-sibs), dam within sire (maternal half-sibs), and sire plus dam (full-sibs) ( Table 4). The dam component, on the other hand, is likely to be bigger in magnitude due to dominance deviations and maternal effects. Among the four generations of selection, heritabilities computed from the sire component (h 2 S) varied from 0.11 to 0.19, whereas those computed from the dam components (h 2 D) are often greater than those computed from the sire components, ranging from 0.26 to 0.33. This could be due to the substantial dam variance components obtained, as well as non-genetic influences such as dominance and maternal, which often cause h 2 D estimates to be much bigger than h 2 S estimates. Heritabilities calculated from full-sib components (h 2 S+D) were moderate to high in magnitude (0.23 -0.26). Debes (2004), Meky (2007), Abou Khadiga, et al. (2016), and Kaye, Joshua, et al. (2016) all noted the same tendency (2016). With a few exceptions, estimations of heritability based on female progeny weights were often greater than those based on male progeny weights. The estimated dam heritability estimates (h 2 D) are significantly bigger than the estimated sire heritability estimates, according to the unweighted mean of heritability estimates for male or female progeny weights (h 2 S).

Realized heritability:
The realized heritability was derived from the LASM line and shown every generation as the ratio between response to selection and selection difference (R/S) (Table,  4).The realized heritability estimates for recent age at sexual maturity dropped from 0.57 in the first generation to 0.30 in the fourth. However, he computed Realized heritability differently than Kaye, Joshua, et al., (2016). (0.23).The current estimates were consistent with those reported by Abu-Mosallam (2014), who calculated higher estimates of realized heritability for HFC4-6 when selection was applied, ranging from 0.70 to 0.77, which were higher than those reported by Aboul-Hassan (1997) and El-Fiky (2005), who reported a range of estimates (0.12-0.21) for realized heritability of BW6.

Correlated responses for lately age at sexual maturity:
Other qualities were harmed as a result of LASM selection. Weights and/or weight gains recorded at various ages and growth phases, as well as egg production, reproductive qualities, and carcass traits, may all be altered.
Tables from show the actual means of the relevant qualities associated with LASM selection (5 to 10).
Hassan (2000), Debes (2004), El-Fiky et al., (2000), and Abdel-Tawab (2006) noticed the same trend when selecting Japanese quail for increased EW, and Abu-Mosallam (2014) noticed the same trend when selecting Japanese quail for increased feed consumption among the first three generations of selection. Aboul-Hassan et al., 2016;El-Fiky et al., 2016. When it came to Japanese quail water and feed intake, they chose carefully. Body weight rose considerably (P0.05) by 0.15, 0.10, and 0.14 g each generation for the RBC line and insignificantly by 0.10 for the BW6 line at hatch, 2 weeks, and 4 weeks of age. However, after four generations of selection, body weight gain features in the LASM line (ADG0-2, ADG4-6, and ADG0-6) increased insignificantly (P0.05) from 2.32, 1.42, and 3.42 g/day in the first generation to 2.62, 1.69, and 3.64 g/day in the fourth. However, after four generations of selection, body weight gain at (ADG2-4) rose significantly (P0.05) from 5.20g/day in the first generation to 5.66g/day. The regression coefficients of the growth variables evaluated (weights and gains) on generation number have been used to determine correlated responses in the line chosen for LASM (table 5).

Egg production and reproductive traits:
After four generations of selection, egg weight (EW) grew considerably (P0.05) from 10.32g in the first generation to 11.62g, whereas BWSM increased insignificantly (P0.05) from 136.32 g in the first generation to 145.02 g. (Table, 6). TEW10 was considerably raised (P0.05) from 470.25g in the first generation to 464.54 g after four generations of selection, while TEN10 was insignificantly lowered (P0.05) from 54.25 egg un the first generation to 53.54g after four generations of selection. However, DEM fell from 9.40 g/day in the first generation to 9.30 g/day after four generations of selection (P0.05) ( Table, 6).
However, reproductive traits showed a significant decrease (P0.05) in FR percent from 82.58 percent in the first generation to 80.00 percent after four generations of selection, while HA percent decreased insignificantly (P0.05) from 70.52 percent in the first generation to 68.11 percent after four generations of selection (Table, 6), After four generations of selection, EEM percent grew significantly (P0.05) from 5.64 percent in the first generation to 6.79 percent. After four generations of selection, the LEM and TEM percents declined dramatically (P0.05) from 16.35 and 21.99 percent in the first generation to 16.69 and 23.48 percent. Farrag et al., (2008) discovered that fertility percentages in G0, G1, and G2 were 69.92, 64.71, 72.61; 59.87, 61.72, 73.21 percent and 60.58, 64.82, 75.90 percent, respectively. According to Samaha et al., (2010), fertility percentages in the base, first, second, and third generation were 94.11, 96.79, 94.05, and 95.34, respectively. Okenyi et al., (2013), reported that fertility percents were 76.29 percent for Go, 70.06 percent for G1, and 65.98 percent for G2 of selection. Furthermore, Abu-Mosallam (2014) found that egg weight (EW) increased significantly (P0.05) from 10.32 g in the first generation to 10.96 g after three generations of selection, and ASM increased significantly (P0.05) from 55.15 day in the first generation to 59.62 day after three generations of selection when he selected Japanese quail for increased feed consumption.
After three generations of selection, TEN10, TEW10, and DEM were reduced significantly (P0.05) from 54.25 eggs, 470.25 grams, and 9.40 grams per day in the first generation to 50.48 eggs, 458.25 grams, and 9.10 grams per day in the third generation.He also mentioned that reproductive traits such as FR percent and HA percent were reduced significantly (P0.05) from 82.58 and 70.52 percent in the first generation to 78.69 and 66.21 percent after three generations of selection, while EEM percent, LEM percent, and TEM percent were increased significantly (P0.05) from 5.64, 16.35, and 21.99 percent in the first generation to 6.95, 17.02, and 23.17 percent after three generations of selection. According to El-Fiky et al. (2016), total embryonic percent for HFC4-6 increased significantly (P0.05) from 6.61, 15.39, and 22.00 percent in the first generation to 8.25, 16.82, and 25.07 percent after four generations of selection. Meanwhile, Aboul-Hassan et al., (2016) found that early, late, and total embryonic mortality rates for HWC4-6 were reduced considerably (P0.05) from 8.25, 15.80, and 24.05 percent in the fourth generation to 6.63, 14.37, and 21.00 percent after four generations of selection. Faruque, Shakila, et al. (2016) calculated the mortality of four quail varieties: Japanese, White,Brown,and Black at 8.43,10.71,6.32,and 17.62, respectively. On the other hand, the control line's comparable values for egg production and reproductive attributes evaluated ranged from generation to generation with no discernible trend (Table, 6).

Carcass traits:
Meat, bone, giblets, and dressing percent attributes in the LASM line grew considerably (P0.05) from 44. 33, 7.57, 13.24, and 60.04 percent in the first generation to 45.41, 8.06, 13.95, and 61.94 percent after four generations of selection in the LASM line (Table, 7).The equivalent values for the control line among carcass qualities studied, on the other hand, were not significant (Table, 7).