The genome sequences of three Orchidaceae species, D. officinale (Zhang et al., 2016), P. equestris (Cai et al., 2014), and A. shenzhenica (Zhang et al., 2017) were downloaded from the NCBI database¹. The Cellulose synthase family proteins contain two different pfam domains: PF03552 for CesAs and PF00535 for Csls. The seed sequences of the two domains were then downloaded from the TIGRFAMs database². HMMsearch from the HMMER suite (version 3.1; Finn et al., 2011) was used to search for Cellulose synthase proteins in D. officinale, P. equestris, and A. shenzhenica with a cutoff E-value of 1e^––4. The molecular weight and theoretical isoelectric point (pI) of the proteins were predicted using the ExPASy website³. Transmembrane domains and subcellular locations were predicted using TMHMM⁴ and WoLF PSORT⁵, respectively.

The amino acid sequences of CesA/Csls were retrieved for nine species, including algae, mosses, lycophytes, monocots, and dicots: Chlamydomonas reinhardtii (Cre), Volvox carteri (Vca), Physcomitrella patens (Ppa), Selaginella moellendorffii (Smo), Oryza sativa (Osa), A. shenzhenica (As), D. officinale (Dof), P. equestris (Peq), and Arabidopsis thaliana (AT). The protein sequences were downloaded from the Pfam database⁶. MEGA 6 (Tokyo Metropolitan University, Tokyo, Japan; Tamura et al., 2013) was used to analyze the protein sequences and construct phylogenetic trees. First, the MUSCLE program of MEGA 6 was used for multiple sequence alignment, and then the maximum likelihood method (ML) method with the Jones--Taylor--Thornton (JTT) model was used to construct a phylogenetic tree with 1,000 bootstrap replicates. A partial deletion with a site coverage cutoff of 70% was used for gap treatment. The phylogenetic trees were visualized and modified using MEGA 6 and iTOL⁷.

The Gene Structure Display Server tool⁸ (v2.0; Hu et al., 2015) was used to analyze the gene structure of all the CesA/Csls identified in the three Orchidaceae species. Gene distribution in the genome of D. officinale was visualized using the TBtools program (Chen et al., 2020), and conserved domains were identified with the NCBI Web CD-Search tool⁹. MEME software¹⁰ (v4.11.0) was used to search for sequences of motifs of CesA/Csl proteins, with a motif window length from 10 to 100 bp (Bailey et al., 2009), with the maximum number of motifs set at 10, and where motifs present in at least three proteins were considered true motifs. Multiple alignments of CesA/Csls in the three orchid species were conducted using DNAMAN software (version 9; Lynnon Biosoft Company, Quebec, QC, Canada).

Two transcriptome sequencing analyses were performed using four different organs of a 3-year-old D. officinale and the stems of three developmental stages of D. officinale. The four different organs of the 3-year-old D. officinale used in the present study include flower, leaf, flower, and root. Three different D. officinale developmental stages were investigated: 1-year old (Y1), 2 years old (Y2), and 3 years old (Y3). The plants were grown in glasshouses at the Mulberry Field Station of Zhejiang Academy of Agriculture Science (Hangzhou, China). The tissues were collected and frozen in liquid nitrogen and stored at –80°C until use. For each tissue, five plants were treated as an independent biological replicate, and three biological replicates were performed. Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s instructions and sequenced on an Illumina HiSeq 2000 platform. Subsequently, the expression profiles of all D. officinale genes were obtained with FPKM (fragments per kilobase of exon per million fragments mapped) values using Cufflinks software¹¹ (v2.2.1) according to annotated gene models with a GFF file. The expression profiles of the CesA/Csl genes from each sample were analyzed using the HemI program¹² with the average hierarchical clustering method.

Stems of D. officinale plants at five different growth stages, i.e., 3 months (3M), 9 months (9M), 1 year (Y1), 2 years (Y2), and 3 years (Y3), were collected (three replicates for each sample) and dried in an oven at 105°C until constant weight. The 3-M and 9-M D. officinale seedlings were cultured on half-strength Murashige and Skoog (MS) (Murashige and Skoog, 1962) media containing 0.1% activated carbon, 2% sucrose, and 0.6% agar (pH 5.4) in a growth chamber (25.5 ± 1°C, 45 μmol m-2 s-1 irradiance, 12-h light/dark photoperiod, 60% relative humidity). Y1, Y2, and Y3 D. officinale plants were grown in glasshouses at the Mulberry Field Station of Zhejiang Academy of Agriculture Science (Hangzhou, China). The samples were shattered into fine powder independently by a mixing mill (MM 400, Retsch). The total polysaccharide was extracted using the water extraction and alcohol precipitation method, and the contents of total polysaccharides were measured using the phenol–sulfuric acid method as described by Wang et al. (2021). Glucose was used as a reference for subsequent calculations in which total polysaccharides (μg/g dry weight) = (A + 0.0037)÷7.981 × V1÷V2 × V3÷W × 1000 = 626. 49 × (A + 0.0037)÷W, where V1 represents the redissolved volume after alcohol precipitation (1 mL), V2 represents the volume of alcohol precipitation (0.2 ml), V3 represents the volume of water added during extraction (1 mL), and W represents sample weight in grams (1,000 g), serving as a coefficient converting milligrams to micrograms.

After drying, the samples from the D. officinale stems of the five growth stages were shattered into fine powder independently by a mixing mill (MM 400, Retsch). The constituting monosaccharides, including mannose, glucose, and galactose, were extracted using the GC–MS/MS method (Sun et al., 2016). Briefly, 20 mg of powder was diluted in a 500 μL solution and the extracts were centrifuged at 14,000 rpm at 4°C for 3 min. Then, the supernatants were mixed, evaporated, and freeze-dried. The residue was used for further derivatization. The small-molecule carbohydrates and 100 μL of a solution of methoxyamine hydrochloride in pyridine (15 mg/mL) were mixed together and incubated at 37°C for 2 h. Then, 100 μL of BSTFA was added and the mixture was maintained at 37°C for 30 min after being vortexed. The mixture was subsequently diluted and analyzed by GC–MS/MS according to the methods of Gómez-González et al. (2010) and Sun et al. (2016) with modifications. Agilent 7890B gas chromatograph coupled to a 7000D mass spectrometer equipped with a DB-5MS column (30 m length × 0.25 mm i.d. × 0.25 μm film thickness, J&W Scientific, Santa Clara, CA, United States) was employed for GC–MS/MS analysis of the monosaccharides.

Total RNA was extracted from all the samples of D. officinale stems at the five growth stages (3M, 9M, Y1, Y2, and Y3) using TRIzol reagent (Invitrogen, Carlsbad, CA, United States). The potential contaminating genomic DNA was eliminated with DNase I. The quality of the RNA samples was checked with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Beijing, China) and 1% denaturing agarose gels. The RNA was used as a template for first-strand cDNA synthesis using PrimeScript reverse transcriptase (TaKaRa Biotechnology, Dalian, China). Gene-specific primers were designed with the Primer Premier 5.0 program. The DnActin (comp205612_c0) gene was used as an internal standard for normalizing the gene expression data. The DofCesA/DofCsl expression levels were analyzed using a quantitative Real-time PCR (qRT–PCR) assay, which was performed by using a SYBR Green qPCR Kit (TaKaRa Biotechnology, Dalian, China) and a Stratagene Mx3000P thermocycler (Agilent, Santa Clara, CA, United States). The PCR program was as follows: 95°C for 5 min followed by 40 cycles at 95°C for 15 s and then 60°C for 30 s. The relative gene expression levels were calculated with the 2^–ΔΔCt method (Livak and Schmittgen, 2001). Three biological replicates were subjected to the analysis, each with three technical replicates. The expression levels in the different tissues were visualized using GraphPad Prism (v. 8.4.3).

Statistical analysis was performed, and the average values and standard errors of the three replicates were calculated. SPSS software (v. 16.0) was used to determine the significant differences in polysaccharide and monosaccharide contents between the different developmental stages using a one-way ANOVA procedure and post hoc analysis. The value p < 0.05 indicates a significant difference and is represented by an asterisk (*) in the figures, and p < 0.01 indicates a very significant difference and is represented by two asterisks (^**) in the figures. (***) indicates p value < 0.001 and (****) indicates p value < 0.0001.

To investigate the potential roles of Cellulose synthase superfamily genes in orchids, we performed a genome-wide identification and characterization of CesA/Csl genes in three sequenced Orchidaceae species: D. officinale (Zhang et al., 2016), P. equestris (Cai et al., 2014), and A. shenzhenica (Zhang et al., 2017). A total of 125 Cellulose synthase superfamily genes were identified in the three orchid species, which were designated ‘DofCesA/DofCsl’ for D. officinale, ‘PeqCesA/PeqCsl’ for P. equestris, and ‘AsCesA/AsCsl’ for A. shenzhenica (Supplementary Table S1). The CesA superfamily members in the Orchidaceae family were classified into seven families: CesA, CslA, CslC, CslD, CslE, CslG, and CslH according to the sequence similarities with that of A. thaliana and O. sativa (Figure 1, Supplementary Figure S1, and Supplementary Table S1).

We identified 54 CesA/Csl proteins from the D. officinale genome, which comprised 11 CesAs and 43 Csls (Table 1). The 54 proteins were further classified into DofCesA (11 members) and six different Csl families, i.e., DofCslA (19 proteins), DofCslC (5 proteins), DofCslD (8 proteins), DofCslE (7 proteins), DofCslG (3 proteins), and DofCslH1, according to the genome annotation information and sequence similarity with those of Arabidopsis and rice (Figure 1 and Table 1). P. equestris was identified with 37 CesA/Csl proteins, which were classified into PeqCesA (8 proteins) and six Csl families: PeqCslA (10 proteins), PeqCslD (10 proteins), PeqCslH (3 proteins), PeqCslE (4 proteins), PeqCslC1, and PeqCslG1. The A. shenzhenica genome was identified with 34 CesA/Csl proteins, including AsCesA (9 proteins), AsCslA (10 proteins), AsCslC (4 proteins), AsCslD (7 proteins), AsCslG (2 proteins), AsCslE1, and AsCslH1. These results indicated that the Cellulose synthase superfamily in D. officinale had greatly expanded, whereas that in the other two Orchidaceae species was similar to Arabidopsis (Persson et al., 2007), tomato (Song et al., 2018), and rice (Wang et al., 2010). Notably, although the number of DofCesAs was similar to that in other species, the DofCsl family had significantly expanded, especially DofCslA.

Three of the 125 CesA/Csl proteins (AsCslA2, PeqCslE4, and PeqCslH3) had extremely short amino acid sequences (less than 300 aa), even after reannotation via Softberry¹³. This may be due to genome sequencing errors or improper assembly. Thus, the three proteins were excluded from the subsequent statistical analysis. The molecular weights of the remaining 122 CesA/Csls ranged from 42.84 kDa (PeqCslA9) to 146.86 kDa (AsCslD5), with pI values ranging from 5.63 (DofCesA11) to 9.29 (AsCslA3 and DofCslA4) (Supplementary Table S1). The longest protein was AsCslD5 (1392 aa), whereas the shortest was PeqCslA9/PeqCslA10 (374 aa). The amino acid sequences of the CesA/Csls ranged from 415 aa to 1,180 aa in D. officinale, 515 aa to 1157 aa in P. equestris, and 438 aa to 1329 aa in A. shenzhenica. Except for four proteins, most of the CesAs and CslDs in the three orchid species were approximately 800–1,000 aa in length with 6–8 transmembrane (TMs) domains (Supplementary Table S1). The DofCslAs and DofCslCs were generally ∼500 aa and ∼700 aa in length, respectively, with 5 or 6 TM domains. Moreover, subcellular localization prediction showed that 94% of the CesA/Csls were localized on the plasma membrane. Nonetheless, two AsCslAs and four DofCslAs were predicted to localize on both the plasma membrane and the Golgi apparatus plasma. Detailed information on the CesA/Csl proteins in the three orchid species, including their gene/protein ID, CDS length, subcellular localization, TM domains, molecular weight, and pI value are shown in Table 1 and Supplementary Table S1.

To investigate the evolution of the Cellulose synthase gene superfamily, phylogenetic analysis was performed for nine representative plant and algal species, i.e., two algal species (Volvox carteri and Chlamydomonas reinhardtii), Physcomitrella patens, Selaginella moellendorffii, Oryza sativa, Arabidopsis, and three orchid species (D. officinale, P. equestris, and A. shenzhenica). A phylogenetic tree was constructed based on 273 amino acid sequences using MEGA 6 software and the maximum likelihood (ML) method with 1,000 bootstraps. The 273 CesA/Csls were classified into CesA and eight Csl families: CslA, CslB, CslC, CslD, CslE, CslF, CslG, and CslH (Figure 2A and Supplementary Figure S2). The lower land plants, i.e., the bryophyte P. patens and the lycophyte S. moellendorffii, are identified with representatives of the CesA, CslA, CslC, and CslD families, which is generally consistent with previous findings (Roberts and Bushoven, 2007; Yin et al., 2014). Each of the two green algae was identified as containing a single copy of the CslA/C-like gene (VcaCslC for V. carteri and CreCslA for C. reinhardtii), implying that the possible duplication of the common ancestor of CslA/CslC occurred in land plants after they diverged from green algae (Yin et al., 2009; Popper et al., 2011). Among the eight Csl families, CslA and CslC are distantly related to the other Csl families, supporting previous suggestions that the CslA and CslC family members evolved from a separate cyanobacterial endosymbiotic event (Yin et al., 2009). CslB was found only in the dicot species Arabidopsis, whereas CslH was specific to the four monocots (Figure 2C). Notably, CslG family members were found in all three Orchidaceae species.

The number of CesAs was rather conserved among angiosperms, with each generally having 8∼13 members (except for S. moellendorffii), whereas the members in the other families greatly varied. P. patens had the highest number of CesAs (13 members) and CslDs (10 members), whereas S. moellendorffii had the lowest (five SmoCesAs and three SmoCslDs). The CslCs were most abundant in S. moellendorffii (10 SmoCslCs); however, the family had greatly contracted in other species, especially in P. equestris (containing only PeqCslC1). Notably, the CslA family had significantly expanded in D. officinale (with 19 DofCslAs) compared with the remaining species (Figure 2B), which may have resulted in functional redundancy or innovation.

The sequences of conserved domains of the 125 CesA/Csl proteins in three Orchidaceae species were searched and analyzed. Five conserved domains: Cellulose synthase domain (cellulose_synt), zf-UDP, zf-RING_4, and two glycosyltransferase family domains (Glyco_transf_2, Glyco_trans_2_3) were found (Figure 3 and Supplementary Figures S3, S4). We found that the Cellulose synthase domain was present in CesA and CslD/E/G/H, whereas all the CslA/CslC proteins were found to contain two glycosyltransferase domains (Table 2), which support the CslA/CslC families evolved from an independent cyanobacterial endosymbiotic event (Yin et al., 2009). The number of CesA and CslD proteins that contain zf-UDP/zf-RING_4 domains varied among the three Orchidaceae species. For example, the zf-UDP domain was found in all nine AsCesA proteins but was absent from three DofCesAs and one PeqCesA. The motif patterns also varied among different species and families; however, we found species-specific motifs and similar conserved motif patterns within the same CesA/Csl family (Table 2, Figure 3 and Supplementary Figures S3, S4). These results suggested that CesA/Csl proteins among different families in Orchidaceae species might have different functional properties. Furthermore, multiple alignments of the predicted Cellulose synthase amino acid sequences showed that 18 of the 125 CesA/Csl proteins had no “D,D,D,QxxRW” integrated active site amino acid sequence (Supplementary Figures S5–S7), implying possible functional redundancy.

To investigate the sequence diversity of the CesA/Csl proteins in the three orchid species, putative motifs were identified using MEME Suite 5.0.2; a total of 10 conserved motifs were found (Figure 3C and Supplementary Figures S3, S4). Although the motifs varied among the different species and families, the CesA/Csl proteins in each family in the same species shared several unique motifs. For example, in D. officinale, motif4, motif8, and motif2 were all observed in the DofCesA, DofCslD, DofCslG, and DofCslE families but not in the DofCesA11 family (Figure 3C and Supplementary Figure S5A). The DofCslA and DofCslC family proteins shared five common motifs: motif4, motif7, motif1, motif6, and motif5; however, DofCslA10 lacked motif 4. Among these motifs, motif7 and motif1 were only observed in DofCslA and DofCslC families. In P. equestris, the conserved motifs in PeqCesA, PeqCslD, PeqCslE, PeqCslH, and PeqCslG were similar; except for four proteins, all contained motif1–motif7 and motif10 (Supplementary Figure S3). Moreover, motif9 was specific to members of the PeqCslA and PeqCslC families. In A. shenzhenica, motif1, motif5, and motif2 were uniformly observed in all CesA/Csl proteins, with the exceptions of AsCesA9 and AsCslA2. Members of AsCesA, AsCslD, AsCslE, AsCslH, and AsCslG also shared similar motifs (Supplementary Figure S4), and motif 10 was specifically present in the AsCslA and AsCslC proteins. Taken together, these results revealed species-specific motifs and similar motifs shared among certain families. Thus, CesA/Csl proteins among different families and species might have different functional properties.

To gain more information on the evolutionary patterns of Cellulose synthase genes, the gene structures of the CesA/Csls in D. officinale, P. equestris, and A. shenzhenica were analyzed (Supplementary Table S1). The intron–exon structure of the CesA/Csls was highly diverse among the different families within the same species (Figures 4A,B). Nonetheless, the same family members have similar exon/intron numbers among the three Orchidaceae species. For example, except for four genes, AsCesAs and PeqCesAs usually have 13 or 14 exons. The exon number in DenCesAs varied more significantly compared to the other two orchids and ranged from 8 to 16 exons. CslAs generally contain 5–10 exons, and most have 9 or 10 exons. The exon–intron structures and intron phases of the genes were relatively conserved within the same family (Figure 4B and Supplementary Figures S3, S4). In D. officinale, all of the DofCslCs had 5 exons and 4 introns. The gene structure varied the most within the CslD genes, from 0–7 or 8 introns/exons. These results suggested that possible functional diversity occurred within the same family and evolutionary conservation in Orchidaceae species.

Based on the physical positions of the Cellulose synthase genes in the D. officinale genome (Zhang et al., 2016), the distribution of 54 DofCesA/DofCsl genes was mapped to 30 scaffolds (Figure 4C and Table 1). Among them, 23 scaffolds had only one Cellulose synthase gene and two scaffolds each had two genes. Apart from the alternative splicing genes (indicated with green scale bars and asterisks), 14 tandemly duplicated genes (with two or more homologous genes ≤ 100 kb apart) were detected on six scaffolds (Figure 4C). A maximum number of 10 genes were detected on scaffold NW_021503689.1, all of which were DofCslA members (DofCslA10-19); after merging alternative splicing genes, there were four tandemly duplicated genes. DofCslA3, -4, -5 and DofCslA6, -7, -8 were also tandem duplicated genes located on scaffold NW_021319309.1 and NW_021319871.1, respectively. In addition, the DofCslEs also contain tandemly duplicated members that clustered on scaffold NW_02131965.1 (Figure 4C). Moreover, we found that the structures of alternative splicing genes, such as the number of exons and the length of introns, varied in different degrees among the splicing variants of the same gene (Supplementary Figure S8 and Supplementary Table S2). These results revealed that the expansion of the DofCesA/DofCsl genes mainly originated through gene duplication.

To analyze the roles of the CesA/Csl genes in D. officinale, their expression patterns were investigated in different organs in D. officinale at different developmental stages. Transcriptome sequencing was performed on four organs of D. officinale (i.e., the leave, stem, flower, and root) and the stems of three different stages (Y1, Y2, and Y3). The expression levels of the DofCesA/DofCsl genes are presented in heatmaps in Figure 5 and Supplementary Table S3). The results revealed diverse expression patterns of the DofCesA/DofCsl genes, suggesting functional divergence had occurred after gene expansion.

In the DofCesA family, three genes (DofCesA1, DofCesA5, and DofCesA8) were consistently expressed across different organs and developmental stages and were all significantly expressed in the flowers. Among them, DofCesA8 also had high expression levels in the roots and low expression in the leaves and stems; both DofCesA1 and DofCesA5 were expressed at low levels in the leaves, roots, and stems. Notably, DofCesA8 was also highly expressed in the Y1, Y2, and Y3 stems of D. officinale. The DofCslD and DofCslC family genes exhibited similar expression patterns, with low expression across the four organs and three developmental stages. The DofCslEs showed organ-specific expression patterns; for example, DofCslE1 and DofCslE4 exhibited moderate expression in the roots and leaves, respectively, whereas DofCslE6 was expressed in the leaves, roots, and stems and across all three growth stages. DofCslH1 was expressed in three organs (flowers, leaves, and roots), with the highest expression level in the roots. Moreover, the DofCslG family genes showed similar expression patterns; all exhibited higher expression levels in the flowers than in the roots, leaves, and stems. DofCslG1 and DofCslG2 also exhibited moderate expression in Y2 and Y3 stems.

CslAs are known for their participation in the synthesis of polysaccharides, especially mannans and glucomannans. Interestingly, the DofCslAs exhibited apparent organ- and development-specific expression in D. officinale. Most of these genes were predominantly expressed in the stems and/or flowers, with very low levels in the roots and leaves. For example, DofCslA2, DofCslA14, and DofCslA15 had significantly higher expression levels in the stems than in the other organs (especially DofCslA15); DofCslA2 and DofCslA14 were also expressed in the flowers. However, several genes had almost no expression in any of the four organs. In addition, DofCslA6 was specifically expressed in the flowers. Among the three growth stages, some genes showed similar expression patterns, such as DofCslA7 and DofCslA8, both of which were expressed at higher levels in the Y2 stems. Notably, DofCslA14 and DofCslA15 had significantly high expression in Y3 stems. However, five of the 19 DofCslAs showed low levels of or no expression in the Y1, Y2, and Y3 stems. DofCslA2 was expressed at moderate levels in the Y3 stems and at low levels in the Y1 stems, whereas DofCslA3 showed a consistent expression pattern across all three stages, with the highest expression levels occurring in Y2 stems. These results revealed the possible functional specialization of DofCslAs in the polysaccharide synthesis in D. officinale.

A previous study showed that polysaccharide accumulation in D. officinale occurred predominantly in the stems (He et al., 2015). To further investigate the polysaccharide synthesis and functions of Cellulose synthase superfamily genes in Orchidaceae species, we measured the contents of total water-soluble polysaccharides, and their component monosaccharides (glucose, mannose, and galactose) in the stems of D. officinale at five different growth stages, i.e., 3 months (3M), 9 months (9M), 1 year (Y1), 2 years (Y2), and 3 years (Y3) (Figure 6). The results showed that the total polysaccharide content varied significantly among the different stages (Figure 6D and Table 3); it was lowest in the 3M stems (∼18.06 mg/g) and highest in the Y3 stems (∼93.02 mg/g). The total polysaccharide contents varied slightly between the 3M and 9M stems and between Y2 and Y3 stems, but increased significantly from 3M/9M to Y2/Y3.

To determine the composition and contents of monosaccharides in the stems, the neutral monosaccharide composition was analyzed at the five growth stages. The results revealed that mannose was the most abundant neutral monosaccharide in the stems of D. officinale across all five growth stages, followed by glucose and galactose (Figure 6). The contents of glucose and mannose ranged from 2.40 to 20.01 mg/g and from 3.09 to 69.13 mg/g, respectively, with the highest detected in the Y3 stems and the lowest in the 3M stems. Mannose and glucose in the stems increased remarkably during D. officinale development, corresponding to the increases in total polysaccharides. Interestingly, the galactose content was highest in the 3M stems; however, this content decreased as D. officinale developed, lowest in the Y1 stems, and then increased slightly from Y2 to Y3.

To better understand the roles of CslAs in polysaccharide synthesis in D. officinale, qPCR analyses of the 19 DofCslAs were performed on D. officinale stems at five different growth stages: 3M, 9M, Y1, Y2, and Y3. The results for Y1, Y2, and Y3 were generally consistent with the transcriptome sequencing data. In general, most of the DofCslAs exhibited high expression levels in Y2 and/or Y3 stems (Figure 7 and Supplementary Table S4), which was in agreement with the accumulation of water-soluble polysaccharides and monosaccharides. Among them, five genes, i.e., DofCslA2, DofCslA6, DofCslA11, DofCslA14, and DofCslA15, exhibited significantly high expression levels in the Y3 stems (especially DofCslA15). Four genes, i.e., DofCslA3, DofCslA4, DofCslA7, and DofCslA8, were highly expressed in the Y2 stems compared with the Y3 stems. Notably, several genes, such as DofCslA2, DofCslA3, and DofCslA7, were also expressed in the 9M stems; DofCslA8 also exhibited moderate expression in the 3M stems. In addition, six DofCslAs, such as DofCslA16-19, showed low expression levels. The primer sequences for qPCR analysis are provided in Supplementary Table S5. Taken together, these results suggest that DofCslA genes may exhibit functional redundancy and may be specialized for activity at different growth stages.

The Cellulose synthase superfamily belongs to the glycosyltransferase-2 superfamily, which is widely distributed throughout the plant kingdom (Kim et al., 2007; Cantarel et al., 2009). In the present study, we identified a total of 125 Cellulose synthase superfamily members from the three sequenced orchids using two Pfam domains (PF03552/PF00535; Table 1 and Supplementary Table S1). Among them, 54 DofCesA/DofCsls, 37 PeqCesA/PeqCsls, and 34 AsCesA/AsCsls were identified in the genomes of D. officinale, P. equestris, and A. shenzhenica, respectively. The results revealed species-specific expansion of the Cellulose synthase superfamily in D. officinale.

In land plants, the Cellulose synthase superfamily can be further classified into a CesA family and nine Csl families: CslA-CslH and CslJ (Fincher, 2009; Yin et al., 2014; Schwerdt et al., 2015). Previous studies indicate that PF03552 encompasses the CesA family and seven Csl families (CslB/D/E/F/G/H/J), whereas PF03552 is absent from the CslA/CslC families (Yin et al., 2009). In the present study, the CesA superfamily genes in Orchidaceae were grouped into one CesA family and six Csl families: CslA, CslC, CslD, CslE, CslG, and CslH (Figure 1 and Supplementary Figure S1). The number of CesAs in the three orchid species was largely the same, while the number of Csls greatly varied, from 25 AsCsls to 43 DofCsls. Interestingly, the Csl proteins exhibited family-specific expansion among the three orchid species; for example, there were 19 DofCslAs, almost twice that in the other two orchids. The results of similar expansion were also observed in Csl families such as PeqCslDs and PeqCslHs. The expansion of CesA/Csl families likely resulted in functional redundancy or innovation. To further investigate the evolution of CesA/Csls, a phylogenetic tree based on 273 amino acid sequences from nine representative plant species was constructed (Figure 2). CesA/Csls were classified into one CesA family and eight Csl families: CslA, CslB, CslC, CslD, CslE, CslF, CslG, and CslH. These results are consistent with those of previous findings in which CslA, CslC, and CslD were found to be conserved in all land plants (Farrokhi et al., 2006; Yin et al., 2014). CslF and CslH may be restricted to monocots, and CslFs were absent from the three Orchidaceae species. CslG members were previously thought to be confined to eudicots (Fincher, 2009). However, CslG family members were found in all three Orchidaceae species, supporting the proposal that CslG should no longer be considered a dicot-specific family (Yin et al., 2014). In addition, our results strongly supported previous findings that the numbers of CesA/Csls vary greatly among different plant species (Paterson et al., 2009; Takata and Taniguchi, 2015; Song et al., 2018; Zou et al., 2018; Li et al., 2020). The number of CesAs is largely conserved among angiosperms, most of which have 8–13 members (except for S. moellendorffii), whereas the members in other families greatly varied among different species, indicating that extensive expansion and diversification have occurred.

We identified three proteins in P. equestris and A. shenzhenica that had extremely short amino acid sequences (less than 300 aa). The molecular weights of the remaining 122 CesA/Csls in the three orchid species ranged from 42.84 to 146.86 kDa. The longest protein was AsCslD5-1392 aa, which is shorter than the 2038 aa sequence in tomato (Song et al., 2018). Previous studies have indicated that cellulose is synthesized by plasma membrane-localized Cellulose synthase complexes with access to the cytosolic GDP-glucose pool (Taylor, 2008; Endler and Persson, 2011; McFarlane et al., 2014; Schneider et al., 2016). In the present study, subcellular localization predictions showed that most of the CesA/Csls were located on the plasma membrane.

The GT2 superfamily is characterized by conserved cytosolic substrate binding and catalytic residues. These residues are positioned in a loop between TM domains 2 and 3 and contain a D,D,D,QXXRW motif (Pear et al., 1996). In addition, plant CesA proteins also contain an extended N-terminal zinc-finger domain and two insertions in the catalytic loop, which are presumably involved in specific functions in higher plants, such as multimerization (Pear et al., 1996; Sethaphong et al., 2013). In the present study, five conserved domains were identified from the 125 orchid CesA/Csls (Figure 3 and Supplementary Figures S3, S4). We found that the Cellulose synthase domain was present in CesA and CslD/E/G/H, whereas all the CslA/CslC proteins were found to contain two glycosyltransferase domains, which support the CslA/CslC families evolved from an independent cyanobacterial endosymbiotic event (Yin et al., 2009). The number of CesA and CslD proteins that contain zf-UDP/zf-RING_4 domains varied among the three Orchidaceae species. For example, the zf-UDP domain was found in all nine AsCesA proteins but was absent from three DofCesAs and one PeqCesA. The motif patterns also varied among different species and families; however, we found species-specific motifs and similar conserved motif patterns within the same CesA/Csl family (Figure 3 and Supplementary Figures S3, S4). These results suggested that CesA/Csl proteins among different families in Orchidaceae species might have different functional properties. Furthermore, multiple alignments of the predicted Cellulose synthase amino acid sequences showed that 18 of the 125 CesA/Csl proteins had no “D,D,D,QxxRW” integrated active site amino acid sequence (Supplementary Figures S5–S7), implying possible functional redundancy.

Gene structure and duplication analysis can facilitate the understanding of gene evolution, structural divergence, and functional conservation within a gene family (Worberg et al., 2008; Xu et al., 2012; Nawaz et al., 2017). We found highly diverse structures among the different families but similar exon/intron numbers in the same family among the three orchid species (Figure 4B). The results revealed that there was functional diversity within the same family and that evolutionary conservation occurred among different orchid species. Moreover, a total of 14 tandemly duplicated genes were detected (Figure 4C and Table 1), indicating that the expansion of DofCesA/DofCsl genes might result from tandem duplication. These results provide insight into the evolutionary divergence and origins of CesA/Csl in Orchidaceae species.

Cellulose synthase superfamily genes are involved in the synthesis of the subunits of cellulose and hemicellulose (Burton et al., 2004; Scheller and Ulvskov, 2010; Takata and Taniguchi, 2015). The expansion and diversification of the plant CesA superfamily are linked tightly with major events in the evolution of plant and algal lineages, including multicellularity, terrestrialization, and vascularization (Popper et al., 2011). To understand the functions of the CesA/Csl genes in orchids, we investigated gene expression profiles in different organs of D. officinale and at different developmental stages (Figure 5). The results revealed diverse expression patterns of the DofCesA/DofCsl genes, suggesting functional divergence had occurred after gene expansion.

Three DofCesAs (DofCesA1, DofCesA5, and DofCesA8) were significantly expressed in the flowers of D. officinale. DofCesA8 also exhibited high expression in the roots and moderate expression in the leaves and stems. Except for DofCesA6, which exhibited low expression in the flowers, the expression levels of the remaining DofCesAs were minimal in all the tested organs. In addition, DofCesA8 was also highly expressed in Y1 stems and moderately expressed in the Y2 and Y3 stems. In the Arabidopsis genome, ten AtCesA genes were identified and found to be responsible for primary (CesA1-3, -5, -6, and -9) and secondary (AtCesA4, -7, and -8) cell wall synthesis (Taylor et al., 2003; Burton et al., 2004; Takata and Taniguchi, 2015). Mutations in secondary cell wall AtCesAs result in collapsed vasculature or irregular xylem phenotypes (Turner and Somerville, 1997; Taylor et al., 2003). Moreover, AtCesA1 and AtCesA3 also play roles in flowering; mutations in atcesa1 and atcesa3 result in gamete lethal phenotypes, whereas single-knockout AtCesA2, -5, -6, and -9 mutants exhibit more moderate phenotypes (Scheible et al., 2001; Caño-Delgado et al., 2003). The AtCesA1 was clustered in the same subgroup with DofCesA1, and AtCesA3 was clustered with DofCesA8. Therefore, the results, combined with the expression patterns, suggested that the three DofCesAs likely play roles in flower organ development, and that DofCesA8 may play additional roles in cellulose deposition in the stems of D. officinale.

In plants, the synthesis of hemicellulose occurs in the Golgi membrane by proteins encoded by Csl genes (Kaur et al., 2017). In land plants, CslAs and CslCs encode proteins that are involved in the synthesis of mannan and the β-1,4-linked glucan backbone of xyloglucan, respectively (Liepman et al., 2005, 2010; Cocuron et al., 2007). Moreover, CslD members have been reported to have mannan and Cellulose synthase activities, especially in tip-growing root hairs and pollen tubes of plant species such as rice and Arabidopsis. In rice, the OsCSLD1 gene is required for root hair morphogenesis (Kim et al., 2007). AtCslD3 also plays a role in Cellulose synthesis, as atcsld3 mutants exhibit defects in polarized growth of root hairs and abnormal distribution of cellulose and xyloglucan (Pauly et al., 2001; Park et al., 2011). CslCs are involved in the biosynthesis of the β-1,4-linked glucan backbone of xyloglucan (Cocuron et al., 2007). However, the expression of the DofCslD and DofCslC family genes was either low or absent from the four organs and across the three developmental stages, indicating minimal effects on polysaccharide synthesis in D. officinale. Interestingly, most of the DofCslAs were predominantly expressed in the stems and/or flowers. DofCslA2, -14, and -15 had significantly high expression levels in the stem; DofCslA6 was specifically expressed in the flowers. Moreover, DofCslAs were differentially expressed among the three growth stages. For example, DofCslA14 and DofCslA15 exhibited significantly high expression in Y3 stems, while DofCslA7 and DofCslA8 were expressed at higher levels in Y2 stems. Five of the 19 DofCslAs exhibited low expression or none in the stems of plants at all three growth stages. These results revealed functional divergence and possible redundancy of the DofCslAs involved in polysaccharide synthesis in D. officinale. The functional role of the CslB/E/G families remains unclear, but CslF/H genes have been suggested to be involved in the synthesis of mixed-linkage glucans (MLGs) (Prasad et al., 2005; Burton et al., 2006). We found that the DofCslGs all presented higher expression levels in the flowers than in the rest of the organs. DofCslG1 and DofCslG2 also had moderate expression in Y2 and Y3 stems. Our phylogenetic study showed that DofCslE had significantly expanded in D. officinale compared with other species (Figure 2). However, the expression of most DofCslEs was very low in different organs and growth stages. Thus, further investigations are needed to elucidate the exact roles of DofCslE/G/H family genes.

Containing approximately 1,450 species, Dendrobium is one of the largest genera of the Orchidaceae family (Zhang et al., 2016). Dendrobium species are characterized by their diverse growth habits and bioactive constituents with immunomodulatory hepatoprotective activities, such as dendrobine and polysaccharides (Ng et al., 2012). D. officinale stems contain abundant polysaccharides, primarily those consisting of GM and GGM (Xing et al., 2014). The high amount of polysaccharides in D. officinale stems supposedly helps maintain osmotic pressure and improve drought tolerance, as this species is epiphytic in natural habitats (Zotz and Tyree, 1996; Wu et al., 2009). CslAs encode proteins with both mannan and glucomannan synthase activity (Liepman et al., 2005, 2010). In Arabidopsis, the stems of csla knockout mutants have no glucomannan, indicating an exclusive role in mannan biosynthesis of the CslA family genes (Goubet et al., 2009). To determine the roles of the 19 DofCslAs in polysaccharide accumulation in D. officinale stems across different growth stages, we performed qRT–PCR and measured the content of total polysaccharides and their monosaccharide components in the stems of 3M, 9M, Y1, Y2, and Y3 plants. The total polysaccharide content was lowest in the 3M stems and highest in the Y3 stems (Figure 6 and Table 3). The content varied slightly between 3M and 9M stems but increased significantly in the mature stems. Monosaccharide composition analysis showed that mannose was the most abundant component across all development stages, supporting previous findings in Dendrobium species (Meng et al., 2013). The amounts of glucose and mannose in the stems also experienced remarkable increases from 3M to Y3. In agreement with the changes in total polysaccharide and monosaccharide component contents, most of the DofCslAs exhibited high expression levels in the Y2 and/or Y3 stems (Figure 7). Among them, five genes exhibited significantly high expression levels in Y3 stems (especially DofCslA15); four genes were highly expressed in Y2 stems. Thus, these genes likely play roles in the biosynthesis of mannan and glucomannan in the stems of D. officinale. In addition, six genes showed low expression levels, suggesting that possible redundancy occurred after gene expansion. These results imply that the DofCslAs may experience functional specialization with respect to polysaccharide accumulation in different growth stages.