Five hybrids from Den. Sect. Spatulata, designated H1–H5, were used to assess floral organ phenotypes. These hybrids were initially acquired from orchid markets or local farms. Hybrids H1–H4 likely originated from genetic crosses involving Den. stratiotes, Den. lasianthera, Den. Panavee, and other varieties. Hybrid H5, also from Den. Sect. Spatulata, has an untraceable parentage and was initially obtained from an orchid farm in Hainan Province, China. All the Dendrobium plants were subsequently cultivated in our greenhouse. Based on greenhouse observations, it takes approximately 25 days for flower buds to develop from their initial separation from clustered growth to blooming. During this period, flower bud length increases from an average of 9.5 mm to 30 mm. Consequently, the developmental stages of Den. Sect. Spatulata were classified as follows: Stage 1, buds shorter than 16.00 mm after separation from clustered growth; Stage 2, buds between 16.00 mm and 23.00 mm; Stage 3, buds exceeding 23.00 mm; Stage 4, blooming initiation. The observation of plant growth and flowering, along with RNA-seq sampling, began in August 2022. For each developmental stage, four replicates of petals and sepals were collected, with approximately 1g of tissue per replicate. Three replicates were used for Illumina short-read sequencing, while the fourth replicate was pooled and used for long-read sequencing. All samples were immediately frozen in liquid nitrogen and stored at − 80 °C.

Total RNA from the sepals and petals of Den. Sect. Spatulata was extracted using Tiangen DP411 (Tiangen, Beijing, China) according to the manufacturer’s instructions. RNA concentration and integrity were analyzed using a NanoDrop 2000 spectrophotometer (Thermo, MA, USA) and an Agilent 2100, LabChip GX (Perkin Elmer, Shanghai, China). Qualified RNA was stored at -80°C for later use.

A total of 1 μg of RNA per sample was used as the input material for library preparation. The mRNA was purified from the total RNA using HieffNGS DNA selection Beads (Yeasen, Shanghai, China). Sequencing libraries were generated from the purified mRNA using the Hieff NGS Ultima Dual-mode mRNA Library Prep Kit for Illumina (Yeasen, Shanghai, China) following the manufacturer’s recommendations with unique index codes. The library quantification and size were assessed using a Qubit 3.0 fluorometer (Life Technologies, CA, USA) and an Agilent 2100, LabChip GX (Perkin Elmer, Shanghai, China). Subsequently, sequencing with a paired‐end sequencing length of 150 bp (PE150) was performed on the Illumina Novaseq 6000 platform (Illumina, CA, USA) by Biomarker Technologies Co., Ltd. (Beijing, China). The raw sequencing data have been deposited in the NCBI Sequence Read Archive (SRA) under BioProject accession number PRJNA1129577.

A total of 600 ng RNA was used to generate the PacBio long-read library. Full-length (FL) cDNA was synthesized using the NEBNext Single Cell/Low Input cDNA Synthesis and Amplification Module (New England Biolabs, MA, USA) and subsequently amplified. Size selection of the FL cDNA (1–6 kb) was performed using the BluePippin system (Sage Science, MA, USA). The ends of the FL cDNA were repaired and the hairpin sequencing adapters were ligated using the SMRTbell^® Express Template Prep Kit 2.0 (Pacific Biosciences, CA, USA). Library quality and quantity were evaluated using an Agilent 2100, LabChip GX (Perkin Elmer, Shanghai, China) and a Qubit 3.0 fluorometer (Life Technologies, CA, USA). The polymerase-bound SMRTbell libraries were sequenced on the Pacific Biosciences Sequel II System (Pacific Biosciences, CA, USA). The raw sequencing data have been deposited in the NCBI SRA under BioProject accession number PRJNA1129577.

The raw subreads were analyzed following the Iso-Seq3 pipeline (https://github.com/PacificBiosciences/IsoSeq). The pipeline included three initial steps: generation of the circular consensus sequence (CCS) subreads, classification of the FL reads, and clustering of the full-length non-chimeric (FLNC) reads. Polished CCS subreads were generated, using CCS v6.2.0, from the subreads bam files with a minimum quality of 0.9 (–min-rq 0.9). The default minimum number of FL subreads (n = 3) required to generate CCS for a zero-mode waveguide (ZMW) was used. FL transcripts were determined when the sequences had the poly(A) and the 5′ and 3′ cDNA primers. Lima v2.1.0 and isoseq3 refine were used to remove the primers and poly(A) tails, respectively. The clustering algorithm ICE was used to obtain high-quality FL consensus sequences. High-quality FL consensus sequences were classified with the criterion of a post-correction accuracy above 0.99. High-quality full-length (FL) transcripts obtained through Iso-Seq were processed using CD-HIT with an identity threshold greater than 0.99 to eliminate redundancy (Fu et al., 2012).

Simple sequence repeats (SSRs) of the transcriptome were identified using MISA (Beier et al., 2017). TransDecoder (https://github.com/TransDecoder/TransDecoder/releases) was used to identify candidate coding regions within transcript sequences based on several criteria. First, a minimum-length open reading frame (ORF) must be present in the transcript sequence. Second, the log-likelihood score, similar to the score computed by GeneID, must be greater than 0. Third, the coding score should be highest when the ORF is scored in the first reading frame compared to the other five frames. Additionally, if one ORF is fully encapsulated by the coordinates of another, the longer ORF is reported, though multiple ORFs may be reported per transcript to account for operons, chimeras, and other scenarios. Optionally, the putative peptide may match a Pfam domain with a score above the noise cutoff. Iso-Seq™ data were used to perform an all-vs-all BLAST with high-identity settings, and alignments that met the following criteria were considered products of candidate alternative splicing (AS) events: the alignment must contain two High-Scoring Segment Pairs (HSPs) in the same forward or reverse direction, one sequence must be continuous or have a small overlap (less than 5 bp), and the other sequence must exhibit a distinct “AS Gap”. The continuous sequence should nearly fully align with the distinct sequence, and the AS Gap must be greater than 100 bp and located at least 100 bp away from the 3’ or 5’ ends.

Raw sequencing data (in fastq format) were initially processed using in-house Perl scripts. During this step, clean reads were obtained by filtering out reads containing adapters, poly-N sequences, and low-quality reads. Concurrently, various quality metrics, including Q20, Q30, GC-content, and sequence duplication levels, were calculated for the clean data. All subsequent analyses were performed using high-quality clean data.

The PacBio sequencing data generated in this study served as the reference transcriptome. Kallisto software was used to construct an index file (Kallisto index transcriptome.fa -i transcriptome.indx), followed by the quantification of gene expression levels for each sample (Kallisto quant -i transcriptome.index -o outfile R1.fastq.gz R2.fastq.gz). Gene expression levels were estimated using the number of transcripts per kilobase million (TPM) mapped. Differentially expressed genes (DEGs) were identified using the DESeq2 package (Love et al., 2014). Gene Ontology (GO) enrichment analysis was conducted with GO annotation data using the clusterProfiler 4.0 package (The Gene Ontology Consortium, 2019; Wu et al., 2021). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using KOBAS 3.0 (Bu et al., 2021).

To validate gene expression, total RNA was extracted from petal and sepal samples collected at the same developmental stages as those used for the transcriptome sequencing. RNA extraction was performed using the Total RNA Extraction Kit for polyphenol- and polysaccharide-rich plants (R4150; Magen Biotech, Guangzhou, China). Complementary DNA (cDNA) was synthesized from the extracted RNA using the Hifair^® II 1st Strand cDNA Synthesis Kit (11119ES60; Yeasen Biotech, Shanghai, China). Quantitative real-time PCR (qRT-PCR) was conducted with Hieff^® qPCR SYBR Green Master Mix (Low Rox Plus) (11202ES08; Yeasen Biotech, Shanghai, China). The primer sequences used for gene expression validation are provided in Supplementary Table S1 .

The samples were organized into dehydration cassettes and processed through a series of ethanol solutions: 75% for 4 hours, 85% for 2 hours, 90% for 2 hours, and 95% for 2 hours, followed by two rounds of absolute ethanol for 1 hour each. They were then immersed in an ethanol-xylene mixture for 5-10 minutes, followed by two changes of xylene for 5-10 minutes each. Paraffin infiltration involved treating the samples in a paraffin-xylene mixture at 62°C for 2 hours and two changes of melted paraffin at 62°C for 4 hours each. The paraffin-infiltrated tissues were transferred to an embedding station, placed into molds with melted paraffin, and cooled on a -20°C plate. The solidified paraffin blocks were trimmed, sectioned into 8-10 µm slices using a microtome, and floated on a 40°C water bath before being mounted on slides and baked at 60°C. The sections were deparaffinized in xylene and rehydrated through a series of ethanol solutions before toluidine blue staining for 2-5 minutes. Finally, the sections were rendered transparent by immersion in clean xylene, mounted with neutral balsam, and observed under a microscope.

Den. Sect. Spatulata usually show spirally twisted floral segments and asymmetrical flowers, resembling antelope horns. Figure 1A shows the flowers of five hybrids (H1-H5) from Den. Sect. Spatulata that were preserved in our germplasm resources. All these hybrids possess curled sepals. The petals of H1 and H2 were curled up into typical spirals, and those of H3 and H4 were gently curled up. However, H5 exhibited contrasting flower organs as the sepal was curled up, whereas the petal remained almost flat ( Figure 1A ). The natural and flat lengths of petals and sepals were measured, and the curvation rate was obtained by dividing flat length by natural length. Except for H3, all the sepals had roughly higher curvation rates than the petals ( Figure 1B ). The curvation rate of the H5 sepals was substantially higher than that of the petals and was also higher than most of the sepals in the other four hybrids. As can be seen from the four developing stages of H5 flower organs, the sepals were curled the most in Stage 4 ( Figure 1C ). The plant height of H5 exceeded 1.2 m ( Figure 1D ), and its biomass supported a higher annual booming frequency. These characteristics make H5 of Den. Sect. Spatulata an ideal plant for studying the molecular mechanism of sepal curvature.

Full-length transcriptome sequencing of H5 in Den. Sect. Spatulata generated a total of 62.62 Gb of data, in which 619,157 CCS reads were obtained. The distribution of CCS read lengths shows that the number of CCS reads was 124,515 for 0-1 kb, 252,760 for 1-2 kb, 118,699 for 2-3 kb, 45,406 for 3-4 kb, and 77,776 for >4 kb ( Figure 2A ; Supplementary Table S2 ). Further, a total of 472,576 FLNC reads were obtained by filtering adaptor and low-quality sequences, accounting for 76.33% of the CCS reads. The number of FLNC reads from the library was 159,049 for 0-1 kb, 188,585 for 1-2 kb, 71,077 for 2-3 kb, 24,050 for 3-4 kb, and 29,864 for >4 kb, respectively ( Figure 2B ; Supplementary Table S2 ). Since the FLNC reads in the cDNA library contained duplicate isoforms, the similar FLNC read sequences were clustered by the IsoSeq module in SMRTLink. This step generated 166,706 consensus isoforms, containing 166,650 high-quality (HQ) consensus isoforms with an accuracy higher than 0.99. The average lengths of HQ consensus isoform reads were 632, 1,461, 2,421, 3,425 and 5,090 bp for the 0-1 kb, 1-2 kb, 2-3 kb, 3-4 kb and >4 kb libraries, respectively ( Figure 2C ; Supplementary Table S2 ). Finally, a total of 94,258 non-redundant transcripts were obtained by removing the redundant transcripts with the CD-HIT software using the clustering threshold of “-c 0.99”.

After removing the redundant transcripts, BUSCO (version 3.0.2) was used to assess the completeness of the non-redundant transcripts. Among the total of 1,440 BUSCO groups, 80.8% (1,163) were complete BUSCOs, including 453 single-copy BUSCOs and 710 duplicated BUSCOs ( Figure 3A ). These results indicated that most of the non-redundant transcripts were assembled with high quality. Alternative splicing is crucial for generating transcripts with different combinations of exons, which can thus be translated to different proteins and result in biological diversity. To identify the potential alternative splicing events, all the non-redundant transcripts were used to perform pairwise sequence alignment using the criteria described in the methods section. Finally, a total of 2,762 potential alternative splicing events were obtained ( Supplementary Table S3 ). These data may support further studies focusing on the identification of robust alternative splicing events.

In order to obtain annotation information for the transcriptome, all the non-redundant transcripts (94,258) were aligned to nucleotide or protein databases. A total of 82,254 transcripts were finally annotated, of which 66,964, 57,071, 23,660, 47,497, 57,628, 55,649, 81,586, 69,340, and 81,477 transcripts were annotated in the GO, KEGG, COG, KOG, Pfam, Swiss-Prot, TrEMBL, eggNOG and nr databases, respectively ( Supplementary Table S4 ). The nr homologous species distribution showed that 86.08% of the nr annotated sequence was consistent with the published sequences of Dendrobium catenatum, and the second had 6.65% sequence identity with Phalaenopsis equestris ( Figure 3B ). These results suggested that the full-length transcriptome sequencing obtained robust data and can be used as a reference transcriptome in the downstream analysis.

SSR is widely applied in the construction of a genetic linkage map and analysis of genetic diversity in orchids. MISA software was used to identify SSR from the transcripts with over 500 bp. A total of 39,157 SSR sites were obtained, including 23,961 mono nucleotides (p1), 6,871 di nucleotides (p2), 4,519 tri nucleotides (p3), 199 tetra nucleotides (p4), 85 penta nucleotides (p5), 145 hexa nucleotides (p6), 3,263 compound SSRs (c), and 114 compound SSRs with overlapping positions (c*) ( Supplementary Table S5 ). The top three high-density SSR sites were p1, p2, and p3 with 147.35, 42.25, and 27.79 SSR sites per Mb, respectively ( Figure 3C ). The coding sequence and protein sequence of all the transcripts were predicted using TransDecoder software (v5.0.0). A total of 81,862 ORFs were obtained, of which 49,209 (60.11%) were complete. The number of proteins with complete ORFs was the highest at the length of 0-100 amino acids, and gradually decreases as the length of protein increases ( Figure 3D ).

Illumina high-throughput RNA sequencing (RNA-seq) was used to quantify the gene expression levels of 12 sepal samples and 12 petal samples, both contain 3 replicates of 4 developing stages of H5 flower organs. Generally, the number of clean reads ranged between 20,360,805 (petal1-1) and 26,910,018 (petal4-2), with an average of 22,982,401. The percentages of Q30 ranged between 92.78% and 94.09%, with an average of 92.78% ( Supplementary Table S6 ). Of these clean read pairs, 45.54% to 48.13% were uniquely mapped to the reference transcriptome, and 40.34% to 43.96 were mapped to multiple loci ( Supplementary Table S6 ). Both the uniquely mapped reads and the reads mapped to multiple loci were counted by Kallisto (v0.50.0) software and the counts were normalized to TPM. The TPM value of each sample was used to assess the total expression status and sample correlations. Most of the samples showed a similar expression distribution pattern ( Figure 4A ). Only the TPM value of the samples in Stage 4 (petal4_1 to petal4_3 and sepal4_1 to sepal4_3) was slightly lower than other samples, possibly due to the phenotypic difference between Stage 4 and other stages. This result was consistent with the principal component analysis (PCA) which revealed a general difference between petal and sepal samples, and the Stage 4 samples obviously clustered together ( Figure 4B ). Meanwhile, the sample correlation analysis showed a higher correlation coefficient between the replicated samples than that between non-replicated samples ( Supplementary Figure S1 ). These results demonstrate that the RNA-seq data were sufficient and reliable for subsequent differential expression analysis.

Although the curling phenotype was most prominent in Stage 4, we hypothesize that this trait is the culmination of regulatory events driven by specific genes across the floral developmental stages. To identify the underlying genes involved in sepal curvature, we conducted a bulk analysis by comparing sepals and petals across 12 samples from the four stages. This analysis revealed a total of 821 DEGs between sepals and petals. Notably, 598 DEGs (72.8%) were upregulated, while 223 DEGs (27.2%) were downregulated ( Figure 5A ). These DEGs displayed four distinct expression clusters: downregulated DEGs were grouped into clusters I and II, while upregulated DEGs formed clusters III and IV ( Figure 5B ). Gene cluster I in petals and cluster IV in sepals showed higher expression during the early stages of flower organ development (Stages 1 to 3). In contrast, gene cluster II in petals and gene cluster III in sepals peaked at Stage 4. Particularly, gene cluster III in sepals demonstrated a gradual increase in expression from early to late stages, suggesting a significant role for these DEGs in sepal morphogenesis.

GO and KEGG pathway enrichment analyses were conducted to examine the functional characteristics of the DEGs. Several GO terms were significantly enriched, including positive regulation of transcription by RNA polymerase II, carbon fixation, and pollen tube guidance ( Figure 6A ; Supplementary Table S7 ). Semantic similarity clustering was applied to group-related GO terms. We identified two subtrees with key terms including carbon peptide cross-linking binding and myosin contractile filament disassembly, which contained the majority of the enriched GO terms ( Figure 6B ). Peptide cross-linking, myosin II filament organization, and actomyosin contractile ring were frequently enriched in these two subtrees. Myosin II filaments play a critical role in regulating both the length and alignment of actin filaments in the contractile rings during cell division (Arima et al., 2023). These findings suggest that the DEGs involved in sepal morphogenesis regulation may also participate in cytokinesis.

KEGG pathway enrichment analysis revealed several significantly enriched pathways, including biosynthesis of secondary metabolites, fatty acid degradation, and starch and sucrose metabolism, all of which play important roles in plant development ( Supplementary Figure S2 ; Supplementary Table S8 ). Metabolites from the phenylpropanoid biosynthesis pathway show organ-specific and developmental patterns (Vogt, 2010). Additionally, pathways related to fatty acid metabolism and amino acid metabolism, including tyrosine metabolism, cysteine and methionine metabolism, cyanoamino acid metabolism, and alpha-Linolenic acid metabolism, were significantly enriched in our data. These pathways and their derivatives are known to be involved in floral organ development (Chakraborty et al., 2022). Notably, the DEGs enriched in pathways such as carbon metabolism, fatty acid degradation, glyoxylate and dicarboxylate metabolism, and starch and sucrose metabolism may work together to regulate plant energy storage, which is crucial for cytokinesis.

Genes can be annotated with multiple GO categories, and identifying key DEGs that are linked to several GO categories can provide deeper insights into their functional roles. To explore this, we constructed linkage network maps to illustrate the connections between the key GO terms and the associated 30 DEGs ( Figures 7A–C ). Among the 30 DEGs, 27 were upregulated, while only 3 were downregulated in sepals ( Figure 7D ). Notably, PB_transcript_51188 and PB_transcript_53627 were associated with diverse GO biological process (BP) categories, including the strigolactone biosynthetic process, reductive pentose-phosphate cycle, and peptide cross-linking ( Figure 7A ). These two DEGs also exhibited similar expression patterns, with high expression observed in Stage 1, Stage 2, and Stage 4 sepals ( Figure 7D ). Recent studies have shown that cross-linking proteins play a role in regulating actin filament network organization and force generation in yeast (Hill et al., 2024). These findings suggest that PB_transcript_51188 and PB_transcript_53627 may contribute to sepal morphogenesis, potentially by regulating the actin filament network.

In addition, two other DEGs, PB_transcript_96143 and PB_transcript_21432, were found to share the same BP, MF, and CC categories ( Figures 7A–C ). The proteins they encode possess calmodulin-dependent protein kinase activity and are involved in myosin II filament organization at the actomyosin contractile ring subcellular location. These two DEGs also showed high expression during sepal development, particularly at Stage 4 ( Figure 7D ). Calmodulin-dependent protein kinase regulates enzymatic activity and cytoskeletal organization by phosphorylating proteins such as myosin light chain kinase (Tokumitsu and Sakagami, 2022). During cytokinesis, the accumulation of actin and myosin II filaments at the cell equator forms a contractile ring that helps divide the cell into two daughter cells (Arima et al., 2023). These results further support the involvement of the key DEGs in the process of cytokinesis.

To validate the RNA-Seq results, we conducted qRT-PCR analysis on six selected genes, with two genes representing each of the three clusters shown in Figure 7D . PB_transcript_44935 and PB_transcript_131531 showed consistent downregulation in sepals throughout all four developmental stages ( Figure 8A ). Conversely, PB_transcript_21432, PB_transcript_44935, and PB_transcript_161752 were significantly upregulated in sepals across all stages. Notably, PB_transcript_96143 displayed a unique expression pattern, with suppression at Stage 3 followed by marked activation at Stage 4 in sepals. Overall, these expression patterns align well with the RNA-Seq data ( Figure 7D ), indicating a high level of concordance between the qRT-PCR and RNA-Seq findings. This consistency supports the reliability of the expression patterns identified across these conditions. To reveal whether the cell number and cell size were different between sepals and petals, the tissue sections of sepals and petals were visualized and compared in Figure 8B . The cell number and size were approximately the same in the upper and bottom layers of the petals. By contrast, the cell number from the upper layers was much higher than those in the bottom layers of the sepals. Meanwhile, the cell size in the upper layers was much smaller than that in the lower layers ( Figure 8B ). A similar cell distribution pattern was also observed in H2 sepals, whose sepals are even more curled than H5 ( Figure 1A ). These results support the hypothesis that the curled sepals are possibly caused by the mechanical force from the increased number of cells.