Comparative anatomical and transcriptomic analyses of the color variation of leaves in Aquilaria sinensis

Plant materials

The golden-vein-leaf and normal-vein-leaf sample (GS and NS) A. sinensis plants were obtained from Dalingshan town, Guangdong province, China (latitude 22°45′43″N, longitude 113°48′45″E). Those plants grow in close proximity with similar environment and at a similar altitude of 70–80 m. The samples were all collected from six A. sinensis trees in which three trees were GS and the others were NS. Then, three leaves and three stems of GS as well as those of NS were collected from each three branches of trees and mixed respectively, with three biological replicates used for each group in DNA barcode sequencing and RNA sequencing. All of the samples were frozen in liquid nitrogen and stored at −80 °C for the next studies.

DNA extraction and determination of chlorophyll content

The leaves were powdered in a mortar with a pestle in liquid nitrogen. The stems were sliced first, and then powered in liquid nitrogen according to the published methods (Jiao et al., 2014). Additionally, those samples were inserted into a 2 mL tube containing 1,000 µL buffer AP1, 8 µL RNase A, and 1% polyvinylpyrrolidone with incubation at 65 °C for 6 h respectively. Cooling to room temperature, the 280 µL buffer P3 was added into the tubes and those samples were incubated for 2 h at −20 °C. Subsequent steps were conducted as the direction supplied by the manufacturer of EASYspin Plus Plant DNA Kit (Aidlab Biotechnologies Co., Ltd, Beijing, China). The determination of chlorophyll content was performed based on the method published before (Sartory & Grobbelaar, 1984).

DNA barcode sequence amplification and phylogenetic tree analysis

ITS2 and trnL-trnF were used as barcode sequences to identify the species and the primers were adopted from a published paper (Lee et al., 2016). The PCR reaction system was also conducted as the previous report (Lee et al., 2016).

The ITS2 and trnL-trnF barcode sequences of GS and NS were aligned with A. crassna, A. hirta, A. malaccensis, A. microcarpa, A. sinensis, A. yunnanensis, Gyrinops versteegii and Gonystylus bancanus (Table S1) using Clustal W algorithm (Larkin et al., 2007). Subsequently, the phylogenetic tree was constructed based on neighbor-joining (NJ) method with 1,000 bootstrap replications. All the above analyses were operated with MEGA-X (10.1.8) (Kumar et al., 2018).

Microscopic observation and measurement

The paraffin microsections were sliced by Leica CM1860 freezing microtome (Leica Microsystems Inc., Wetzlar, Germany) with Safranin O/Fast green staining subsequently. The microscopic observations of those microsections were conducted by ZEISS AX10 microscope (ZEISS corporation, Jena, Germany). ZEISS ZEN 2 lite software was applied to measure the parameters of the photomicrographs. The statistical differences were calculated by T-test using IBM SPSS Statistics 22. The ultrastructure was observed through transmission electron microscope H-7500 (Hitachi, Ltd., Japan).

RNA extraction and transcriptome sequencing

The RNA was extracted from the A. sinensis samples using the EASYspin Plus Plant RNA Kit (Aidlab Biotechnologies Co., Ltd, Beijing, China) according to the manufacturer’s instructions. The quality of extracted RNA was examined by NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA). Then, 0.5 µg RNA per sample was reverse-transcribed into single-stranded cDNA using cDNA Synthesis Kit (TaKaRa, Japan) as what the manufacturer’s instructions suggested. Subsequently, the RNA template was eliminated and the second-strand cDNA synthesis was conducted with DNA polymerase, dNTPs, and RNase. After end repair, A-tailing and indexing ligation, the products were purified with PCR extraction kit, and amplified to be cDNA libraries. Finally, the sequencing was conducted on Illumina NovaSeq 6000. The data have been deposited in NGDC’s Genome Sequence Archive (Wang et al., 2017), under accession number CRA002670 that is publicly accessible at https://bigd.big.ac.cn/gsa/browse/CRA002670. The de novo assembly was conducted by the Trinity method (Grabherr et al., 2011).

Annotation and differentially expressed gene analysis

The annotations of transcriptome data were conducted by Diamond (Buchfink, Xie & Huson, 2015), KAAS (Moriya et al., 2007) and Blast2GO (Gotz et al., 2008). The gene expression level was estimated with FPKM (Trapnell et al., 2010). Using DESeq2 (Love, Huber & Anders, 2014), when the parameters of genes expressed in different groups at log₂(Fold Change) larger than 1 or less than −1, and false discovery rate (FDR) less than 0.05, they could be accepted as differentially expressed genes (DEGs). The GO and KEGG enrichments were conducted via R package clusterProfiler (Yu et al., 2012) with the custom data sets from the annotation results. For GO enrichment, 15,811 annotated genes were used as the background gene set. For KEGG, 4,865 annotated genes were used as the background gene set. In GO enrichment, the FDR screening threshold was set to be less than 0.05, while the value was set to be less than 0.1 in KEGG enrichment. The cluster analysis was conducted via the STEM clustering algorithm (Ernst, Nau & Bar-Joseph, 2005) and the cluster whose P value was less than 0.05 was considered as significant difference. All the plots were generated using one public on-line tool BIC (http://www.ehbio.com/Cloud_Platform/front/) or R scripts.

Quantitative real-time PCR

The qRT-PCR was carried out with the TB Green^@ Fast qPCR Mix (TaKaRa, Japan) and LightCycler^@ 480 Real-Time PCR System (Roche, Switzerland). The reaction mixture contained 3 µL of nuclease-free water, 5 µL of TB Green Premix Ex Taq II (TaKaRa, Japan), 1 µL of cDNA, 0.5 µL of forward primer (10 µM), and 0.5 µL of reverse primer (10 µM). The program was set to run for 1 min at 94 °C, 40 cycles of 10 s at 94 °C, and 34 s at 60 °C. Melt curves were obtained by heating the samples from 60 °C to 95 °C at a rate of 1.0 °C/s (Table S9). The primers of genes were listed in Table S9 with a reference gene GAPDH (Xu et al., 2016). Each RNA sample was conducted with three biological repeats and three technical replicates. The relative gene expression levels were calculated through 2^−ΔΔCT method (Livak & Schmittgen, 2001). The correlation analysis was conducted by Pearson correlation coefficient.

DNA barcode sequencing and phylogenetic tree analysis

The samples used in this study were all collected from Guanxiang Intangible Cultural Heritage Protection Park, while their phenotypes differed from each other. The golden-vein-leaf sample (GS) showed a color variation of leaf vein compared to the normal-vein-leaf sample (NS) and had deeper green and harder texture in leaves (Figs. 1A and 1B), and the contents of chlorophyll including chlorophyll a, chlorophyll b and carotenoid in the leaves of GS (LGS) were all higher than those in the leaves of NS (LNS) (Fig. 1C). Though the selected samples were cultured as the same species in the cultivation place, the phenotypes differences were so remarkable that the DNA molecular identification was employed to confirm the phylogeny relationships between GS and NS (Fig. S1). According to the previous report, the combination of ITS2+trnL-trnF was used as DNA barcode for Aquilaria species identification (Lee et al., 2016). The phylogenetic tree was constructed based on ITS2+trn L-trnF DNA barcode sequences in 8 closely related species, and one distantly related species Gonystylus bancanus as outgroup according to the previous study (Lee et al., 2016) (Fig. 1D). Most of the selected reference species had their separate clades respectively indicating success of the phylogenetic tree construction. The high-level homology of GS and NS sequences with A. sinensis verified that they were all A. sinensis, but GS samples could not be discriminated clearly with NS samples, which suggested that the variety was not the cause of the GS formation.

Anatomical comparison of leaves and stems

The photomicrographs of the LGS and LNS were acquired after Safranin O/Fast green staining, of which the anatomical parameters were measured (Fig. 2; Table S2). In the transverse section of the midrib, the bicollateral bundles were apparent (Figs. 2A and 2C). It was visible that the total area of included phloem in LGS midrib surpassed which in LNS midrib for more than 5 times (Fig. 2E). In addition, the palisade tissue was thicker in LNS whereas the spongy tissue was thicker in LGS (Figs. 2B, 2D, 2F and 2G). Furthermore, the upper epidermis cell thickness in LGS is significantly smaller than LNS, while the difference of lower epidermis cell between LGS and LNS was not significant (Figs. 2H and 2I; Table S2).

Using transmission electron microscopy, the ultrastructure of chloroplast was observed (Fig. 3). It was obvious that the structure of chloroplast in the midrib of LNS was more complete than that in LGS in which the thylakoids were seriously broken (Fig. 3B). Apart from the thylakoid, the starch could be observed in the midrib of LNS (Fig. 3A), however, it was seldom in the midrib of LGS. In addition, we found some particles which were suspected as virus (Figs. 3C and 3D), and some vesicles (Fig. 3C).

With respect to the stems of GS and NS (SGS and SNS), the differences were also significant and the most obvious characteristics observed in transverse surface had been measured and digitized (Fig. 4; Table S3). One of the most obvious features in SGS was that the included phloem tended to connect with each other and their shapes were more irregular than that of SNS (Figs. 4A and 4E) which was typical island-shaped (Liu et al., 2018). In addition, the area, length and width of included phloem of SNS were all greater than that of SGS (Fig. 4I; Table S3). Because most of the included phloem distributed in xylem, we calculated the percentage of included phloem in xylem of SGS and SNS and found that the number of the former was a little bigger though it was less than 2% (Fig. 4J), indicating that the difference of distribution level of included phloem in SGS and SNS was not as significant as the area. In addition, many crystals could be observed clearly in the pith of SNS (Fig. 4G) while much less in the pith of SGS in which there were more colored grease-like granular substances (Fig. 4C).

The parameters about vessel of SGS were all smaller than SNS (Figs. 4K and 4L; Table S3). In consideration of the features that single vessel area of SNS was larger than that of SGS and the difference between the percentage of vessel to xylem of SGS and SNS was less than 2%, the density of SGS vessel could be higher (Figs. 4B and 4F). Furthermore, the single fiber cell of SGS was significantly smaller than that of SNS (Fig. 4M) and it could be observed intuitively (Figs. 4D and 4H).

Transcriptome sequencing and differential gene expression profile

Totally 437 million high quality reads (65 G bases) were generated for four samples, LGS, LNS, SGS, and SNS (each with three biological duplicates) (Table S4). Due to the lack of gene annotation information, de novo assemble were performed to get 26,381 unigenes. The contig N50 of all transcripts was 2,168 nt, the average transcript length was 1418 nt, and the complete BUSCOs (Benchmarking Universal Single-Copy Orthologs) was more than half of the total (Table S5; Fig. S2). Functional annotation showed that 18,786 genes could be annotated to at least one of Pfam, SwissProt, TrEMBL, KEGG and GO databases (Fig. S3).

Differential gene expression analyses were performed between LGS and LNS, SGS and SNS using negative binomial models. These genes with multiple test corrected P-value less than 0.05 and absolute fold change no less than two were defined as DEGs. In summary, 1525 genes were upregulated while 752 genes were downregulated in LGS compared to LNS (Fig. 5A), and 808 genes were upregulated while 245 genes were downregulated in SGS compared to SNS (Fig. 5B). In addition, the samples from same superset could be clustered together using principal component analysis (PCA) of DEGs which were more varianced compared to the no-differential expressed genes (Fig. S4).

Photosynthesis related genes were enriched in color variation of A. sinensis by GO term enrichment analysis of DEGs

Genes upregulated in LGS compared to LNS were significantly enriched in photosynthesis related biological processes including photosynthesis, light harvesting in photosystem I, response to light stimulus (Fig. 5C). Besides, defense response and carbohydrate metabolic process were two other enriched items occupying the highest proportion of DEGs. Those upregulated photosynthesis related terms may explain the deeper green in LGS and carbohydrate metabolic process was another important term participating in photosynthesis (Atkin et al., 2000). In addition, some terms related to cell wall were also enriched, such as pectin catabolic process and plant-type cell wall organization. It was worth considering that the phenotype of midrib may in connection with pectin catabolism. As for cellular component, genes upregulated in LGS enriched in components referred to photosynthesis. Except for “extracellular region”, the GO term “photosystem I” and “photosystem II” showed the most significance among these items. Meanwhile, “chloroplast thylakoid membrane” and “chloroplast envelope” covered a large amount of genes. GO terms “cell wall” and “plant-type cell wall” were also enriched which coincides with the terms about cell wall in biological process. For molecular function, the most significantly enriched terms about upregulated genes in LGS were related to chlorophyll such as “chlorophyll binding” and “pigment binding” as well as many genes were enriched into glycosyl hydrolysis relevant terms like “hydrolase activity, hydrolyzing O-glycosyl compounds”.

With regard to the downregulated genes (Table S6), the biological process GO terms “response to abscisic acid” and “ethylene biosynthetic process” coincided with phytohormone playing important roles in plant development (Emenecker & Strader, 2020; Iqbal et al., 2017). Some terms in molecular function involved transporting, such as “xenobiotic transmembrane transporting ATPase activity” and “ATPase activity, coupled to transmembrane movement of substances”. In summary, the enrichment terms of downregulated genes were much less than upregulated genes so that there were only three molecular function terms, and two cellular component terms enriched into “plasmodesma” and “plant-type vacuole”.

Comparing the GO enrichment analysis for DEGs between SGS and SNS, the significantly enriched GO terms of upregulated genes in SGS were shown in Fig. 5D. For biological process, a lot of upregulated genes in SGS were categorized into the terms which were bound up with translation, energy transporting and metabolism such as “translation”, “cytoplasmic translation”, “ATP synthesis coupled proton transport” as well as some terms related to ribosome. Focusing on cellular component, many upregulated genes were about ribosome among which the term “ribosome” was the most significant whose gene number was also the most. Moreover, some terms about mitochondrion were enriched too. For molecular function, the most significantly enriched term “structural constituent of ribosome” occupied a dominant quantitative advantage over others echoing with the situation in biological process and cellular component.

The number of downregulated genes which could be enriched was small (Table S7). The terms in cellular component category were not enriched because of the lack of significance and only one term “nitrate transport” was enriched in biological process. The three most significant terms in molecular function were “ADP binding”, “sucrose alpha-glucosidase activity” and “double-stranded RNA binding”.

Photosynthesis related pathways were characterized in color variation of A. sinensis by KEGG pathway enrichment analysis of DEGs

To analyze the gene functional pathways of DEGs in leaves and stems, the KEGG pathway enrichment analyses were conducted in leaves and stems samples respectively. It was found that only upregulated genes of GS compared to NS were enriched successfully together with the small amount (Fig. 5E). Although the number of genes enriched was much smaller than that of GO enrichment, it also attached importance to analyze the meaningful pathways such as the “photosynthesis - antenna proteins”, “photosynthesis”, “porphyrin and chlorophyll metabolism” pathways which were related to photosystem obviously. Furthermore, the other pathways associated with saccharide and steroid were worth observing too.

Compared with SNS, 74 genes were enriched in “ribosome” pathway which occupied the largest number in upregulated genes of SGS. The remaining three pathways enriched were “oxidative phosphorylation”, “cutin, suberin and wax biosynthesis” and “glyoxylate and dicarboxylate metabolism” respectively.

Analysis of transcription factors in DEGs

Transcription factors play important roles in plant development in response to biotic or abiotic stimulus through regulating gene expression of abundant genes (Singh, Foley & Onate-Sanchez, 2002). To explore the expression of transcription factors, the transcription factors expressed differentially in leaves or stems of samples were screened as shown in Fig. 5F. There were 46 transcription factors expressed differentially in leaf samples and the number in stem samples was 18. In total transcription factors, bHLH took up the largest number, and they are important in phytochrome signaling according to the previous studies (Duek & Fankhauser, 2005). With respect to other transcription factor family genes, they also played roles in the formation of the phenotype of our samples. For example, the HB family gene HAT5 (ATHB1) could emit its effect in the conversion of palisade to spongy tissue cells in the leaves of Arabidopsis thaliana (Aoyama et al., 1995). In addition, the G2-like gene GLK1 was another important transcription factor which could regulate the chloroplast development (Gang et al., 2019a) involved in the leaf etiolation (Xie et al., 2018). Additionally, since suberin biosynthetic process in SGS could be in connection with the formation of the phenotypes in Fig. 4 because wood fiber tissue accounted for the majority of stem and cutin and suberin play important roles in the formation of wood (Kolattukudy, 2001), some NAC and MYB family genes may be active in the regulation of those process (Zhong & Ye, 2007).

Pathogenesis-related (PR) genes significantly upregulated in golden-vein leaves and stems

The yellow vein symptom of leaves is widely distributed in the plants infected with pathogens. For instance, the genus Begomovirus is the largest geminivirus genus containing more than 200 species or members (Fauquet et al., 2003; Fauquet et al., 2008) of which some viruses have infected a lot of plants with yellow vein symptoms in south China (Jiao et al., 2013) which is the place we got our samples from. Dissimilar to the common phenotypes of the infected showing only yellow vein symptom or accompanied with fading, LGS showed a deeper green beside the yellow vein symptom compared to LNS. In consideration of the enrichment of “defense response” and “response to biotic stimulus” especially the latter in LGS, most of the genes in “response to biotic stimulus” belonged to pathogenesis-related (PR) gene family (Table 1) which was associated to infection and defense (Van Loon, Rep & Pieterse, 2006). Using the annotation information, 25 PR family genes expressed differentially in leaves or stems were screened. In those genes, 19 genes were upregulated while one gene was downregulated in LGS compared to LNS as well as five upregulated and two downregulated genes in SGS compared to SNS (Fig. 6A; Table S8). Meanwhile, as shown in Fig. 6B, most of the genes distributed in the third quadrant so that it was clear that most of the PR genes were upregulated in LGS and SGS. Using cluster analysis, 7 genes (PRB1, PR-1, E137, PER3, PER47, STH-2, NLTP) were grouped under one cluster in which the gene expressions were both higher in LGS and SGS than their corresponding groups (Fig. 6C). The upregulation of those genes in same expression pattern may be in connection with the infection and the formation of the phenotypes in our samples.

Quantitative real-time PCR verification of transcriptome data

The expression of those genes in different tissues were identified via qRT-PCR (Figs. 6E and 6F). Based on the enrichment data, those 23 selected genes used in qRT-PCR referred to photosynthesis, cell wall constituent, and some about secondary metabolism (Fig. 6D). Our data revealed that most of the gene expression data in qRT-PCR matched well with transcriptome results (Fig. S5), which verified that photosynthesis related pathways were involved in the color variation phenotypes of A. sinensis.