Search
Submit Manuscript
2.8CiteScore
Advanced Search
MAP  >  Emergency Management Science and Technology
2023 Volume 3
Article Contents
Abstract
Introduction
Materials and methods
Sample collection and genome sequencing
Genome survey
Genome assembly and annotation
Phylogenomic analysis and whole-genome duplication (WGD) identification
Evolution of BIA biosynthesis-related gene family
Results
Genome assembly and annotation
Genome evolution
Chromosomal location and gene duplication of BIAs biosynthesis genes in Corydalis
Phylogeny of genes involved in BIAs biosynthesis
Discussion
Extremely complex tetraploid genome of Corydalis sheareri
Tandem duplications drive the diversity of BIA biosynthetic genes in Corydalis
Conclusions
Author contributions
Data availability
Acknowledgments
Conflict of interest
Rights and permissions
References
About this article
Next Previous
ARTICLE   Open Access    

Effect of stabilizer EDTA on the thermal hazard of green synthesis process of adipic acid and development of microchannel continuous flow process

  • 1.

    Jiangsu Key Laboratory of Hazardous Chemicals Safety and Control, College of Safety Science and Engineering, Nanjing Tech University, Nanjing 210009, China

  • 2.

    BASF Corporation, 1609 Biddle Avenue, Wyandotte, MI 48192, USA

More Information
  • In the green synthesis reaction of adipic acid, the oxidant H2O2 undergoes exothermic decomposition easily and the high exothermic amount of the reaction can easily lead to thermal runaway. This study carries out calorimetric experiments on the green synthesis reaction of adipic acid, studies the effect of the stabilizer ethylenediaminetetraacetic acid (EDTA) on the safety of the synthesis reaction, and further explores the effect of EDTA on the Na2WO4-catalyzed decomposition reaction of H2O2 under the conditions of continuous flow, and develops the microchannel continuous flow process of green synthesis of adipic acid containing EDTA, in order to provide theoretical support for the improvement of intrinsic safety of green synthesis of adipic acid.

  • Introduction

    Corydalis DC. is the largest and most diverse genus of Papaveraceae, comprising more than 500 species, which are mainly divided into four subgenera (subg. Bipapillatae, subg. Cremnocapnos, subg. Sophorocapnos, and subg. Corydalis) and 39 sections[1]. Corydalis is widely distributed in the Northern Hemisphere with a few species extending into East Africa[2]. Moreover, this genus demonstrates remarkable diversity in China, particularly in the Himalaya–Hengduan Mountains and the adjacent regions[2,3]. As a genus with spectacular radiation, Corydalis exhibits an extremely high level of diversity in its leaves, subterranean organs, fruits, seeds, flower color, and the length of spurs, and can adapt to diverse habitats, such as riversides, forests, shrubs, grasslands, screes, or even cliffs. Previous studies speculated that Corydalis had an ancient origin, and underwent rapid radiation after the middle Miocene, which was most likely promoted by the continuous orogenesis and climate change associated with the uplift of the Qinghai-Tibetan Plateau (QTP)[4−7].

    The sequencing of plant genomes has greatly advanced our understanding of the underlying diversification mechanisms[8−10]. Although considerable progress has been made recently in the classification and phylogeny of Corydalis, the underlying diversification mechanisms remain poorly understood. Particularly, polyploidy is common in this genus, suggesting that polyploidization possibly contributed to the diversification of Corydalis. Polyploidy could drive dramatic changes in the genome landscapes, such as genome size and structural variation, providing a genetic basis for the adaptation and diversification of species-rich groups. While more and more studies link polyploidization events with speciation, fewer studies have documented chromosomal variation in plant groups that underwent rapid radiations, thus the extent to which it may have contributed to radiation is still elusive, especially in biodiversity hotspots. In Corydalis, only the genomes of C. tomentella and C. yanhusuo have been released to date[11,12], which has hindered our understanding of its diversification and adaptation.

    Moreover, a large number of Corydalis species have been frequently used in folk medicine due to their antibacterial, antiviral, and anticancer activities. For instance, C. yanhusuo, C. bungeana, and C. decumbens are the most famous medicinal plants recorded in the Pharmacopoeia of China (http://db.ouryao.com/yd2020/). High medicinal value of C. sheareri, C. hendersonii, C. incisa, C. repens, C. edulis, C. racemosa, and C. pallida, also have been reported[2]. Previous phytochemical investigations have isolated various components from Corydalis, including alkaloids, coumarins, flavonoids, anthraquinones, triterpenes, steroids, and organic acids[13]. Of particular interest are alkaloids, especially benzylisoquinoline alkaloids (BIAs), which have been reported to play a crucial role in sedation, releasing pain, promoting blood circulation, and inhibiting cancer cells[13]. Compared with other genera in Papaveraceae, Corydalis can produce multiple types of BIAs, such as cavidines, apocavidine, tetrahydropalmatine, corydalis, protopine, dehydroapocavidine, and dehydrocavidine. Notably, it was also reported that species-specific BIAs were also isolated in different species of Corydalis[13−16], which indicates the remarkable diversity of BIAs in these plant groups. However, the limited genomic resources hampered an in-depth understanding of the diversity of BIAs in Corydalis. Fortunately, biosynthesis pathways of BIAs have been proposed in Corydalis[12,17], as well as in other Papaveraceae species, such as Macleaya cordata[18], Eschscholzia californica[19], and Papaver somniferum[20−22]. The biosynthesis of BIAs in Corydalis involves a series of enzymes, namely berberine bridge enzyme (BBE), berberine bridge enzyme-like (BBEL), C-methyltransferase (CMT), cytochrome P450 (CYP), norcoclaurine synthase (NCS), N-methylcoclaurine 30-hydroxylase (NMCH), coclaurine N-methyltransferase (CNMT), 4-hydroxyphenylpyruvate decarboxylase (HPPDC), O-methyltransferase (OMT), tyrosine aminotransferase (TAT), tetrahydroprotoberberine N-methyltransferase (TNMT), tyrosine decarboxylase (TYDC), and tyramine 3-hydroxylase (TYR).

    In this study, PacBio long read sequencing, chromosome conformation capture (3C)-based Hi-C sequencing, and Illumina short-read sequencing were used to assemble a high-quality genome of Corydalis sheareri, one species from subg. Corydalis, the largest and most diverse lineages of Corydalis. Genomic comparison was then carried out between C. sheareri and C. tomentella, one previously reported diploid genome from subg. Sophorocapnos. Additionally, we identified the candidate BIAs biosynthesis genes in representative species in Papaveraceae and traced their evolution history by combing phylogenetic reconstruction, chromosomal location, and gene duplication analyses. Our study will provide more insights into the genome evolution as well as the BIAs diversity in Corydalis.

    Materials and methods

    Sample collection and genome sequencing

    Fresh leaves, stems, rhizomes, and flowers at different developmental stages of C. sheareri were collected from Zhongshan Botanical Garden (Nanjing, China) (Fig. 1a). After collection, these samples were immediately frozen in liquid nitrogen or dried in silica gel followed by preservation at −80 °C in the laboratory. The silica gel-dried leaves were used for the flow cytometry measurement, and the material stored in liquid nitrogen was used for genome and transcriptome sequencing. Genomic DNA was extracted using a modified CTAB method. Total RNA was extracted from leaves, stems, rhizomes, and flowers at different developmental stages using RNAprep Pure Plant Kit (Tiangen, Beijing).

    Figure 1.  Overview of Corydalis sheareri genome. (a) Photo of C. sheareri. (b) Genome survey of C. sheareri. (c) Heatmap for coverage pattern of heterozygous k-mer pairs, in which X axis indicates the normalized minor k-mer coverage, while Y axis indicates the total k-mer pairs coverage. (d) Genomic features of eight pseudochromosome. The outermost circle (blue) represents each chromosome of the genome. The bar charts of the second to fifth circles suggest gene density, LTR density, Copia, and Gypsy density, respectively. The inner circular shows inter-chromosomal synteny.
    DownLoad: Full-Size Img PowerPoint

    For the genome survey, the 350 bp paired-end library was constructed according to the Illumina protocols and sequenced on the Illumina Novaseq 6000 platform. For PacBio HiFi sequencing, the DNA SMRT bell library with an insert size of approximately 15 kb was prepared using SMRTbell® Express Template Prep Kit 2.0 (Pacific Biosciences, PN 101-853-100), and subsequently sequenced on the PacBio Sequel II platform (Pacific Biosciences, USA). For high-throughput 3C-based Hi-C sequencing, the library was generated using the standard procedures and sequenced on the Illumina Novaseq 6000 platform. For transcriptome sequencing, RNA from different tissues was pooled equally for library construction to obtain more expressed genes. Thereafter, the cDNA library with an insert size of 300−500 bp was prepared using VAHTS mRNA-seq v2 Library Prep Kit for Illumina (Vazyme) and sequenced on the Illumina Novaseq 6000 platform to generate paired-ends reads. All sequencing was carried out in Berry Genomics Company (www.berrygenomics.com), Beijing, China. The raw Illumina data for genome survey, Hi-C, and RNA-seq reads, was trimmed using Fastp v0.23.2[23] to remove the adaptors and low-quality paired reads. The HiFi reads were generated through CCS software (https://github.com/PacificBiosciences/ccs) with the following parameters: --min-passes = 3 --min-rq = 0.99.

    Genome survey

    To guide genome sequencing and assembly, flow cytometry[24], and genome survey were used to estimate the genome size. For flow cytometry, the dried material was chopped with a razor blade and then the DNA content was measured following the protocol of CyStain®PI Absolute P (Sysmex-Partec, Germany) with an Elite flow cytometer (BD FACSCalibur, USA) at the Institute of Botany, Chinese Academy of Sciences (China). Finally, the genome size was inferred based on the external reference standards (Glycine max = 1,100 Mb). For genome survey, the obtained clean reads were used to estimate the genome size based on the k-mer method with KMC v3[25], and Genomescope v2.0[26]. Given that multiple peaks were detected in the k-mer spectrum, Smudgeplot[26] was further employed to visualize and evaluate the ploidy and genome structure through the analysis of heterozygous k-mer pairs.

    Genome assembly and annotation

    Hifiasm v0.14.2[27] was used for de novo assembly of the draft contig genome of C. sheareri, which contained unphased contigs from two homolog genomes (csh v1.0). Purge_dups v1.2.3[28] was used to improve the assembly by removing duplications. The filtered Hi-C reads were aligned to the draft genome (csh v1.0) using Juicer v1.6.2[29], and scaffolding was performed on the contigs with 3D-DNA v180922[30]. Juicebox v1.9.8[31] was used to visualize the Hi-C heatmap, and the scaffolds were manually adjusted to get the chromosomal-level assembly of C. sheareri (csh v2.0). BWA v0.7.15[32] and minimap2[33] were used to assess the quality of the genome by aligning the Illumina short reads and the long HiFi reads to the genome, respectively. BUSCO v4.14[34] was performed to evaluate the integrity of the assembled genome by searching against the 1,614 conserved single-copy genes obtained from the embryophyta_odb10 database.

    LTR_Finder v1.07[35], MITE-Hunter v1.0[36], and Repeat-Masker v4.1.0 (www.repeatmasker.org) were used to predict the repeat sequences. Protein-coding genes were predicted by combining the results of ab initio-based, homology-based, and RNAseq-based predictions. For ab initio prediction, Augustus v3.2.2[37], Snap v6.0[38], Glimmer hmm v3.0.4[39], and GeneMark-ET v4.57[40] were utilized to predict the gene structure in the repeat-masked genome. GeMoMa v1.7.1[41] was used to perform homology prediction with Arabidopsis thaliana, Oryza sativa, Macleaya cordata, and Papaver rhoeas as references. For RNAseq-based prediction, transcriptomic data from different tissues were assembled de novo using trinity v2.2.0[42]. PASA r20140417[43] was used to predict the gene structure based on the obtained transcript. All predicted protein-coding genes were annotated by blast against five databases, including KEGG (www.genome.jp/kegg/brite.html), Gene Ontology (GO) terms, NR (https://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz), SwissProt (https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz), and eggNOG (http://eggnog5.embl.de/#/app/home). The tRNA (transfer RNA) was predicted by tRNAscan-SE v2.0[44]. Other types of noncoding RNAs (ncRNAs), including rRNA (ribosomal RNA), miRNA (microRNA), and snRNA (small nuclear RNA), were annotated by BLAST against the Rfam database (https://ftp.ebi.ac.uk/pub/databases/Rfam/14.1/).

    Phylogenomic analysis and whole-genome duplication (WGD) identification

    Eight species from Papaveraceae were selected for phylogenomic analysis, including C. sheareri, C. tomentella, Eschscholzia californica, Macleaya cordata, Papaver somniferum, Capnoides sempervirens, Ceratocapnos vesicaria, and Hypecoum procumbens, with Aquilegia coerulea (Ranunculaceae) as outgroup. The genomes or transcriptomes were directly retrieved from Genome Warehouse in National Genomics Data Center (https://ngdc.cncb.ac.cn/gwh), National Center for Biotechnology Information (www.ncbi.nlm.nih.gov/datasets/genome), and ONEKP database (https://db.cngb.org/onekp). Gene families were clustered using protein sequences by OrthoFinder v2.2.7[45] with the parameter '-S diamond' . Each gene family was aligned with MAFFT[46]. Single or low-copy gene families (one or two copies in the polyploid species Papaver somniferum, while one copy in the other eight species) with more than 50 amino acids were retained to reconstruct the phylogenetic tree by RAxML v8.1.17[47] with 1,000 bootstrap replicates under PROTCATWAG model using protein sequences. Divergence time was estimated under a relaxed molecular clock model by the MCMCTree implement in the PAML v4.8[48]. Calibration points retrieved from TimeTree (http://timetree.org) were used as priors in divergence time estimation, including the split of Papaveraceae and Ranunculaceae (103~117 Ma), the crown age of Papaveraceae (65~111 Ma), the split age of Hypecoideae and Fumarioideae (63~96 Ma) and the crown age of Fumarioideae (38~44 Ma). Gene family expansion and contraction were inferred by CAFE v4.1[49], with an input species tree constructed from the single or low-copy orthologs. Meanwhile, we performed functional enrichment analysis for the expanded gene families by BLAST against the KEGG and GO databases.

    The whole-genome duplication (WGD) and whole-genome triplication (WGT) events were identified by Ks method, syntenic analyses, and phylogenomic methods. For Ks method, the homologous gene pairs were firstly identified by the all-against-all BLASTP search (e-value cutoff < 1e-5). Then, YN00 in PAML v4.8[48] was called by WGDI[50] to calculate the synonymous substitution rate (Ks) of each gene pair between two species or within a single species. For syntenic analyses, collinear blocks for intra- and interspecies comparisons were detected using MCScanX v0.8[51] with '-s 15', meaning that each block contained at least 15 collinear gene pairs. JCVI v0.8.12[52] was used to draw dotplots of C. sheareri, C. tomentella, V. vinifera, and A. trichopoda with the default parameters. TBtools[53] was used to visualize the synteny between C. sheareri and C. tomentella. Additionally, structural variants between C. sheareri and C. tomentella, i.e., inversion, translocation, and duplication, were identified with SyRI v1.6.3[54]. For phylogenomic methods, gene families were firstly identified by OrthoFinder v2.2.7[45] and multiple sequence alignments were performed by MAFFT[46]. Then, Maximum-Likelihood (ML) trees were constructed using RAxML[47], with bootstrap values estimated from 100 replicates using the PROTCATWAG model. The WGD/WGT event was identified by tree2GD (https://sourceforge.net/projects/tree2gd/) with the default parameter and WGD/WGT events were considered to have occurred according to any of the following conditions: (1) gene duplication (GD) > 500, of which the number of (AB)(AB) type is over 250; (2) GD > 1,500, of which (AB)(AB) type is over 100, and at the same time, the sum of (AB)(AB) type and (AB)A or (AB)B type are over 1,000[55].

    Evolution of BIA biosynthesis-related gene family

    To gain more insights into the diversity of BIAs in Corydalis, 12 gene families involved in BIAs biosynthesis, i.e., BBE, BBEL, CMT, CYP719, CYP80B, CYP82N, NCS, NMT, OMT, TAT, TYDC, and TYR, were identified. All BIA biosynthesis genes reported in C. tomentella genome were firstly retrieved by the gene ID reported previously[12]. Then, TBLASTN was conducted to identify the candidate BIA-related genes in C. sheareri and seven other species of Papaveraceae, including Chelidonium majus, Capnoides sempervirens, Ceratocapnos vesicaria, Hypecoum procumbens, Eschscholzia californica, Macleaya cordata, Papaver somniferum, with Aquilegia coerulea as the outgroup. A sequence is regarded as a candidate gene if it encompasses the entire domain region and the pairwise amino acid identity between the queries and the targets exceed 40%. The annotation for each candidate gene was manually checked according to the blast result. All candidate sequences were confirmed by BLAST against the InterPro (www.ebi.ac.uk/interpro/) and the NCBI Conserved Domain Database (CDD, www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) with e < 1e-10.

    To infer the expansion mechanism of the BIAs biosynthesis genes, phylogenetic reconstruction, chromosomal location, and gene duplication analyses were performed. ML analyses were conducted based on the protein sequences by IQ-TREE v1.6.8[56] using the JTT model and 100 bootstrap replications. Gff files, gene files, and targeted BIA biosynthesis gene IDs of C. sheareri and C. tomentella were downloaded or extracted. The chromosome location information was obtained from the genome annotation file, and visualized by TBtools[53]. Gene duplication events were analyzed using MCScanX[51] with a BLASTp search (e < 1e-10). The synonymous (Ks) and nonsynonymous (Ka) values of the identified species-specific tandem duplicated gene pairs were calculated by the M0 model in PAML v4.8[48].

    Results

    Genome assembly and annotation

    Flow cytometry indicated that C. sheareri had a genome size of approximately 580 Mb (Supplementary Fig. S1a). Based on the obtained 66.29 Gb short paired-end reads (Supplementary Table S1), the 23-mer distribution showed that the estimated genome size of C. sheareri was 261 Mb (Fig. 1b), and the k-mer spectrum exhibited four distinct peaks at ~40, 80, 120, 160, which was highly similar to the results for autotetraploids (M. sativa and S. spontaneum)[26]. Moreover, nucleotide heterozygosity analysis showed 1.73% aaab and 1% aabb, which was consistent with the expectation that the heterozygous rate of autotetraploid AAAB would be greater than that of AABB[26]. The Smudgeplot analysis also revealed that more heterozygous k-mer pairs concentrated at 1/4 for the normalized coverage of minor k-mer and 4n for the total coverage of k-mer pairs, and the prevalence of AAAB (53%) was considerably greater than AABB (21%) (Fig. 1c). All these results suggested that the genome of C. sheareri exhibited complex genome structure and possibly an autotetraploid.

    PacBio Sequel II sequencing yielded approximately 33.02 Gb of high quality Hifi reads, and the average length and N50 of filtered subreads (1,990,989) were 16,584 bp and 16,614 bp, respectively (Supplementary Table S1). De novo assembly generated the draft genome of C. sheareri (csh v1.0, ~550 Mb) with contig N50 of 9.18 Mb, and the longest contig was approximately 25.09 Mb (Table 1). The assembly was further scaffolded with 67.96 Gb Hi-C data (Supplementary Table S1). Finally, the high-quality genome of C. sheareri (csh v2.0, ~282 Mb) comprised eight pseudochromosomes and 36 contigs (Supplementary Fig. S1b, Supplementary Table S2), with contig N50 of 11.39 Mb (Table 1). The minimum length of the chromosome was greater than 30 Mb (Supplementary Table S2). Approximately 97.62% of DNA reads and 88.32% of RNA-seq reads could be mapped to the assembly genome. BUSCO analysis revealed that 97.1% (1,567/1,614) of the core eukaryotic genes were completely present in the C. sheareri genome, of which 91.4% were single copy (1,475) and 5.7% were duplicated (92), while 0.9% (15), and 2.0% (32) were partially present or missing, respectively.

    Table 1.  Genome assembly and annotation of Corydalis sheareri.
    Analytical process Characteristic C. sheareri
    Genome survey Genome size (flow cytometry) (Mb) 580
    Genome size (k-mer spectrum) (Mb) 261
    Assembly_csh v1.0 Number of contigs 2,713
    Assembly size (Mb) 550
    Contig N50 (Mb) 9.18
    Shortest contig (bp) 16,037
    Largest contig (bp) 25,090,490
    Assembly_csh v2.0 Total number of contigs 179
    Assembly size (Mb) 282
    Contig N50 (Mb) 11.39
    Number of pseudochromosomes 8
    GC content 37.11%
    Annotation Number of protein-coding genes 26,287
    Mean gene length (bp) 5,092
    Mean CDS length (bp) 1,440
    Complete BUSCOs (C) 1,567
    Percentage of repeat sequences (%) 44.47
     | Show Table
    DownLoad: CSV

    Approximately 44.47% (125,404,725 bp) of the C. sheareri genome was annotated as transposable elements (TEs), of which 26.47% were retrotransposons and 7.37% were transposons (Supplementary Table S3). For retrotransposons, a total of 84,322 LTR elements were identified, of which 53,097 (7.27%) belonged to the Copia superfamily and 29,747 (5.05%) belonged to the Gypsy superfamily (Supplementary Table S3). Based on a combination of homology search, de novo prediction, and RNA-seq based prediction, a total of 26,287 protein-coding genes were confidently annotated for C. sheareri, and the mean lengths of the predicted gene and coding sequence were 5,092 and 1,440 bp, respectively (Table 1). Approximately, 92.39% (24,286/26,287) of the genes were functionally annotated, of which, 24,257, 19,162, 8,566, 5,531, and 23,330 genes showed high similarity to known proteins in the NR, SwissProt, KEGG, GO, and eggNOG databases, respectively (Supplementary Fig. S2). In addition, a total of 2,746 noncoding RNA genes were identified in the genome of C. sheareri, including 1,014 tRNA genes, 919 rRNA genes, 609 snRNA genes, and 86 miRNA genes (Supplementary Table S4).

    Genome evolution

    A total of 28,545 orthologous groups (OGs) were identified in nine selected species (Fig. 2b, Supplementary Table S5). Of them, 7,769 gene families were shared by all species, while 1,714 gene families were specific to Corydalis. Additionally, 1,025 single or low-copy nuclear gene families were identified in these species. After removing the OGs less than 50 aa, 1,003 single or low-copy nuclear gene families were retained to infer a high-confidence species tree. As expected, three subfamilies were recovered in the phylogenomic analysis, and Hypecoideae was strongly supported as a sister to Fumarioideae. Within Fumarioideae, C. sheareri and C. tomentella formed a highly supported clade, while the relationship of Ceratocapnos vesicaria, Capnoides sempervirens, and Corydalis was not fully resolved. Within Papaveroideae, Eschscholzia californica diverged firstly, and Papaver somniferum was sister to Macleaya cordata (Fig. 2a). Molecular dating indicated that the divergence of C. sheareri and C. tomentella was dated to 24.92 Ma, with a 95% confidence interval (95% CI) of 16.35~33.22 Ma (Fig. 2a). Gene family expansion and contraction analyses revealed that 674 and 566 gene families expanded and contracted in Corydalis, respectively (Fig. 2a). GO and KEGG enrichment analyses showed that the significantly expanded gene families were mainly related to response to stimulus, membrane, cell periphery, response to chemical, and cellular response to stimulus, which primarily enriched in secondary metabolites pathways such as phenylpropanoid biosynthesis, pentose, and glucuronate interconversions, photosynthesis, flavonoid biosynthesis, and monoterpenoid biosynthesis (Supplementary Fig. S3). Specifically, cytochrome P450 (CYP) and photosynthesis proteins were also enriched in the KEGG analysis (Supplementary Fig. S3).

    Figure 2.  Genome evolution of Corydalis sheareri. (a) Phylogenetic tree and divergence time estimation. Pal., Palaeogene; Neo., Neogene; Hyp., Hypecooideae. (b) Gene families identified in each species. (c) Collinearity within the genome of C. sheareri. (d) Ks distributions of anchor pairs for the paralogous genes of C. sheareri, C. tomentella and Aquilegia coerulea, and for orthologous genes between C. sheareri and C. tomentella, A. coerulea, respectively. (e) Collinearity analysis between C. sheareri and C. tomentella. (f) Structural variation detection between C. sheareri and C. tomentella performed by SyRI.
    DownLoad: Full-Size Img PowerPoint

    WGD events play a crucial role in the duplication and retention of genes. Intra-genomic colinearity analyses uncovered remnants of one WGD event in C. sheareri (Fig. 2c). Obviously, one signature peak of synonymous substitutions per synonymous site (Ks) distribution was detected for the C. sheareri genome at approximately 1.0 (Fig. 2d), indicating an ancient WGD event. Similarly, previously sequenced genomes of Ranunculales species, including C. tomentella, and A. coerulea, also showed a signature peak at 1.0~1.2 in their genomes (Fig. 2d), suggesting this WGD event was probably shared by all Ranunculales species. Comparison of the C. sheareri paralogue Ks distribution between Corydalis and Aquilegia coerulea also indicated a WGD occurred in Ranunculales (Fig. 2d). Previous studies reported that the Vitis vinifera genome had an ancestral hexaploidization[57] and the Amborella trichopoda genome shows no evidence of more recent lineage-specific genome duplications[58]. The inter-genomic comparison and syntenic depths analyses among Corydalis sheareri, Vitis vinifera, and Amborella trichopoda further confirmed that this ancient WGD event, in which two paralogous segments in the Corydalis sheareri genome correspond to one and three orthologous regions in the Amborella trichopoda (with no WGD), and Vitis vinifera (with one WGT) genomes, respectively (Supplementary Fig. S4). The tree2GD analyses revealed 1,281 GDs with 934 (73%) (AB)(AB) retention type, which also strongly suggested the identified WGD event was shared by the ancestor of Ranunculales (node 5) (Supplementary Fig. S5). The genome of Corydalis sheareri and C. tomentella exhibited good collinearity (Fig. 2e, f, Supplementary Fig. S6). A number of large structural variations were detected, including multiple invasions across the genome, and a small fraction of translocations on chromosomes 6 (Fig. 2f). To further verify the authenticity of these structural variations, several specific inversion regions were randomly checked by mapping the HiFi reads to the assembled genome with Integrative Genomics Viewer (IGV)[59]. Our results confirmed that these inversions are genuine and not assembly artifacts (Supplementary Fig. S7).

    Chromosomal location and gene duplication of BIAs biosynthesis genes in Corydalis

    A total of 172 genes involved in the BIAs biosynthesis have been identified in Corydalis sheareri, which is notably larger than that found in C. tomentella (125), Macleaya cordata (128), and Aquilegia coerulea (138), while smaller than that observed in Papaver somniferum (364) and Eschscholzia californica (197) (Fig. 3b; Supplementary Table S6), which has undergone two WGD events[12]. For each gene family, the upstream BIA biosynthetic genes, such as TAT, TYDC, TYR, BBE, and CYP80B (NMCH), generally possess small copy numbers in Corydalis. Conversely, the downstream BIA biosynthetic genes, including OMT, CMT, NMT, and BBEL, typically have a higher copy number (Fig. 3b; Supplementary Table S6). The only exception is NCS, which encodes enzymes to catalyze the condensation of dopamine and 4-hydroxyphenylacetaldehyde (4-HPAA) to produce norcoclaurine, in the upstream of BIA biosynthesis, which also has a greater number of copies (Fig. 3a, b).

    Figure 3.  BIAs biosynthesis related gene families identified in this study. (a) The biosynthetic pathway of benzylisoquinoline alkaloids (BIAs) in Corydalis. The species-specific tandem duplicated gene pairs or clusters in Corydalis identified in Fig. 5 and Table 2 were highlighted with red color. (b) The copy number of BIAs biosynthesis related gene families. Four species with transcriptome sequence are represented by asterisks. (c) Gene duplication type identified in C. sheareri (left) and C. tomentella (right). Abbreviations: BBE, berberine bridge enzyme; BBEL, berberine bridge enzyme-like; CMT, C-methyltransferase; CNMT, coclaurine N-methyltransferase; CYP, Cytochrome P450; HPPDC, 4-hydroxyphenylpyruvate decarboxylase; NCS, norcoclaurine synthase; NMT, N-methyltransferase; OMT, O-methyltransferase; TAT, tyrosine aminotransferase; TNMT, tetrahydroprotoberberine N-methyltransferase; TYDC, tyrosine decarboxylase; TYR, tyramine 3-hydroxylase. Aco, Aquilegia coerulea; Cma, Chelidonium majus; Cse, Capnoides sempervirens; Csh, Corydalis sheareri; Cto, C. tomentella; Cve, Ceratocapnos vesicaria; Eca, Eschscholzia californica; Hpr, Hypecoum procumbens; Mco, Macleaya cordata; Pso, Papaver somniferum. SD, Singleton; DD, dispersed duplications; PD, proximal duplications; TD, Tandem duplications; WGD, segmental/WGD duplications.
    DownLoad: Full-Size Img PowerPoint

    Chromosomal location analysis indicated that some BIA biosynthetic genes are unevenly distributed among chromosomes. For instance, all 22 NCS genes of C. sheareri are located on chromosome 6, whereas, in C. tomentella, all members are mapped on the chromosome 5. Additionally, in C. sheareri, 34 out of the 39 BBEL genes are located on chromosome 4, 17 out of the 22 NMT genes are mapped on the chromosome 6, and nine out of 15 CYP82N genes are found on chromosome 3 (Fig. 4a). Similarly, in C. tomentella, 25 out of 31 BBEL genes are found on chromosome 6, and 15 out of the 20 NMT genes are located on chromosome 5 (Fig. 4b). On the contrary, some other genes are widely distributed throughout genomes but are uneven among chromosomes in both two species. For instance, seven out of eight chromosomes (except for chr02) harbor the 22 CMT genes in C. sheareri, and 11 of them are located on chromosome 4. Furthermore, five chromosomes (1, 2, 4, 6, and 8) contain the OMT genes, and both chromosome 2 and 4 harbor nine members (Fig. 4a). In C. tomentella, 20 OMT genes are distributed across five chromosomes (1, 2, 4, 5, and 6). and chromosome 6 harbors eight members, chromosome 1 harbors four (Fig. 4b). Nineteen CMT genes are located on seven chromosomes (except for chr01), and 9 of them are located on chromosome 8 (Fig. 4b). Interestingly, a gene cluster including seven NCS genes and ten NMT genes, was identified in chromosome 6 within a 270-kb region of C. sheareri genome (Fig. 4a).

    Figure 4.  Chromosome location and duplication event analyses of BIA biosynthesis genes in (a) C. sheareri and (b) C. tomentella. The chromosome number is indicated at the top of each chromosome. Tandem duplicated genes are indicated with orange color. The scale bar on the left indicates the length (mb) of chromosomes. The gene and species abbreviations were the same as Fig. 3. Chr, chromosome.
    DownLoad: Full-Size Img PowerPoint

    Gene duplication analysis revealed that the majority of BIAs biosynthesis genes were generated through gene duplication events in Corydalis. Over half of these genes were identified as tandem duplications, with 55.2% in C. tomentella and 64.5% in C. sheareri, respectively (Fig. 3c; Supplementary Table S6). Additionally, 18% to 27.2% of genes were identified as dispersed duplications, and one-tenth of genes, corresponding to 10.5% to 12.0%, were recognized as proximal duplications (Fig. 3c; Supplementary Table S6). Meanwhile, a small fraction, ranging from 4.7% to 5.6%, were identified as segmental/WGD (Fig. 3c; Supplementary Table S6). Specifically, BBEL genes exhibit an extremely high probability of tandem duplication in both C. sheareri and C. tomentella. For instance, 31 BBEL genes (79.4%) of C. sheareri were found to be tandemly duplicated on chromosome 4, while 24 members (80%) of C. tomentella tandemly duplicated on chromosome 6 (Figs 3c, 4; Supplementary Table S6). Moreover, NCS, CYP719, and TYR genes also demonstrated an extremely high rate of tandem duplication in C. sheareri, with 19 out of 22 NCS genes (86%), five out of six CYP719 genes (83%) and nine out of 11 TYR genes (82%) were observed to be tandemly duplicated (Figs 3c, 4a; Supplementary Table S6).

    Phylogeny of genes involved in BIAs biosynthesis

    Among the genes associated with the BIAs biosynthesis, the largest is the BBEL family, with 39 members in Corydalis sheareri and 30 members in C. tomentella, respectively (Fig. 3b; Supplementary Table S6). Phylogenetic analysis of BBE and BBEL revealed three monophyletic clades (BBE, I, and II). Clade II has undergone a significant expansion in Corydalis, encompassing 34 members from C. sheareri and 23 from C. tomentella. Seven C. sheareri or C. tomentella-specific monophyletic groups were identified in clade II, and one of them comprised a notably high number of 14 members from C. sheareri (Fig. 5a).

    Figure 5.  Phylogenetic trees of BIAs biosynthesis genes. (a) BBE and BBEL, (b) CYP719, CYP80B (NMCH), and CYP82N, (c) NCS, (d) CMT, (e) NMT, (f) OMT, (g) TAT, (h) TYDC, and (i) TYR. Corydalis sheareri and C. tomentella – specific monophyletic groups are highlighted with green and pink color, respectively. The gene and species abbreviations were the same as Fig. 3.
    DownLoad: Full-Size Img PowerPoint

    CYPs are important for determining chemical diversity in metabolism, in which, CYP719, CYP80B (NMCH), and CYP82N, have been identified as key components of BIA biosynthesis[12,60,61]. Phylogenetic analysis of these three CYP subfamily suggested that 13 CYP82N genes of C. sheareri (a total of 15 members) were clustered into four monophyletic groups, and two members of C. tomentella (CtoCYP82N5 and CtoCYP82N6) were clustered into one highly supported monophyletic group (Fig. 5b). For CYP719 genes, two C. sheareri specific monophyletic groups were identified, consisting of three (CshCYP719_01, CshCYP719_02, and CshCYP719_03) and two members (CshCYP719_04 and CshCYP719_06), respectively. Similarly, four members of C. tomentella clustered into two monophyletic groups in the phylogenetic tree (Fig. 5b).

    Phylogenetic analysis of the NCS genes revealed four highly supported monophyletic clades (I−IV). In clade III, two C. sheareri specific monophyletic groups were identified, comprising seven (CshNCS02, CshNCS03, CshNCS04, CshNCS05, CshNCS06, CshNCS16 and CshNCS17), and three (CshNCS19, CshNCS20 and CshNCS21) members, respectively (Fig. 5c). In clade IV, one C. sheareri specific monophyletic group, and one C. tomentella specific monophyletic group were identified, containing three (CshNCS07, CshNCS08 and CshNCS09) and two members (CtoNCS4 and CtoNCS9), respectively (Fig. 5c).

    Phylogeny the CMT genes revealed seven monophyletic clades (I−VII), and the majority of Corydalis members are scattered throughout the phylogenetic tree (Fig. 5d). Conversely, phylogenetic tree of the NMT genes revealed five highly supported clades (I−V) and eight species-specific monophyletic groups were recognized in Corydalis (Fig. 5e). In clade II, three species-specific monophyletic groups were identified, comprising four C. sheareri members (CshNMT04, CshNMT05, CshNMT07, and CshNMT08), four C. tomentella members (CtoNMT14, CtoNMT15, CtoNMT16, and CtoNMT17), and two C. sheareri members (CshNMT19 and CshNMT20), respectively. In clade IV, four species-specific monophyletic groups, with three C. sheareri members (CshNMT17, CshNMT21, and CshNMT22), three C. tomentella members (CtoNMT4, CtoNMT5, and CtoNMT6), two C. sheareri members (CshNMT14 and CshNMT16) and two C. tomentella members (CtoNMT2 and CtoNMT3), respectively. In clade V, one monophyletic group, containing five C. sheareri members (CshNMT10, CshNMT11, CshNMT12, CshNMT13, and CshNMT15), was identified. Similarly, phylogenetic analysis of the NMT genes revealed eight well-supported clades (I−VIII). In clade III, three C. tomentella members (Cto9OMT, CtoOMT10, CtoOMT11) formed one well-supported monophyletic group. In clade VIII, another three C. tomentella members (Cto6OMT, CtoOMT7, and CtoOMT8) also formed one well-supported monophyletic group. While in clade V, two C. sheareri - specific monophyletic groups were identified, containing four (CshOMT03, CshOMT04, CshOMT05, and CshOMT09), and three (CshOMT08, CshOMT10, and CshOMT11) members, respectively.

    Notably, the close relationship of the species-specific tandem duplicated gene pairs or clusters with high sequence similarity was confirmed in the phylogenetic analyses, such as CshBBEL03 and its paralogs (CshBBEL04 and CshBBEL05), CshBBEL30 and its paralog CshBBEL31, CshCYP719_01 and its paralogs (CshCYP719_02 and CshCYP719_03), CshNCS07 and its paralogs (CshNCS08 and CshNCS09), CshNCS19 and its paralogs (CshNCS20 and CshNCS21), CshNMT19 and its paralog CshNMT20, CshTYR03 and its paralogs (CshTYR04, CshTYR05, and CshTYR06), CtoCYP82N5 and its paralog CtoCYP82N6, and CtoCYP719A2 and its paralog CtoCYP719A2-2 (Fig. 5). The Ka/Ks ratio of these gene pairs ranged from 0.15543 to 0.43823 with an average of 0.292909 (Table 2), suggesting that purifying selection was the primary evolutionary force on these species-specific tandem duplicated gene pairs or clusters.

    Table 2.  Selective pressure and sequence similarity of the species-specific tandem duplicated gene pairs or clusters in Corydalis.
    Gene pairs or clusters Ka Ks Ka/Ks Similarity
    CshBBEL03, CshBBEL04, CshBBEL05 0.1656 0.4321 0.38315 97.54%
    CshBBEL30, CshBBEL31 0.1090 0.3349 0.32543 98.73%
    CshCYP719_01, CshCYP719_02,
    CshCYP719_03    		 0.0841 0.3688 0.22803 98.81%
    CshNCS07, CshNCS08, CshNCS09 0.0546 0.2236 0.24430 99.11%
    CshNCS19, CshNCS20, CshNCS21 0.1325 0.8527 0.15543 94.61%
    CshNMT19, CshNMT20 0.1435 0.4799 0.29900 98.68%
    CshTYR03, CshTYR04, CshTYR05, CshTYR06 0.1635 0.373 0.43823 97.79%
    CtCYP719A2, CtCYP719A2-2 0.0841 0.3808 0.22096 97.73%
    CtCYP82N5, CtCYP82N6 0.0913 0.2672 0.34165 97.79%
     | Show Table
    DownLoad: CSV

    Discussion

    Extremely complex tetraploid genome of Corydalis sheareri

    In this study, by combining data from PacBio long-read sequencing, 3C-based Hi-C sequencing, and Illumina short-read sequencing, we assembled the genome of Corydalis sheareri, one species from subg. Corydalis, the largest and most diverse lineages of Corydalis. Genome survey showed that its genome is extremely complex (Fig. 1b). Both the GenomeScope and Smudgeplot analyses implied that the genome structure of C. sheareri might be an autotetraploid (a special tetraploid with three homologous chromosomes and one non-homologous chromosome) with a special AAAB karyotype (Fig. 1b, c). Intriguingly, the karyotype of C. sheareri is remarkably similar to that of C. yanhusuo[11], which also belongs to subg. Corydalis.

    As previously reported, the diploid C. tomentella has a genome size of 258 Mb[12], and the estimated tetraploid C. sheareri genome size was 580 Mb (Supplementary Fig. S1a), which is nearly twice that of C. tomentella. This indicates that polyploidization might have played a significant role in the genome evolution within this genus. It is noteworthy that both C. sheareri and C. yanhusuo are tetraploid, yet their genome sizes differ by more than threefold. In Corydalis, both the genome size (https://cvalues.science.kew.org/search/angiosperm) and the chromosome number (https://ccdb.tau.ac.il/; http://legacy.tropicos.org/Project/IPCN) varied considerably. We deduce that diploidization following polyploidization, chromosomal rearrangements, including inversions, translocations, and changes in chromosome number via fusion and fission, as well as gene loss might be common and could thus trigger the genome size diversity in Corydalis. In our study, the detection of large-scale chromosomal structural variants, especially multiple inversions, between the genomes of C. sheareri and C. tomentella (Fig. 2e, f), seems to provide strong evidence in support of this hypothesis.

    Tandem duplications drive the diversity of BIA biosynthetic genes in Corydalis

    Polyploidization and structural variations might not only lead to changes in genome size but also substantially affect the gene content. In this study, we identified 172 candidate genes involved in the BIAs biosynthesis within the C. sheareri genome and traced their evolutionary history in Papaveraceae. As previously reported, BBELs, CMTs, NMTs, OMTs, and CYPs (CYP719, CYP82N) are the key enzymes downstream of the BIA biosynthetic pathway (Fig. 3a). In the Corydalis genome, the genes that encode these downstream enzymes, such as BBEL, CMT, NMT, and OMT, tend to possess a larger number of copies than the upstream genes of the BIA biosynthetic pathway (Fig. 3a, b), which enables plants to synthesize various BIAs. Interestingly, NCS, one upstream gene family, also has a relatively high copy number (Fig. 3b). NCS catalyzes the condensation of dopamine and 4-HPAA to produce norcoclaurine, which was identified as one of rate - limiting enzymes in BIA biosynthesis in opium poppy[62].

    Furthermore, the chromosomal location and gene duplication analyses demonstrated that these genes frequently appeared as tandem duplications scattered throughout the genome (Figs 3c, 4). All these findings corroborated that tandem duplications might play a key role in the diversity of BIA biosynthetic genes in Corydalis. What's more, our comprehensive phylogenetic analyses, which cover representative species throughout Papaveraceae, revealed that clade I and II of NCS, clades I, II, and VI of CMT, clades I and III of NMT, clades II, VI, and VII of OMT, two members of TAT, four members of TYDC, and clade I of TYR, each retain a single copy in Corydalis (Fig. 5). In addition, the fact that their homologs are shared by the majority of species suggests that they likely originate before the divergence of Papaveraceae or even Ranunculales. By contrast, the identification of C. sheareri or C. tomentella - specific monophyletic groups in BBEL, CYP719, CYP82N, clades III and IV of NCS, clades V and VI of CMT, clades II, IV and V of NMT, clades V and VIII of OMT, clade III of TYR (Fig. 5) indicate that recent gene duplications might occur frequently for these members in Corydalis. Particularly, the close relationship of some tandem duplicated gene pairs or clusters was confirmed in our phylogenetic analyses (Fig. 5), strongly suggesting that these genes may be generated from the recent duplication events that occurred within C. sheareri, or C. tomentella. Additionally, more C. sheareri-specific tandem duplication events are identified compared to those in C. tomentella, which is likely associated with the significant expansion of BIAs biosynthetic genes in C. sheareri.

    Until recently, with the burst of plant genome sequencing projects and the advancement of bioinformatic tools, biosynthetic pathways for many natural products have been elucidated[63−65]. Numerous studies have stated that tandem duplication is a major factor contributing to the diversity of secondary metabolites biosynthesis, by recruiting novel genes and potentially introducing new metabolic pathways. For instance, Papaver somniferum has undergone significant tandem duplication events, which result in the emergence of morphinan and noscapine biosynthesis pathways[22,66]. The divergence and expansion of CYP genes strongly contribute to the alkaloid diversity in Coptis[67]. Tandem duplications are also common for the triterpene biosynthetic genes in Aralia elata, especially for CYP72A, CSLM, and UGT73, which may drive the diversity of triterpenoids[68]. In the case of Scutellaria, tandem duplications of the CYP82D subfamily shape the flavonoid diversification[69]. In Corydalis, each species was found to contain a particular set of BIAs, some of which are common to other species but not in the same combinations. BBEL genes have been reported to be expanded and considered to be related to the appearance of cavidines and coptisine in C. tomentella[12]. In this study, tandem duplications are also found to be evident in other genes involved in the biosynthesis of BIAs, particularly NCS, CMT, NMT, OMT, and CYP. Despite not conducting the functional experiment, we still have reason to believe that tandem duplications may drive the expansion of BIA biosynthesis genes, which is conducive to the complexity of biosynthesis pathway and further contributes the BIAs diversity in Corydalis. The incorporation of sequencing data from more Corydalis species, in conjunction with the conduction of multi-omics landscape investigations and the performance of functional experiment studies related to these recently duplicated genes in the future, may lead to an updated canvas for illustrating the genetic mechanisms underlying the diversity of BIAs.

    Conclusions

    In this study, we proposed that C. sheareri might be a complex autotetraploid with a special AAAB karyotype. Genomic comparison detected large syntenic blocks between C. sheareri and its relatives C. tomentella, and also uncovered large-scale chromosomal structural variations (particularly inversions) between these two genomes, which might have profound effects on the divergence of Corydalis. Furthermore, we identified 172 candidate genes involved in BIAs biosynthesis in C. sheareri and traced their evolution history in Papaveraceae. We deduce that tandem duplication has played a prominent role in the expansion of the BIAs biosynthesis genes, especially for BBEL, CMT, NMT, OMT, and NCS. Moreover, the identification of the species-specific tandem duplication events implies that the growth in the number of gene members is likely to complicate the metabolic pathway, and consequently, has further contributed to the diversification of BIAs biosynthesis, ultimately leading to the diversity of BIAs in Corydalis. Our study provides more insights into the genome evolution for species-rich taxa with radiation, as well as the mechanisms underlying the diversity of the BIA biosynthetic pathway. It is also of great value for future genetic studies and medicinal applications of Corydalis.

  • The authors confirm contribution to the paper as follows: study conception and design: Liu YY, Peng D, Li JM; data collection: Liu YY, Peng D, Yu CL, Liu YJ, Chen M, Kan SL, Cao YN; analysis and interpretation of results: Liu YY, Peng D, Yu CL, Liu YJ, Kan SL, Cao YN; draft manuscript preparation: Liu YY, Peng D, Wang HW, Li JM, Peng D. All authors reviewed the results and approved the final version of the manuscript.

  • The whole-genome sequence data, including Illumina short reads, PacBio HiFi reads, Hi-C interaction reads, transcriptome data, and genome annotation files, have been deposited in The National Genomics Data Center (NGDC), under the project number: PRJCA035358.

    • Acknowledgments

    We thank Dr Yuan-Yuan Feng (Institute of Botany, Chinese Academy of Sciences) for her help in the flow cytometry measurement. We also thank Mr Hai-Kuan Zhang (Berry Genomics Company) for his help in the genome assembly and annotation. This research was funded by National Natural Science Foundation of China (32000170) and Xinyang Academy of Ecological Research Open Foundation (2023XYMS05).

  • The authors declare that they have no conflict of interest.

  • [1]

    Noyori R, Aoki M, Sato K. 2003. Green oxidation with aqueous hydrogen peroxide. Chemical Communications 2023(16):1977−86

    doi: 10.1039/b303160h

    CrossRef   Google Scholar

    [2]

    Van de Vyver S, Román-Leshkov Y. 2013. Emerging catalytic processes for the production of adipic acid. Catalysis Science & Technology 3(6):1465−79

    doi: 10.1039/C3CY20728E

    CrossRef   Google Scholar

    [3]

    Jin P, Zhao Z, Dai Z, Wei D, Tang M, et al. 2011. Influence of reaction conditions on product disTribution in the green oxidation of cyclohexene to adipic acid with hydrogen peroxide. Catalysis Today 175(1):619−624

    doi: 10.1016/j.cattod.2011.04.041

    CrossRef   Google Scholar

    [4]

    Chen KY, Lin CM, Shu CM, Kao CS. 2006. An evaluation on thermokinetic parameters for hydrogen peroxide at various concentrations by DSC. Journal of Thermal Analysis and Calorimetry 85(1):87−89

    doi: 10.1007/s10973-005-7362-8

    CrossRef   Google Scholar

    [5]

    Eto I, Akiyoshi M, Miyake A, Ogawa T, Matsunaga T. 2009. Hazard evaluation of runaway reaction of hydrogen peroxide – Influence of contamination of various ions. Journal of Loss Prevention in the Process Industries 22(1):15−20

    doi: 10.1016/j.jlp.2008.09.009

    CrossRef   Google Scholar

    [6]

    Wu L, Chen K, Cheng S, Lee B, Shu C. 2008. Thermal decomposition of hydrogen peroxide in the presence of sulfuric acid. Journal of Thermal Analysis and Calorimetry 93(1):115−120

    doi: 10.1007/s10973-007-8829-6

    CrossRef   Google Scholar

    [7]

    Zhao H. 2003. Research progress on H2O2 stabilizer. Science & Technology in Chemical Industry 11(2):55−59

    doi: 10.16664/j.cnki.issn1008-0511.2003.02.015

    CrossRef   Google Scholar

    [8]

    Zhang X, Wang L, Li Y, Lei L. 2016. Decomposition of hydrogen peroxide by manganese in the presence of lignin. Appita Journal 69(3):241−46

    Google Scholar

    [9]

    Jung YS, Park JY, Ko SO, Kim YH. 2013. Stabilization of hydrogen peroxide using phthalic acids in the Fenton and Fenton-like oxidation. Chemosphere 90(2):812−819

    doi: 10.1016/j.chemosphere.2012.09.089

    CrossRef   Google Scholar

    [10]

    Li X. 2019. Study on process safety and reaction hazard in peracetic acid production. Thesis. Beijing Institute of Petrochemical Technology, Beijing.
    [11]

    Liu S. 2021. Study on thermal hazard of synthesis process of N-methyl morpholine-N-oxide based on calorimetry. Thesis. Beijing Institute of Petrochemical Technology, Beijing.
    [12]

    Parker T, Mao Y, Wang Q. 2022. Application of response surface methodology for hazard analysis of 2-butanol oxidation to 2-butanone using RC1 calorimetry. Journal of Loss Prevention in the Process IndusTries 75:104703

    doi: 10.1016/j.jlp.2021.104703

    CrossRef   Google Scholar

    [13]

    Zhao H, Liu S, Shang M, Su Y. 2021. Direct oxidation of benzene to phenol in a microreactor: Process parameters and reaction kinetics study. Chemical Engineering Science 246:116907

    doi: 10.1016/j.ces.2021.116907

    CrossRef   Google Scholar

    [14]

    Constantinou A, Wu G, Corredera A, Ellis P, Bethell D, et al. 2015. Continuous heterogeneously catalyzed oxidation of benzyl alcohol in a ceramic membrane packed-bed reactor. Organic Process Research & Development 19(12):1973−79

    doi: 10.1021/acs.oprd.5b00220

    CrossRef   Google Scholar

    [15]

    Jiang W, Ni L, Jiang J, Chen Q, Chen Z, et al. 2021. Thermal hazard and reaction mechanism of the preparation of adipic acid through the oxidation with hydrogen peroxide. AIChE Journal 67:e17089

    doi: 10.1002/aic.17089

    CrossRef   Google Scholar

    [16]

    Zang N, Qian XM, Liu ZY, Shu CM. 2016. Themal hazard evaluation of cyclohexanone peroxide synthesis. Journal of Thermal Analysis and Calorimetry 124:1131−39

    doi: 10.1007/s10973-015-5209-5

    CrossRef   Google Scholar

    [17]

    Stoessel F. 2008. General Aspects of Reactor Safety. Thermal Safety of Chemical Processes: Risk Assessment and Process Design. New York: Wiley-VCH Verlag GmbH & Co. KGaA. pp. 101−17. https://doi.org/10.1002/9783527621606.ch5
    [18]

    Liu Y, Ni L, Yao H, Su J, Cheng Z. 2023. The effects of acid and hydrogen peroxide stabilizer on the thermal hazard of adipic acid green synthesis. Sustainability 15:6530

    doi: 10.3390/su15086530

    CrossRef   Google Scholar

    [19]

    Shang M, Noël T, Su Y, Hessel V. 2017. Kinetic study of hydrogen peroxide decomposition at high temperatures and concentrations in two capillary microreactors. AIChE Journal 63:689−97

    doi: 10.1002/aic.15385

    CrossRef   Google Scholar

    [20]

    Kuhn S, Hartman RL, Sultana M, Nagy KD, Marre S, et al. 2011. Teflon-coated silicon microreactors: impact on segmented liquid-liquid multiphase flows. Langmuir 27:6519−27

    doi: 10.1021/la2004744

    CrossRef   Google Scholar

    [21]

    Santacesaria E, Di Serio M, Tesser R, Turco R, Tortorelli M, et al. 2012. Biodiesel process intensification in a very simple microchannel device. Chemical Engineering and Processing: Process Intensification 52:47−54

    doi: 10.1016/j.cep.2011.12.001

    CrossRef   Google Scholar

    [22]

    Ross J. 1986. Ullman's encyclopedia of industrial chemistry. Applied Catalysis 27(2):403−4

    doi: 10.1016/S0166-9834(00)82943-7

    CrossRef   Google Scholar

  • Cite this article

    He W, Li Y, Ni L, Zhu W. 2023. Effect of stabilizer EDTA on the thermal hazard of green synthesis process of adipic acid and development of microchannel continuous flow process. Emergency Management Science and Technology 3:22 doi: 10.48130/EMST-2023-0022
    He W, Li Y, Ni L, Zhu W. 2023. Effect of stabilizer EDTA on the thermal hazard of green synthesis process of adipic acid and development of microchannel continuous flow process. Emergency Management Science and Technology 3:22 doi: 10.48130/EMST-2023-0022
    shu

Figures(10)  /  Tables(7)

Download PDF

Article Metrics

Article views(3919) PDF downloads(493)

Other Articles By Authors

ARTICLE   Open Access    

Effect of stabilizer EDTA on the thermal hazard of green synthesis process of adipic acid and development of microchannel continuous flow process

  • 1.

    Jiangsu Key Laboratory of Hazardous Chemicals Safety and Control, College of Safety Science and Engineering, Nanjing Tech University, Nanjing 210009, China

  • 2.

    BASF Corporation, 1609 Biddle Avenue, Wyandotte, MI 48192, USA

  • Received: 23 November 2023
  • Accepted: 16 December 2023
  • Published online: 27 December 2023
Emergency Management Science and Technology  3,  Article number: 22  (2023)  |  Cite this article

Abstract: In the green synthesis reaction of adipic acid, the oxidant H2O2 undergoes exothermic decomposition easily and the high exothermic amount of the reaction can easily lead to thermal runaway. This study carries out calorimetric experiments on the green synthesis reaction of adipic acid, studies the effect of the stabilizer ethylenediaminetetraacetic acid (EDTA) on the safety of the synthesis reaction, and further explores the effect of EDTA on the Na2WO4-catalyzed decomposition reaction of H2O2 under the conditions of continuous flow, and develops the microchannel continuous flow process of green synthesis of adipic acid containing EDTA, in order to provide theoretical support for the improvement of intrinsic safety of green synthesis of adipic acid.

    HTML

    Introduction

    • Adipic acid is a basic raw material to produce chemical products such as chemical fiber, nylon 66 and engineering plastics, and is widely used in the pharmaceutical, food and chemical industries. The synthesis of adipic acid by oxidation of cyclohexene with hydrogen peroxide (H2O2) catalyzed by sodium tungstate (Na2WO4) is a common green synthesis process of adipic acid, as shown in Fig. 1[1−3]. Since H2O2 is easy to decompose and exothermic and the high exothermic amount of this reaction can easily lead to thermal runaway of the reaction[4], improvements in the safety of the process is of great significance for the popularization and application of this green process.

      Figure 1. 

      Reaction equation for the synthesis of adipic acid by oxidation of cyclohexene with H2O2 catalyzed by Na2WO4.

      In H2O2, unstable -O-O- exists, which can easily decompose to produce O2 under certain conditions, and many researchers have conducted studies related to the stability and thermal hazards of H2O2. Eto et al.[5] studied and evaluated the influence law of various metal compounds on the reaction runaway hazard of H2O2 aqueous solution, and found that: iron, copper and chloride ions accelerated the thermal runaway of the decomposition of H2O2, however, nickel, potassium, sulfate and nitrate ions played an inhibitory role in the thermal runaway of the decomposition of H2O2. Wu et al.[6] mixed 20 wt.% H2O2 aqueous solution with different concentrations of sulfuric acid (H2SO4), utilized various thermal analyzers to study the thermal decomposition characteristics of H2O2 and obtained the thermodynamic parameters. It was found that the decomposition of H2O2 could be inhibited to a certain extent when the concentration of H2SO4 was 0.2, 1, and 2 mol/L. However, the addition of 96 wt.% H2SO4 to the solution led to a rapid increase in the rate of H2O2, and the whole exothermic reaction of the decomposition of H2O2 was violent. Adding stabilizers can improve the stability of H2O2. Domestic and foreign research on H2O2 stabilizer is mainly applied in the field of paper bleaching. The mechanism of H2O2 stabilizer includes complexation theory and colloidal adsorption theory. Complexation theory suggests that metal ions in solution and multivalent chelating agents can chelate to form stable water-soluble complexes, so that metal ions do not catalyze the decomposition of H2O2. And the theory of colloid adsorption suggests that heavy metal ions can be adsorbed through the polymer colloid nationality electrostatic or hydrogen bonding to achieve stability[7]. Zhang et al.[8] investigated the effect of adding different combinations of H2O2 stabilizers on the decomposition of H2O2 by free manganese ions and manganese ions in the presence of precipitated lignin, and the results of the study found that the addition of the H2O2 stabilizer, ethylenediaminetetraacetic acid (EDTA), alone significantly reduced H2O2 decomposition under both conditions. Jung et al.[9] evaluated the stability of H2O2 in the Fenton reaction using phthalic acid as a stabilizer. It was found that at low pH of the solution, iron ions (Fe3+) in the solution complexed with phthalic acid to form stable complexes, which inhibited the catalytic activity of Fe3+ and prolonged the service life of H2O2.

      In recent years, domestic and international scholars, based on a variety of technical means from the material characteristics, reaction process and mechanism analysis of the safety of the oxidation process related research, and according to the results of the study put forward corresponding improvement measures. Li[10] used a reaction calorimeter to investigate the effect of different charging methods on the thermal parameters during the reaction of H2O2 oxidation of acetic acid to prepare peroxyacetic acid. The 'one-pot' charging method was found to have a lower reaction risk than the 'drop-in' method. In addition, the secondary decomposition reactions of the products and feedstocks of the process were investigated using adiabatic accelerated calorimetry. The results of the study show that the possibility of secondary decomposition of peroxyacetic acid and H2O2 is higher, and the control should be focused on in the actual process. Liu[11] investigated the thermal hazards in the synthetic process for the production of N-methylmorpholine oxide by H2O2 oxidation using a simultaneous thermal analyzer (STA), an adiabatic accelerated calorimeter (ARC) and a fully automated reaction calorimeter (RC1e). It was found that the product N-methylmorpholine oxide has a certain instability, and the decomposition reaction has a large amount of heat release, which is potentially explosive. In the reaction process, compared with the intermittent charging, the semi-intermittent charging method can effectively avoid the accumulation of material in the system and reduce the probability of thermal runaway. Parker et al.[12] investigated the exothermic properties of the process of oxidation of 2-butanol by hydrogen peroxide for the preparation of 2-butanone using RC1e in terms of initial concentration, reaction temperature and catalyst dosage. The results showed that the initial concentration of 2-butanol had a more significant effect on the total exothermic amount of the reaction, on the basis of which it was concluded that this synthesis process should be carried out at a relatively low concentration of 2-butanol and moderate temperature and catalyst dosage.

      Microchannel reactors are structures in which chemical reactions are accomplished within a lateral dimension of less than 1 mm, which is different from the 'three passes and one reversal' of conventional reactors and offers significant advantages in chemical applications. Such as large specific surface area, good heat transfer efficiency, good mass transfer, precise control of reaction time and flexible production methods. In recent years, more and more domestic and international scholars have been studying it. Zhao et al.[13] used a microchannel reactor to develop a continuous flow process for the preparation of phenol from benzene oxidized by H2O2, considering the effects of material molar ratio, reaction temperature, catalyst concentration and residence time on phenol selectivity and yield, respectively. Under the conditions of H2O2 to benzene molar ratio of 1.7:1, reaction temperature of 55 °C, catalyst concentration of 0.01 mol∙L−1, and residence time of 60 min, 94% phenol selectivity was obtained, however, due to the low conversion of benzene the final phenol yield was only 15%. Constantinou et al.[14] designed and constructed a combined microchannel reactor combining a microchannel reactor and a membrane reactor for the synthesis of benzaldehyde by oxidizing benzyl alcohol with oxygen. In the experiments, by adjusting the residence time and reaction pressure, a benzyl alcohol conversion of 44.1% and a selectivity of 73.0% for benzaldehyde were obtained at the reaction temperature of 120 °C. At present, there are many studies on the thermal decomposition and oxidation reaction of hydrogen peroxide, but there is a lack of research on how to improve the stability of hydrogen peroxide in the reaction and how to improve the safety of the process of adipic acid preparation by oxidizing cyclohexene with hydrogen peroxide. For this reason, this study investigates the effect of stabilizer EDTA on the stability of hydrogen peroxide in the adipic acid preparation reaction and develops a microchannel continuous flow synthesis process for adipic acid to improve the intrinsic safety of the process.

    Experimentation and methods

      Experimental medications

    • This study involved 30% hydrogen peroxide solution, sodium tungstate, concentrated sulfuric acid, and ethylenediaminetetraacetic acid as shown in Table 1; all medications were used directly without further purification.

      Table 1.  Material information used in decomposition experiments.

      MaterialMolecular formulaFinenessCAS No.Manufacturer
      H2O2 aqueous solutionH2O230.0%7722-84-1Sinopharm Chemical Reagent Co., Ltd.
      Sodium tungstateNa2WO4≥ 99.5%10213-10-2Shanghai Macklin Biochemical Co., Ltd.
      Concentrated H2SO4H2SO4≥ 98%7664-93-9Sinopharm Chemical Reagent Co., Ltd.
      EDTAC10H16N2O8AR60-00-4TCI (Shanghai) Development Co., Ltd.

      Experimental instrumentation

      Reaction calorimeter

    • Relevant sensors can record the temperature, exothermic rate and other parameters during the whole reaction process in real time, and through the determination of the heat transfer coefficient and specific heat capacity, the thermodynamic parameters such as reaction enthalpy can be calculated by integration.

      Homemade high pressure microchannel reactor

    • The high-pressure microchannel reactor constructed in this study is shown in Fig. 2, which mainly consists of a premixing module, feeding module, reaction module, pressurization module and discharge module. The premixing module includes a digital thermostatic water bath and the three-necked flask to realize the premixing of the reaction materials. The feed module consists of a high-pressure advection pump ,which flow rate range is 0.1−80 ml/min and its pressure range is 0−7 MPa, a check valve (to prevent backflow of material during the reaction process) and a ball valve. The reaction module includes a 316 L stainless steel reaction coil (tube I.D. (1 mm) * O.D. (3 mm), the volume is 13.4 mL), and an oil bath (temperature range of −20 ~ 100 °C to provide reaction temperature for chemical reactions). The pressurization module consists of a nitrogen cylinder and a pressure-reducing valve (to provide the required working pressure for the reaction system). The discharge module consists of a filter (to adsorb fine particles in the reaction tube and prevent clogging of the HP backpressure valve), and a HP backpressure valve (to control the pressure in the tube).

      Figure 2. 

      High pressure microchannel reactor.

      Experimental methods

      Reaction calorimetric experimental method

    • In the experiment, the molar ratio of H2O2, cyclohexene and Na2WO4 was 4.4:1:0.06, the concentration of H2SO4 was 0.12 mol/L, the mass concentration of SDS was 9.15 g/ L and the mass concentration of EDTA was 1.83 g/L. 728.68 g of 30 wt.% H2O2, 28.908 g of Na2WO4, 7.96 g of concentrated H2SO4, 6 g of SDS, 12 g of EDTA and 120 g of cyclohexene were weighed. After all materials were mixed at room temperature, stirring was turned on, and the isothermal mode was selected to increase the temperature to 73 °C and reflux reaction for 2 h. After that, the reaction was continued to increase the temperature to 90 °C for 5 h, and then the experiment was finished.

      Experimental methodology for Na2WO4-catalyzed H2O2 decomposition reaction in microchannel reactor

    • The specific experimental setup is shown in Fig. 2, and the experimental steps are as follows:

      (1) Close the feed port ball valve 1, open the nitrogen port ball valve 2, nitrogen into the stainless steel coil then adjust the high temperature and high pressure back pressure valve, so that the pressure inside the pipe to maintain a stable 5MPa after closing the nitrogen port ball valve 2;

      (2) Prepare H2O2 decomposition reaction solutions for reaction conditions 2 and 3, respectively. Stirring was turned on in a water bath and an ice bath was used to minimize H2O2 decomposition;

      (3) Adjust the temperature of oil bath 90, 95, 100 °C and the flow rate of high-pressure advection pump 2.68, 1.34, 0.893, 0.67 ml/min, and the volume of microchannel reactor is 13.4 ml, corresponding to the theoretical retention time of 5, 10, 15, and 20 min, respectively. Open the inlet ball valve 1 to start the H2O2 decomposition experiments under 2 reaction conditions. opening the inlet ball valve 1, and start the H2O2 decomposition experiments under 2 reaction conditions, see Table 2.

      Table 2.  Information table of different reaction conditions.

      Reaction conditionn (H2O2) : n (Na2WO4)Molar concentration of H2SO4 (mol/L)Mass concentration of EDTA (g/L)
      14.4:0.060.22−
      24.4:0.060.2218.3
      Adipic acid continuous flow experiment

    • The experimental steps of the continuous flow synthesis process of adipic acid are as follows:

      (1) The molar ratio of H2O2 : C6H10 : Na2WO4 = 4.4:1:0.06, the concentration of H2SO4 is 0.12 mol/L, the mass concentration of surfactant SDS is 9.15 g/L, and the mass concentration of H2O2 stabilizer EDTA is 18.3 g/L were configured. The reaction is carried out at reflux for 120 min by turning on the stirrer in a water bath and raising the temperature to 73 °C.

      (2) Close the feed port ball valve 1, and open the nitrogen port ball valve 2, the nitrogen will be passed into the microchannel stainless steel reaction tube, adjust the high-temperature backpressure valve to control the pressure in the tube for 5, 6, 7 MPa waiting to stabilize after closing the nitrogen port ball valve.

      (3) Adjust the temperature of oil bath 100 °C and the flow rate of the high-pressure advection pump to ensure that the theoretical residence time is 5, 10, 15, 20, 25, and 30 min, respectively. Open the inlet ball valve 1 for the experiment.

    Reaction calorimetry experiment results and discussion

    • Table 3 shows the thermal parameters of the adipic acid synthesis reaction in the RC1e experiment, and Fig. 3 shows the graphs of temperature (Tr), jacket temperature (Tj) and heat flow rate (qr) in the reactor during the pilot scale-up of the adipic acid green synthesis reaction RC1e. Cyclohexene is gradually oxidized to epoxycyclohexane by H2O2 under reflux conditions starting at 73 °C, followed by rapid hydrolysis to 1,2-cyclohexanediol. As the cyclohexene is gradually consumed, the cooling effect due to reflux becomes less and less effective, so the temperature difference and heat flow rate of the reaction system gradually increase. The maximum temperature difference in stage A reach 14.69 °C. When the reflux reaction is finished and the temperature of the reaction system is increased to 90 °C, 1,2-cyclohexanediol is gradually oxidized by H2O2 and hydrolyzed to the final product adipic acid[15]. There is no reflux phenomenon in this process, the exothermic heat of the reaction is very large, and the maximum temperature difference in the B stage reaches 50.19 °C, therefore, in the process of industrialized application of green synthesis of adipic acid, the emergency cooling down of the B stage should be done.

      Table 3.  Thermal parameters of adipic acid synthesis reaction in RC1e experiment.

      Phase−ΔHr
      (kJ)      		Cp
      (J·K−1·g−1)      		ΔTad,r
      (°C)      		MTSRr (°C)
      (Adding EDTA)      		MTSRr (°C)
      (No EDTA)
      A103.213.837.23104.76209.77
      B504.774.1168.74268.81345.61

      Figure 3. 

      The Tr, Tj, and qr curves of adipic acid synthesis reaction in RC1e experiment.

      In this thesis, the thermal hazard analysis is performed for the two stages A and B of the green synthesis process of adipic acid, respectively. The severity of reaction loss of control is usually measured by the adiabatic temperature rise due to the exothermic nature of the target reaction, and the yield YAA of adipic acid was calculated to be 80.75% according to HPLC, which indicates that there are still unconverted intermediates, and when the system temperature is warmed up to the target temperature here, the reaction will still continue to exothermic. Therefore in order to assess the most hazardous case of the reaction process, the adiabatic temperature rise ΔTad needs to be corrected using the yield YAA[16]. ΔTad and ΔTad,r are calculated by Eqns (1) and (2). RC1e calorimetric test data and results are shown in Tables 4 and 5, where ΔHr is the heat of the oxidation reaction, kJ; mr is the total mass of the oxidation reaction system, g; and Cp is the specific heat capacity of the system after the oxidation reaction, J/(K·g).

      Table 4.  Average values of reaction rate constants at different temperatures in a stainless steel capillary microreactor calculated based on the first order reaction assumption under different reaction conditions.

      Reaction temperature (°C)Reaction condition 2 calculated rate constant average k1 (10−3 s−1)Reaction condition 2 calculated rate constant average k1 (10−3 s−1)
      900.25950.2082
      950.77720.5695
      1001.33040.9690

      Table 5.  Average values of reaction rate constants at different temperatures in a stainless steel capillary microreactor calculated based on the zero order reaction assumption under different reaction conditions.

      Reaction temperature (°C)Reaction condition 2 calculated rate constant average k0 (10−3 mol/(L·s)Reaction condition 3 calculated rate constant average k0 (10−3 mol/(L·s)
      902.33801.9032
      956.16594.7752
      1009.47327.38445
      ΔTad=ΔHrCpmr (1)
      ΔTad,r=ΔHrCpmrYAA (2)
      Tcf=Tp+ΔHr−∫t0qr(t)dtCpMt (3)
      Tcf,r=Tp+ΔHrYAA−∫t0qr(t)dtCpMt (4)

      In Eqns (3) and (4): Tp is the temperature in the kettle during the synthesis process, °C; Mt is the total mass of the oxidized reactive substance; Tcf is the temperature that can be reached at a certain moment after the cooling failure of the reaction system, and the maximum value of Tcf is the MTSR, and the Tcf is also corrected by the yield YAA; MTSR is the maximum temperature of the reaction and characterizes the severity of hardening out of control[17]. Figure 4 shows the variation of Tcf, Tcf,r in stages A and B, yielding a MTSRr of 104.76 °C in stage A and 263.81 °C in stage B. The MTSRr of stage A is 104.76 and 263.81 °C in stage B, respectively. However the process without the addition of the H2O2 stabilizer EDTA resulted in an MTSRr of 209.77 °C for stage A and 345.61 °C for stage B[15]. The values of MTSRr in this experiment are all reduced in comparison. Combined with the literature[18] to illustrate that the addition of H2O2 stabilizer EDTA can effectively reduce the severity of adipic acid green synthesis reaction out of control, with adipic acid green synthesis application prospects.

      Figure 4. 

      Tcf curves and MTSR of adipic acid synthesis reaction (a) Phase A; (b) Phase B.

    Experimental results and discussion of Na2WO4-catalyzed H2O2 decomposition reaction in the microchannel reactor

    • In the continuous flow process where H2O2 is involved, the higher reaction temperature causes the decomposition of H2O2 to produce oxygen. Gas expansion pushes the fluid through the microreactor quickly, which makes the gap between actual and theoretical residence times too large. The significant reduction in residence time leads to incomplete reaction, which can be addressed by increasing the reaction pressure to compress the gases produced by decomposition in order to minimize their effect on the residence time[19]. In this section, the effect of adding and not adding the H2O2 stabilizer EDTA on the decomposition reaction of H2O2 in a capillary microreactor is investigated. The H2O2 decomposition kinetics at high temperature and pressure in the microchannel reactor are studied to obtain the reaction rate constant k, the activation energy of the decomposition reaction Ea, and the finger-forward factor A.

      Figure 5 shows the relationship between the catalytic decomposition rate and the theoretical residence time of H2O2 for reaction conditions 1 and 2 with different reaction temperatures. The conversion of H2O2 under both conditions increase with increasing temperature and theoretical residence time. The conversion of H2O2 after 100 °C and theoretical residence time of 20 min is 60.72% and 51.07%, respectively. It can be seen from Fig. 5 that the addition of EDTA result in a decrease in H2O2 conversion at both the same temperature and flow rate. The higher the temperature, the more pronounced the inhibition of H2O2 decomposition by EDTA. Therefore, the addition of EDTA can also improve the utilization of H2O2 in the microchannel continuous flow process.

      Figure 5. 

      Effect of reaction temperature and residence time on hydrogen peroxide ecomposition in a stainless steel capillary microreactor. (a) Reaction condition 1; (b) Reaction condition 2.

      As the O2 is produced during the decomposition of H2O2, the actual residence time of the reaction cannot be determined in the absence of kinetic information. In this study, a plunger flow model is used to describe the extent of H2O2 decomposition in a stainless steel capillary microreactor, focusing on the effect of gas production on residence time (reaction time). For the continuous flow system, the average volumetric flow rate (v) of the fluid, in unit m3/s, can be expressed by Eqn (5):

      v=dVdτ (5)

      τ is the reaction time in unit s. Integration yields Eqn (6):

      ∫Vreactor0dV=∫τ0vdτ (6)

      For the decomposition of H2O2, v can be expressed by Eqn (7):

      v=vgas+vliquid (7)

      vgas is the volume flow rate of the gas phase (O2) produced; vliquid is the volume flow rate of the liquid. Equation (7) is brought into Eqn (6) to obtain Eqn (8):

      ∫Vreactor0dV=∫τ0(vgas+vliquid)dτ (8)

      For the resulting gas phase flow rate, the ideal gas-state equation is used to convert to obtain Eqn (9):

      vgas=nRTP=n0(H2O2)XRT2P (9)

      n0(H2O2) denotes the initial number of moles of H2O2 in the reaction mixture, mol, and the n0(H2O2) value in this experiment is 0.1312 mol; X is the conversion of H2O2 andRis the ideal gas constant of 8.314 J/(mol-K), PandTrepresent the reaction pressure (Pa) and the absolute temperature (K), respectively, and Eqn (9) is brought to Eqn (8) to obtain Eqns (10) and (11):

      ∫Vreactor0dV=∫τ0(vliquid+n0(H2O2)XRT2P)dτ (10)
      Vreactor=∫τ0(vliquid+n0(H2O2)XRT2P)dτ (11)

      Vreactor indicates that the volume of the microchannel reactor can be calculated from the inner diameter and length of the stainless steel tubing in mL. The volume of the microreactor in this decomposition experiment is 13.4 mL. The rate of reaction of H2O2 can be described by the stage rate equation, as shown in Eqn (12):

      −dcdτ=kcx (12)

      Where k is the reaction rate constant relative to H2O2, including the effect of the catalyst on it, and is the concentration of H2O2 in the liquid phase in mol/L. In the literature, H2O2 decomposition is mostly considered as a primary reaction equivalent to H2O2[20,21]. If the first order rate law is followed (Case 1) Eqn (12) can be transformed into Eqn (13):

      dcc=−k1dτ (13)

      Integrating Eqn (13) yields,

      ln(c0c)=k1τ (14)

      The concentration of H2O2 can be calculated based on the initial concentration and the conversion rate which is shown in Eqn (15), The H2O2 in this decomposition experiment is 30 wt.% converted to give a c0of 9.79 mol/L. Defaulting to the fact that there is no change in the volume of the aqueous phase during H2O2 consumption, Eqn (15) is brought into Eqn (14) to obtain Eqns (16) and (17). Taking the derivative of Eqn (17) with respect to X yields Eqn (18):

      c=c0(1−X) (15)
      ln(c0c0(1−X))=k1τ (16)
      1k1ln(11−X)=τ (17)
      dτ=1k1(1−X)dX (18)

      Eqn (19) is obtained by combining Eqns (18) and (11). With known temperature T, H2O2 conversion X and inlet liquid volume flow rate, the integral term is calculated from experimental data and Eqn (20).

      Vreactor=∫X0(vliquid+n0(H2O2)XRT2P)[1k1(1−X)]dX (19)
      k1Vreactor=∫X0(vliquid+n0(H2O2)XRT2P)[1(1−X)]dX (20)

      The function of k1*Vreactor and H2O2 conversion rate X is calculated according to Eqn (20), and the results are shown in Fig. 6. From Fig. 6 it can be seen that k1*Vreactor depends strongly on the H2O2 conversion rate. Thus for the catalytic decomposition of H2O2 under reaction conditions 1 and 2, the first-order rate fixing rate in H2O2 does not seem to be observed. The average values of the reaction rate constants (k1) at the three reaction temperatures are shown in Table 4.

      Figure 6. 

      Relationship between k1*Vreactor and H2O2 conversion rate. (a) Reaction condition 1; (b) Reaction condition 2.

      In the literature, the decomposition of H2O2 is zero order compared to H2O2 when the initial concentration of H2O2 is 35.5 wt.%[22]. The initial concentration used for the experiment is 30 wt.%, so assuming that the decomposition reaction in H2O2 is zero-order, following the zero-order rate law (Case 2), Eqn (12) can be transformed into Eqn (21), and Eqn (15) is brought into Eqn (21) to give Eqn (22). Equation (22) and (11) are combined to obtain Eqn (23).

      −dcdτ=k0 (21)
      dτ=c0dXk0 (22)
      k0Vreactor=c0∫X0(vliquid+n0(H2O2)XRT2P)dX (23)

      According to Eqn (23), the function of k0*Vreactor and H2O2 conversion rate X is calculated, and the results are shown in Fig. 7. The value of k0*Vreactor is not related to the H2O2 conversion within the experimental error. The results show that H2O2 decomposition under reaction conditions 1 and 2 is a zero-order reaction. The average values of the reaction rate constants (k0) at the three reaction temperatures are shown in Table 5.

      Figure 7. 

      Relationship between k0*Vreactor and hydrogen peroxide conversion rate. (a) Reaction condition 1; (b) Reaction condition 2.

      To further verify that Case 2 (zero-stage reaction of H2O2) is more reasonable than Case 1 (one-stage reaction of H2O2) in the catalytic decomposition of H2O2 under reaction conditions 1, 2, the standard deviation (SD) between k for each residence time and its corresponding mean value at the same reaction temperature is calculated using Eqn (24):

      SD=√1n−1n∑i=1(ki−km)2=√134∑i=1(ki−km)2 (24)

      From Tables 6 & 7 it can be seen that for different temperatures based on zero level reaction the SD values are below 5%. These SD values are significantly lower than the SD assumed for the first-order reactions. The values obtained for the 95% confidence intervals are also closer to the reaction rate constant k values. Thus it can be concluded that in a stainless steel capillary microreactor the initial H2O2 concentration of 30 wt.%. The H2O2 decomposition reaction with Na2WO4 as catalyst under reaction conditions 1 and 2 is a zero-order reaction compared to H2O2.

      Table 6.  Comparison of the calculated reaction rate constant k, standard deviation (SD), and confidence interval of the calculated reaction rate constant k based on the assumption of first and zero order reactions in a stainless steel capillary microreactor under reaction condition 2.

      Reaction temperature (°C)k1(10−3 s−1)SD95% confidence interval of k1k0 (10−3 mol/L·s)SD95% confidence interval for k0
      900.25950.10450.2471−0.27192.33800.04932.2806−2.3953
      950.77721.07540.6495−0.90496.16590.02256.1397−6.1920
      1001.33042.90350.9856−1.67529.47320.02719.4417−9.5048

      Table 7.  Comparison of the calculated reaction rate constant k, standard deviation (SD), and confidence interval of the calculated reaction rate constant k based on the assumption of first and zero order reactions in a stainless steel capillary microreactor under reaction condition 3.

      Reaction temperature (°C)k1 (10−3 s−1)SD95% confidence interval of k1k0 (10−3 mol/L·s)SD95% confidence interval for k0
      900.20820.07190.1967−0.21961.90320.00981.8875−1.9188
      950.56950.48710.4920−0.64701.90320.00684.7644−4.7861
      1000.96900.12010.7778−1.16014.77520.02057.3518−7.4171

      According to the Arrhenius equation, the relationship between the reaction rate constant and temperature can be expressed as Eqn (25).

      lnk=lnA−EaRT (25)

      A is the prefactor, Ea is the activation energy of the reaction, and the fit of lnk to 1/T is shown in Fig. 8.

      Figure 8. 

      Experimental value and fitting curve of lnk to 1/T in hydrogen peroxide decomposition. (a) Reaction condition 1; (b) Reaction condition 2.

      The activation energy Ea of the H2O2 decomposition reaction under reaction condition 1 is calculated to be 157.665 kJ/mol from the slope of the straight line in Fig. 8, and the intercept is calculated to give a pre-finger factor A of 1.128 × 1023 mol· L−1·s−1. The value of this activation energy is significantly higher than that obtained in Gaussian calculations and reaction calorimetry of H2O2 decomposition in the reaction system, because the addition of H2SO4 inhibits the H2O2 decomposition cycle, so that the value of the activation energy of its reaction is high. The activation energy Ea of the H2O2 decomposition reaction under reaction condition 2 is 164.428 kJ/mol, which is higher than that derived from reaction condition 1, indicating that the addition of the H2O2 stabilizer EDTA could be able to reduce the decomposition of H2O2 to a certain extent, and the intercept calculation yielded the finger forward factor A is 0.848 × 1023 mol· L−1·s−1.

    Experimental results and discussion on adipic acid microchannel continuous flow synthesis

    • From the results in the previous section, it can be seen that although the addition of H2O2 stabilizer EDTA can improve the stability of hydrogen peroxide and reduce the temperature of MTSR, the reaction is still prone to thermal runaway, resulting in the violent decomposition of H2O2 and thus triggering the phenomenon of flushing material. Therefore, in order to further enhance the intrinsic safety of green synthesis process of adipic acid, the batch process of green synthesis of adipic acid is improved into a continuous flow synthesis process by utilizing the characteristics of microchannel reactor with small liquid holding capacity and excellent mass and heat transfer, with the aim of realizing intrinsic safety of the process.

      In the process of the preliminary experimental exploration, it was found that under the conditions of atmospheric pressure and reaction temperature, the reaction material in the reaction tube is violently decomposed and gasification of the material, resulting in the residence time of the reaction is far less than the theoretical residence time. Increasing the residence time by increasing the reaction pressure and thus compressing the volume of gas produced by decomposition is a common solution. The ratio of bubbles to liquid in the tubes for different reaction pressures is shown in Fig. 9. Therefore, in this study, the reaction pressure of 5−7 MPa and the reaction temperature of 100 °C were chosen for the microchannel continuous flow synthesis experiments with different residence times.

      Figure 9. 

      Schematic diagram of gas to liquid ratio in pipelines.

      It can be seen from Fig. 10 that with the reaction pressure increasing so that the gas produced in the pipeline is gradually compressed, the actual residence time of the reaction increases, and the yield of adipic acid also gradually increases. In addition, the actual residence time of the reaction is also influenced by the theoretical residence time, which is constantly increased to make the yield of adipic acid also increase gradually. Finally, the yield of adipic acid is 63.25% at 7 MPa reaction pressure and 30 min theoretical residence time, and the reaction time of the microchannel continuous flow experiment was greatly shortened compared with the reaction calorimetry time.

      Figure 10. 

      Green synthesis yield of adipic acid under different pressures and residence times.

    Conclusions

    • (1) Through the calorimetric experiments of adipic acid synthesis, the MTSR of phase A (isothermal stage at 73 °C) is 104.76 °C, and that of phase B (concentrated exothermic stage in the first 3 h at 90 °C) is 263.81 °C, which is reduced compared with that of the process without adding EDTA.

      (2) Through the study of EDTA in continuous flow Na2WO4-catalyzed H2O2 decomposition, it is found that the decomposition of H2O2 increases with the increase of residence time and temperature, and the addition of EDTA can effectively reduce the decomposition. The activation energies for the decomposition of 30 wt.% H2O2 with and without EDTA in stainless steel capillary microchannels are 164.428 and 157.665 kJ/mol, respectively.

      (3) Continuous flow synthesis process of adipic acid in a homemade high-pressure microchannel reactor, obtaining the highest yield of adipic acid up to 63.25%. The microchannel continuous flow process shortened the reaction time of the intermittent experiment in terms of reaction time, but the yield is significantly lower than that of adipic acid in the intermittent experiment of 82.79%, mainly because of the waste of reaction materials and shortening of reaction time due to the violent decomposition of H2O2. It is recommended to pay attention to H2O2 decomposition when subsequently optimizing the continuous flow process.

    Author contributions

    • The authors confirm contributions to the paper as follows: study conception and design: He W, Ni L; data collection: Li Y; analysis and interpretation of results: He W, Ni L; draft manuscript preparation: Li Y, Zhu W, Ni L. All authors reviewed the results and approved the final version of the manuscript.

    Data availability

    • All data generated or analyzed during this study are included in this published article.

      Acknowledgments

      • This work is supported by the National Natural Science Foundation of China (52274209).

      Conflict of interest

      • The authors declare that they have no conflict of interest.

      Rights and permissions
      • Copyright: © 2023 by the author(s). Published by Maximum Academic Press on behalf of Nanjing Tech University. This article is an open access article distributed under Creative Commons Attribution License (CC BY 4.0), visit https://creativecommons.org/licenses/by/4.0/.
    Figure (10)  Table (7) References (22)
  • About this article
    Cite this article
    He W, Li Y, Ni L, Zhu W. 2023. Effect of stabilizer EDTA on the thermal hazard of green synthesis process of adipic acid and development of microchannel continuous flow process. Emergency Management Science and Technology 3:22 doi: 10.48130/EMST-2023-0022
    He W, Li Y, Ni L, Zhu W. 2023. Effect of stabilizer EDTA on the thermal hazard of green synthesis process of adipic acid and development of microchannel continuous flow process. Emergency Management Science and Technology 3:22 doi: 10.48130/EMST-2023-0022
    shu

Catalog

  • HTML

  • Introduction
  • Experimentation and methods

  • Experimental medications

  • Experimental instrumentation

  • Reaction calorimeter
  • Homemade high pressure microchannel reactor

  • Experimental methods

  • Reaction calorimetric experimental method
  • Experimental methodology for Na2WO4-catalyzed H2O2 decomposition reaction in microchannel reactor
  • Adipic acid continuous flow experiment
  • Reaction calorimetry experiment results and discussion
  • Experimental results and discussion of Na2WO4-catalyzed H2O2 decomposition reaction in the microchannel reactor
  • Experimental results and discussion on adipic acid microchannel continuous flow synthesis
  • Conclusions
  • Author contributions
  • Data availability
  • Acknowledgments
  • Conflict of interest
  • Rights and permissions
  • About this article
logo

Contents

Policies

Services

Partnerships

© 2025 Maximum Academic Press, except Open Access articles. MAP is a member of ALPSP, COPE, CrossRef, cOAlition S, OASPA and STM.

Privacy Policy  |  Site Map  |  Terms and Conditions

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return