The Chromosome-Level Assembly of Ramie (Boehmeria Nivea L.) Genome Provides Insights into Molecular Regulation of Fiber Fineness

ABSTRACT Ramie (Boehmeria nivea L.), belonging to Urticaceae, is principally used for fabric production. It is a well-known natural fiber material for ancient clothing. Despite its important position and application value, the understanding on genetic regulation mechanism of fiber quality is limited. Here, we generate a chromosome-scale, high-quality reference genome of ramie, in which, approximately 90.2% of the assembled sequences have been anchored to 14 pseudochromosomes. Totally 27,664 protein-coding genes are predicted which cover 268.24 Mb region of the genome. Comparative genomic analysis reveals that 2,047 and 796 gene clusters expand and contract, respectively, underlying significant genes in plant hormone signal transduction and cellulose/lignin biosynthesis pathways. An integrative analysis combining quantitative trait loci (QTL), comparative transcriptomic data, and cytological experiments unravels the molecular regulatory mechanism of ramie fiber fineness, especially the critical regulating role of ethylene. This study would lay a solid foundation for the research of molecular biology in ramie and provide valuable reference for the improvement of high-quality fiber varieties.


Introduction
Bast fiber, the earliest textile used by the ancients for making clothes, is a symbol of civilization. Ramie (Boehmeria nivea L.), also called China grass, is one of the oldest and most important natural fiber crops worldwide. The superior fiber obtained from ramie stem is long and highly durable, pure white in color, silky in texture, and has a high degree of hygroscopicity and superior heat dissipation properties, which makes ramie being a highly versatile and useful natural raw textile material. The products of bast fiber are affected by its quantitative characteristics including length, strength, and fineness. At present, the ramie fiber fineness, which might be affected by multiple QTL (Topdar et al. 2013;Zeng et al. 2019) and numerous genes (Xie et al. 2020;Zhang et al. 2020), can't meet the requirements of high-quality textiles. It is well known that ramie fiber fineness is inversely proportional to fiber yield, which hinders the improvement of fiber fineness by traditional breeding technology. For a long time, only a small number of genetic markers and one genetic map of ramie have been published. The molecular mechanism of ramie fiber fineness regulatory is still to be revealed.
Previously, we have assembled a ramie genome of about 341 Mb, while another publication showed another version of ramie genome with about 266 ~ 270 Mb (Luan et al. 2018;Wang et al. 2021). More precise genome sequence assembly is conducive to the functional analysis of ramie genes. Here, we present a high-quality chromosome-scale genome assembly of B. nivea using a combination of shortread sequencing, single-molecule real-time sequencing, and a high-density genetic map. The detailed information of the ramie genome would be helpful for understanding the molecular basis of its fiber fineness development. We further performed a comparative transcriptomic analysis with high and low fiber fineness ramie varieties to investigate the molecular regulatory mechanism of fiber fineness.

Plant materials and growth conditions
The ramie variety (Zhongzhu No. 1) used for the earlier draft assembly was used for whole-genome sequencing in this study (Luan et al. 2018). Two ramie varieties with high (H) and low (L) fiber fineness were used for transcriptomic analysis. In the stage of 14 (T1), 24 (T2), 34 (T3), 44 (T4), 54 (T5) days after germination, the bast of middle plant stem was picked for transcriptome sequencing. Three independent biological replicates were set for each variety. Each replicate was collected from the same ramie plant at each sampled stage.

Genome sequencing, assembly, quality assessment, and annotation
In total, approximately 10 mg high-quality genomic DNA was extracted from young leaves of Zhongzhu No. 1. Two 20 Kb independent PacBio single-molecule real-time (SMRT) bell libraries were sequenced by the PacBio RSII sequencer (Pacific Biosciences, CA, USA) and 75.5 Gb of previous short insert size paired-end Illumina clean data (under BioProject ID: RJNA417329) was incorporated into this genome assembly for error correction. The PacBio long subreads was assembled by using FALCON (version 0.5.0) (Chin et al. 2016). The assembled consensus sequences were polished with Quiver algorithm using PacBio reads. Pilon (version: 1.21) was then adopted for genome correction with the published short paired-end reads generated from Illumina HiSeq Platforms (Luan et al. 2018). Reads mapping, Benchmarking Universal Single-Copy Orthologs (BUSCO V5) evaluation, transcripts alignment, and LTR assembly index (LAI) were used to evaluate the quality of the genome assembly, respectively. The gene sets of the genome were annotated combined using denovo, homology and RNA assisted approaches by the Maker pipeline. The predicted genes were functionally annotated through BLAST searches against seven public protein databases, including NCBI non-redundant, Swissprot, KEGG, KOG, TrEMBL, Interpro, and GO databases.

Phylogenetic analysis
We used 1,155 single-copy identified gene families for phylogenetic analysis. Protein sequences of single-copy gene families were performed with multi-alignment using MUSCLE (Edgar 2004). Supergenes of each species were produced by concatenating multi-alignment gene families. Oryza sativa was the out-group, supergenes were fed into RAXml with parameters (-m PROTCATJTT) for phylogenetic tree construction. We then used MCMC module of the Phylogenetic Analysis by Maximum Likelihood (PAML) package to estimate the divergence time of species. MCMC runs were conducted for 1,000,000 generations (1 sample/100 generations) and the first 250 samples were burned in. Calibration times were marked by the data from the TimeTree database.

Expansion and contraction of gene families
We evaluated the changes of gene families based on phylogenetic tree and gene families using CAFÉ (version 2.1) (De Bie et al. 2006). Expanded genes were then functional annotated by KEGG enrichment.

Genetic mapping and QTL mapping
A ramie F 1 hybrid progeny population (253 lines) and its parents Zhongzhu No. 1 and Hejiangqingma were used for QTL analysis. The extracted DNA from leaf samples was digested with EcoR I (New England Biolabs, USA), and then sequenced on an Illumina Hiseq 4000 platform (Illumina, San Diego, CA, USA) using paired-end reads (150 bp) sequencing strategy in BGI-Shenzhen, China. After filtering the low-quality reads and adapters, the first (left) paired-reads with a length of 136 nt were clustered and genotyped by the software package of Stacks (version: 1.23). The double pseudo-test cross strategy was used for genetic map construction using LepMap2 software. The RAD markers were mapped to the ramie genome sequence using Blat v.35 with E value <1e × 10 −10 . After filtering with parameters: identity ≥95% and aligning rate >95%, markers with more than one best hit in different contigs were eliminated. The remained markers were used to anchor the contigs to the consensus genetic linkage map with the ALLMAPS software. The Multiple QTL Mapping (MQM) algorithm in MapQTL6.0 software was used for QTL analysis. The QTLs were identified according to genome-wide logarithm of the odds (LOD) scores, and a significant level was marked at a level of p < .05 with the threshold values calculated using 1,000 permutations.
Plant height, stem diameter, bark thickness, fiber fineness, and ramet number were measured in two harvest seasons (October 2017 and June 2018). Fifteen ramets per plot were selected randomly to record the phenotypic data using band tape or vernier caliper. Fiber fineness was measured using the following method. The fiber is derived from the bast stripped, degummed, and fiber fineness was examined with optical fiber diameter analyzer systems (OFDA 100 and OFDA 2000; Wald, Switzerland) according to methods described by Brims and Hornik (2000).

Ethylene treatment experiment
Zhongzhu No. 1 was used in an ethylene treatment experiment. Ramie shoots of 15 cm in length were planted in a hydroponic device and cultured in 1/2 Hoagland solution under the following conditions: 500 μmol m −2 s −1 light; 12 h light/12 h dark period; 30°C in the day and 25°C during the night; 70% relative humidity. When roots of the seedlings grew to 8 cm, ethylene treatment (treated with 0, 50, 100, 200, 300, and 400 mg L −1 ethephon, respectively) was performed. Ethylene solution was sprayed on leaves and stems of ramie seedlings, once every one day, for three consecutive times. Three replicates were set for each treatment concentration, and each replicate contained 15 individual plants. After 7, 13, 23, 33, 43 days of treatment (counted from the first spray of ethephon), stems were collected for making slices. Paraffin sections methods according to Cai and Lashbrook (2006), and fibrocyte measurement followed as: five different regions of each splice were chosen to capture pictures using OLYMPUS DP73. Length and width of fiber cells were measured using ImageJ software.

Transcriptome analysis
Total RNA was isolated from the phloem tissue in different developmental stages. The cDNA libraries were constructed with an Illumina Paired End Sample Prep kit (Illumina, San Diego, CA, USA), quantified by TBS380 (Picogreen, Invitrogen, Carlsbad, CA, USA) and sequenced by BGISEQ-500 platform (MGI, Shenzhen, China). Paired-end reads with 100 bp length were obtained. Then reads that contained a high proportion of N, low-quality and contaminated by adaptors were filtered. Clean reads were mapped against gene models using Bowtie2 (version: v2.2.5). RSEM (v1.2.12) was used to evaluate and calculate expression quantity with default parameters. We detected differentially expressed genes using DEGseq with Fold Change > = 2 and Adjusted P value < = 0.001 among different groups. Differentially expressed genes were performed function enrichment and metabolic pathway analysis. In addition, expression patterns of differential expression groups were clustered by pheatmap, a package of R.

Genome sequencing, assembly, and annotation
We obtained a total of 30.6 Gb (~90-fold coverage with the estimated genome size) single-molecule real-time sequences from 22 Pacbio RS II cells (Supplemental Table S1). The final assembly of ramie genome yielded 268.24 Mb contigs and contained 522 contigs (N50: 1.36 Mb) (Table 1) with the largest contig length of about 6.14 Mb, a 47-fold improvement (Luan et al. 2018). To anchor and orient the contigs, we constructed a linkage map, with 253 F 1 mapping population and two parents using RAD sequencing (Supplemental Figures S1-2, Supplemental Table S2), which contains 13149 RADbased DNA markers consisting of 5,796 loci covering 2,934.3 cM on 14 consensus genetic linkage groups ( Figure 1, Supplemental Table S3). A total of 11289 markers were used to anchor the contigs to a consensus genetic linkage map. In total, 305 contigs with 241.85 Mb length were anchored to 14 chromosomes with Pearson's correlation coefficients ranging from 0.942 to 0.997, illustrating the high congruence between genome assembled sequences and genetic map data (Table 1; Supplemental Table  S3, Supplemental Figure S3-9).
To evaluate the genome assembly quality, the assembled transcript sequences were aligned against the assembled genome. A total of 97.57% of transcripts were covered by genome and the mapping rates were 93.42% and 99.13% with >90% and >50% identity in one scaffold, respectively (Supplemental Table 4). The BUSCO evaluation showed high quality of genome assembly with 98.1% completeness of the 1,614 genes in the Plantae BUSCO dataset (Simao et al. 2015) (Supplemental Table 5), which was more complete than that of 96.0% in previously reported genomes (Wang et al. 2021). In addition, the LAI score of ramie genome is 14.47 (Supplemental Figure 10), which also proved the high quality of this updated genome. A total of 27,664 protein-coding genes (Supplemental Table 6) were predicted. On average, the genes encode transcripts of 1,233.8 bp with 5.09 exons and exon length is 242.15 bp, which is longer than the C. capsularis, C. olitorius, and M. notabilis genome (Supplemental Table 7). Among these predicted genes, 88.99% could be functionally annotated (Supplemental Table 8) and 94 miRNA, 578 tRNA, 197 rRNA were identified (Supplemental Table 9). Features of ramie genome and differentially expressed genes distributions. Intra-genome syntenic blocks (a), gene expression between high fiber fineness and low fiber fineness in different stages (T1-T5 respectively represent b-f), GC content (g), repeat content (h), LTR density (i), gene density (j), and the chromosomes (k).

Genome evolution and gene family expansion
Comparative genomic and evolutionary analysis was performed between ramie and 14 other sequenced plant genomes, including representatives from the Monocots (Oryza sativa), Asterids (Solanum tuberosum), Rosids (Vitis vinifera, Linum uitatissimum and so on) (Supplemental Figure  11). The phylogenetic analysis was applied to 450 genes including single-copy and one-to-one orthologous genes, the results supported the placement of ramie with mulberry and cannabis in the Urticales and showed that B. nivea diverged withM. notabilis ~48.2 million years ago (Mya) and C. sativa (the nearest family species), ~64.0 Mya (Supplemental Figure 12). The gene families were identified as unique, shared, and expansion among the 15 plants (Supplemental Figures 11-12). We found 2,047 and 795 gene families that, respectively, expanded and contracted in the B. nivea genome compared with other 14 plant genomes (Supplemental Figure 12). Among these gene families, 69 gene families exhibited significant (p-value ≤.01) expansions, which are mainly enriched in physiological processes, such as biosynthesis of secondary metabolites and so on (Supplemental Figure 13, Supplemental Table 10).

QTL analysis of fiber fineness and fiber yield related traits
Fiber fineness is one of the key factors to evaluate the quality of ramie fiber, but it is inversely proportional to fiber yield (including plant height, ramet number, stem diameter, and bark thickness), which is the bottleneck in the breeding of high-quality fiber ramie varieties. Based on the high-density linkages mentioned above and the phenotypic data (Supplemental Figure 14), QTLs of fiber fineness and fiber yield-related traits were analyzed. In total, 58 QTLs of fiber fineness were identified with LOD value >3.0, which explained more than an average of 10% phenotype variance (Supplemental Figure 15, Supplemental Table 11). There were four QTLs that simultaneously controlled fiber fineness and ramet number, while two QTLs that simultaneously controlled fiber fineness and bark thickness (Supplemental Figure 15). Ten QTLs that controlled plant height and seven QTLs that controlled stem diameter, were found to simultaneously control fiber fineness. The linked QTLs synchronously controlled fiber quality and fiber yield in ramie identified in our study may explain why ramie fiber fineness is inversely proportional to fiber yield, which provided valuable reference information for the breeding of fiber quality improvement varieties in ramie.

Identification of key pathways and related genes regulating ramie fiber fineness
To reveal the molecular mechanism of fiber fineness regulatory, we did RNA-seq for 30 samples of ramie phloem fibers from two ramie varieties with high fiber fineness (H) and low fiber fineness (L) at 5 developmental stages (T1, T2, T3, T4, and T5) (Figure 2a, Supplemental Table 12). Differentially expressed genes (DEGs) were identified (Supplemental Figure 17A-D) and KEGG analysis was performed to investigate the detailed function of DEGs (Supplemental Figures 18-22). The plant hormone signal transduction pathway was specific in T1, T2, and T3, which mainly included auxin, brassinosteroid, ethylene, and cytokinine signal transduction (Figure 3a-c, Supplemental Figures 18-20). Auxin was more active in H at the early stage of fiber development judging by its higher number of up-regulated genes. Surprisingly, many genes in photosystem I, photosystem II, cytochrome b6/f complex, photosynthetic electron transport, and F-type ATPase were down-regulated in T2 and T3 (Figure 3d-h), suggesting that there was a weak photosynthetic system in H, which may also lead to a weak photosynthesis. Previously, there was evidence showing that photosynthetic rate was inversely proportional to fiber fineness in ramie (Peng and Yang 1986), and our study underlies the molecular regulatory mechanisms at the transcriptomic level for the first time.
The pathway of starch and sucrose metabolism was specifically enriched in T2, T3, and T4 ( Figure  3i, Supplemental Figures 18-20 Table 13), suggesting the major role of these genes in regulating ramie fiber fineness. We also predicted 18 cellulase genes in ramie genome, which might be an important factor in degradation of ramie fiber. Among these genes, cellulase.7, and cellulase.9 were reduced expressed in H, while cellulase.13, cellulase.15, and cellulase.17 were induced, underlying the main role of these genes in regulating ramie fiber fineness.
Lignin is also a main component of the fiber cell wall, which determines the size of the fiber cell so as to determine the fiber fineness. Compared with L, C3 H, HCT.8 and HCT.9 were silence expressed in H (Figure 4c). Differently, 4CL.13, 4CL.14, 4CL.9, F5 H.6, COMT.16, COMT.17, COMT20, and COMT.9 were dominantly expressed in H, highlighting possible high-lignin bast fiber in ramie which may suppress the expansion of fiber cells and result in high fiber fineness. There were 20 lignin-related genes could map on chromosome 3 distributed near the fiber fineness QTLs Ecor116209 and Ecor124469, but only four genes (BnID-D21, BnID-D2, BnC3 H-D190, and BnC4H2-D2) exhibited differential expression between H and L (Figure 2c), indicating that these lignin-related genes showed functional redundance in ramie fiber development. BnC3 H-D190 and BnID-D21 was, respectively, up-and down-regulated in H at the five stages, while BnC4H2-D2 was up-regulated at T2-T4 and BnID-D2 was up-regulated at T1-T4, indicating gene functional division in ramie fiber development. (d)

Ethylene involves in the regulating of ramie fiber fineness
Previous studies have suggested the critical regulating role of ethylene in cotton fiber cell elongation (Qin et al. 2007;Qin, Hu, and Zhu 2008;Shi et al. 2006), and there is also evidence that ethylene can increase ramie fiber fineness (Ouyang 2016;Wang et al. 2012). The core synthases in ethylene biosynthesis include S-adenosyl-L -methionine synthase (SAMS), 1-aminocyclopropane-1-carboxylic acid synthase (ACS), and 1-aminocyclopropane-1-carboxylic acid oxidase (ACO) (Figure 5a). We identified two SAMS, seven ACSs, and four ACO genes in the ramie genome (Supplemental Table 14, Figure 5e). The conversion of SAM to 1-aminocyclopropane-1-carboxylic acid (ACC) by ACS is generally considered to be the first committed step (Yang and Hoffman 1984) and also regarded as a rate-limiting step in the ethylene biosynthetic pathway (Prescott and John 1996). The dominant expression of BnACS1, BnACO2, BnACO1, BnACO2, BnACO3, and BnACO4 in H suggested the higher ethylene biosynthesis capacity in H, which resulted in high fiber fineness (Figure 5b, Supplemental Figure 23A, B). The result of the ethylene treatment experiment showed that the size of fiber cell was inhibited under ethylene treatments (Figure 5f-h, Supplemental Figure 24). These results confirmed that ethylene played a role in regulating the fineness of ramie phloem fiber. After ethylene synthesis, it is perceived by receptors located in the membrane and its signal transduced through transduction machinery to trigger specific biological responses (Wang, Li, and Ecker 2002). The main nodes and related genes of the ethylene signal transduction pathway have been well characterized (Figure 5c). We identified 5 ETR, 3 CTR, 2 EIN2, 6 EIN3, 2 EBF1, 2 EBF2, 7 ERF1, and 12 ERF2 genes in the ramie genome (Supplemental Table 14, Figure 5e). BnETR2, BnCTR1, BnEIN3-3, BnEIN3-5, BnERF1-1, BnERF1-2, BnERF1-3 was significantly higher expressed in H ( Figure  5d; Supplemental Figure 23C, D). As important ethylene receptors, ETR gene family members showed similar expression models with the BnACO1, BnACO3, and BnACO4 genes. CTR1 is a negative regulator and became inactivated after it activated the downstream components including EIN2, EIN3, and ERF (Shakeel et al. 2013), thus the reduced expression of BnEIN2-1 and BnEIN2-2 in H may be resulted from the high expression of BnCTR-3. However, our data suggested that some of the EIN3 and ERF gene family members may present different regulatory mechanisms in H, as the transcripts of BnEIN3-2,  in this study accumulated rapidly at T3, which may activate the downstream genes to regulate ramie fiber fineness. Overall, our data suggested that transcripts of ethylene biosynthetic genes specifically accumulated at the elongation and thickening stage of ramie fiber, and then stimulated ethylene signaling-related genes to trigger downstream genes, which may finally regulate ramie fiber fineness.

Conclusions
In conclusion, we presented a high-quality genome assembly of ramie at chromosome level, revealed the molecular mechanism of fiber fineness regulation at transcriptome level, and found that ethylene plays a critical role in regulating of ramie fiber fineness. Our study will facilitate the molecular biology research in ramie and provide valuable reference for the improvement of high-quality fiber varieties as well.

Highlights
• This study presented a high-quality genome assembly of ramie at chromosome level, which was more complete in gene content than that of the previous assembly. • This study identified linked QTLs simultaneous controlling fiber quality and fiber yield in ramie, which explain why ramie fiber fineness is inversely proportional to fiber yield. • This study revealed the molecular mechanism of fiber fineness regulation at transcriptome level for the first time, and confirmed that ethylene plays a critical role in regulating of ramie fiber fineness.