Healthy “Reyan No. 1” were collected from the areca germplasm resources nursery of the Coconut Research Institute, Chinese Academy of Tropical Agricultural Sciences (Wenchang, Hainan Province; longitude 110gitude,a E, latitude 19titudee, N) as experimental material. Three plants with relatively uniform growth were selected from different directions of the site. For each selected plant, 20-30 fruits on the bunches were bagged, while the untreated portion was used as CK. Harvesting was conducted 120 days after pollination for subsequent experiments.

The double-layer inner black and outer yellow paper bags selected for the test measured 31 × 23cm and were purchased from Green Farmers Fruit Bags.

About 45 days after pollination of areca nut female flowers, spray insecticide to kill insects, the insecticide used is thiamethoxam and chlorpyrifos mixed application, thiamethoxam concentration of 0.1%, chlorpyrifos concentration of 0.15%. Then hang the pest avoidance agent bag, the pest avoidance agent is p-dichlorobenzene (content 99%), hanging the pest avoidance agent into the bag tied to the near set of fruit stalks, bagging, the pest avoidance agent bag together in the bag. Cut off the male flowering branches at the front end of the spike before bagging, and the fruits need to be harvested about 120 days after pollination when the bags can be removed.

The longitudinal diameter, transverse diameter, and weight of the freshly harvested areca nut, as well a color difference of exocarp, were measured and presented as average values.

The chromaticity color parameters of the exocarp include L^*, a^* and b^*. The fruit surface color composite index (CCI) was utilized to indicate color saturation, and was calculated by the formula CCI =(1000 × a*)/(L* × b*). L* denotes brightness, a^* value represents the red-green color difference of the exocarp, and b^* value represents the yellow-blue color difference of the exocarp (Wang et al., 2023).

The extraction method of chlorophyll and carotenoid content of areca nut exocarp was done according to the method described by previous research (Li S. et al.,2017). Three portions of fresh areca nut exocarp were weighed from both the YP and CK (mixed sample of fruit exocarp, 0.5 g each). Initially, calcium carbonate and quartz sand were used for grinding. Then the powder was further ground with liquid nitrogen and extracted using 80% acetone.

Following the freeze-drying of areca nut fruit exocarp samples in a vacuum freeze dryer (Scientz-100F), the lyophilized exocarp samples were ground (30 Hz,1 min) to powder using a ball mill (MM 400, Retsch). Subsequently, 50 mg of each sample was weighed and extracted with 0.5 mL of a hexane/acetone/ethanol mixture (1:1:1, v/v/v) containing 0.01% BHT (g/mL). After vortexing (room temperature) for 20 min, the mixture was centrifuged (4 °C, 12000 rpm) for 5 min. The supernatant was collected and the extraction was repeated once, mixing the resulting supernatant. The collected supernatants were concentrated, re-dissolved with methanol/methyl tert-butyl ether mixture (1:1, v/v, 100 μL), passed through filter membrane (0.22 μm), and consequently stored in a brown sample bottle for LC-MS/MS analysis. The data acquisition instrumentation system primarily consisted of Ultra Performance Liquid Chromatography (UPLC) (ExionLC™AD, https://sciex.com.cn) and Tandem Mass Spectrometry (MS/MS) (QTRAP^®6500+, https://sciex.com.cn/) (Inbaraj et al., 2008; Zhao et al., 2022).

LC was performed according to the following conditions: column: YMC C30 (3 μm, 100 mm × 2.0 mm i.d.); mobile phase: A phase = methanol/acetonitrile (1:3, v/v) with 0.01% BHT and 0.1% formic acid; B phase = methyl tert-butyl ether with 0.01% BHT; gradient elution program: 100:0 (V/V) for 0 min A/B, 100:0 (V/V) for 3 min, 30:70 (V/V) for 5 min, 5:95 (V/V) for 9 min, 100:0 (V/V) for 10 min, and 100:0 (V/V) for 11 min. The flow rate, column temperature, and injection volume were 0.8 mL/min, 28 °C, and 2 μL, respectively.

The mass spectrometry conditions consisted primarily of an atmospheric pressure chemical ionization source (APCI, 350 °C) and curtain gas (CUR, 25 psi). In the Q-Trap 6500+, per ion pair was scanned and detected based on optimized de-clustering potential (DP) and collision energy (CE) (Mase et al., 2022).

Carotenoid standards at various concentrations were prepared, and the mass spectral peak intensities of the quantitative signals of each concentration were obtained for each concentration. Standard curves were plotted using concentration as the horizontal coordinate and area as the vertical coordinate. The content of the substance in the samples were calculated using the linear equation of the standard curve. The difference in the amount of each substance in different groups was shown by drawing bar charts. DEMs with criteria of |log2 (FoldChange)|≥1 and p-value<0.05.

The total RNA of the areca exocarp was used to obtain the required mRNA for library construction. Enrichment of mRNA containing polyA tails was performed with Oligo (dT) magnetic beads. To construct the library (NEBNext^® Ultra™ RNA Library Preparation Kit for Illumina^®), the first strand of cDNA was synthesized in the M-MuLV reverse transcriptase system using fragmented mRNA as a template and random oligonucleotides as primers. Subsequently, the RNA strand was degraded with RNaseH and the second strand of cDNA was synthesized with dNTPs in the DNA polymerase I system. Thereafter, the purified double-stranded cDNA was subjected to end repair, A-tailing and ligation to sequencing junctions. The cDNA (250-300 bp) were screened with AMPure XP beads and PCR amplified, and the products were purified again using AMPure XP beads, ultimately resulting in the libraries.

Initial quantification was performed using a Qubit 2.0 fluorometer, and the library was diluted to 1.5 ng/μl. Later, the insert size of the library was detected using an Agilent 2100 Bioanalyzer. After confirming size compliance, the effective concentration of the library was accurately quantified utilizing qRT-PCR to ensure that the quality of the library met the required standards for an effective concentration higher than 2 nM. Sequencing was performed on an illumina NovaSeq 6000 (illumina, USA).

Quantitative analysis was performed using the featureCounts tool in the subread software, based on information about the position of the gene ratio on the reference genome, thus counting the number of reads covered by each gene from the start to the end of the range, and filtering out reads with a ratio quality value of less than 10. Differential analysis of gene expression was performed using DESeq 2 software. GO functional enrichment analysis and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis of differential gene sets using clusterProfiler software. DEGs with criteria of |log2 (FoldChange)|≥1 and padj<0.05.

Correlation analysis was performed by comprehensive screening of exocarp carotenoid-related DEMs and DEGs using quantitative values of genes and metabolites in all samples. The correlation method was to calculate Pearson correlation coefficients of genes and metabolites using cor function in R. Correlation coefficients greater than 0.80 with a pvalue of less than 0.05 correlation result was selected for further discussion and analysis. According to the results of KEGG enrichment analysis, the annotated DEGs and DEMs were simultaneously mapped to the carotenoid-related KEGG pathway.

The total RNA from the areca exocarp was extracted by a RNAprep Pure polysaccharide polyphenol plant total RNA extraction kit (DP441, TianGen, China). Total RNA concentration, quality and integrity were assessed using NanfoDrop 2000 (Thermo, USA) and gel electrophoresis. High-quality total RNA as a template for cDNA synthesis using the PastKing One-step Synthesis of Premixed Reagents by Removing the First Chain of Genomic cDNA (KR118-02, TianGen, China). After the initial screening of DEMs and DEGs, we combined correlation coefficients (>0.8) and KEGG annotations, focusing on the upstream and downstream of the major DEMs, and screened the key genes in the carotenoid biosynthesis pathway as candidate genes for qPCR validation. We used the housekeeping gene (Acactin) as an internal reference gene and the relative expression of the genes was calculated using the relative quantification method (2^-ΔΔCT). The primers involved in qPCR are shown in Supplementary Table 3.

In this study, SPSS 26.0 (IBM Corporation, USA) was used for statistical analysis of data, and Excel 2016 (Microsoft, USA) was used for data collation. GraphPad Prism 9.5 (GraphPad Software, USA), the online bioinformatics platform (http://www.bioinformatics.com.cn) and Adobe Illustrator CS6 (Adobe, USA) software were used for graphical representation.

Measurements of areca nut fruit weight, transverse diameter and longitudinal diameter and coloration are shown in the Figure 2A . And measurements of areca nut exocarp pigment content are shown in the Figure 2B . The single fruit weight, longitudinal diameter, transverse diameter and fruit shape index of the bagged treatments (YP) were 33.38 g, 5.35 cm, 3.60 cm and 1.49, respectively, while several indices of CK were 28.43 g, 5.22 cm, 3.33 cm and 1.57, respectively, and the differences were not significant. However, the results of indicators related to chroma are completely different, and all indicators (L*, a*, b* and CCI) have reached highly significant levels compared with CK (p < 0.01).

We compared the contents of CK and YP pigments. Based on the results shown regarding the pigment content of the exocarp both chlorophyll and carotenoid contents exhibited a sharp decrease after bagging treatment, reaching to highly significant levels (p < 0.01). It is worth mentioning that most pronounced decrease was observed in chlorophyll A content, dropping from an average value of 471.52 μg/g to 1.99 μg/g. The change between YP and CK even reached 237 times. The total chlorophyll content showed a significant reduction (695.43 μg/g vs. 14.41 μg/g), chlorophyll A (471.52 μg/g vs. 1.99 μg/g), chlorophyll B (224.06 μg/g vs. 12.48 μg/g) and carotenoid (405.62 μg/g vs. 24.50μg/g) changing by 48, 236, 18 and 16 times, respectively ( Figure 2B ).

In terms of correlation between indicators ( Figure 2C ), the correlation of single fruit mass with transverse diameter and fruit shape index was highly significant (p < 0.01). Additionally, several indicators, including L*, a*, b* and CCI, all parameters related to chromaticity, as well as the total chlorophyll, chlorophyll A, chlorophyll B, and carotenoids pigments, exhibited highly significant correlation coefficients with each other (p < 0.01).

In this study, we examined 68 carotenoids in areca nut exocarp like α-carotene, β-carotene, and lutein. Furthermore, the principal component analysis (PCA) of two groups, the untreated exocarp and the exocarp of bagged samples are denoted as CK and YP, respectively. The variances of the first and second principal components were 65.48% and 17.25%, respectively. This clearly revealed a significant separation of areca nut exocarp metabolites before and after bagging treatment, with a certain consistency perceived between the three replicates ( Figure 3A ). Hence, the data involved in this study can be used for further analytical exploration.

The results of the screening showed that there were 21 carotenoid-like DEMs in the epidermis of areca nut fruits between the CK and YP, the numbers of up-regulated and down-regulated metabolites were 2 and 19, respectively ( Figure 3B ).

Interestingly, the results show that the vast majority of DEMs were down-regulated after bagging, with all three screened carotenoids (α-carotene, β-carotene, ϵ-carotene) being highly significantly down-regulated (P<0.01) compared to CK. The two up-regulated species (lutein dioleate, violaxanthin palmitate), both belonging to the lutein class, also exhibited highly significantly different compared with CK. Except for lutein dimyristate and lutein distearate, which showed significant differences (P<0.05) from the control. The differences between the YP and CK for the rest of the screened metabolites of lutein class reached highly significant levels (P<0.01) ( Figure 3C ).

Among the screened DEMs, it is noteworthy that lutein had the highest metabolic expression in both CK and YP exocarp when harvested 120 days after pollination. However, concerning the mean value, the expression after bagging was narrowed to 21.5-fold that of CK, with the mean expression value decreasing from 844.62 to 39.20. Although the expression of neoxanthin was second only to that of lutein, the numerical difference between the two substances was still very intuitive. The decrease of neoxanthin after bagging was the most prominent among all DEMs, and the difference between the mean values of the CK group and those of the bagged YP group could reach 45.8-fold, with the mean expression value decreasing from 47.71 to 1.04. In addition, the metabolites such as β-carotene, α-carotene, and violaxanthin were relatively highly expressed in the CK group, with mean expression values of 29.72, 16.07, 9.36, respectively. However, after bagging, the reductions in these examples reached 8.7-fold, 11.1-fold, and 10.4-fold, respectively. Among the rest of the screened DEMs, antheraxanthin decreased by more than 10-fold ( Figure 3D ). Annotation results for the DEMs were categorized by pathway type in KEGG, and KEGG pathway enrichment analysis showed that four pathways (metabolic pathways, biosynthesis of secondary metabolite, carotenoid biosynthesis and biosynthesis of cofactors) were involved in two comparison groups(CK and YP) ( Figure 3E ).

In this study, we aimed to utilize the transcriptomics to comprehend the genetic mechanisms underlying the differences in metabolite content of areca nut exocarp after bagging. As observed from the correlation heat map, which was analyzed for all samples based on transcriptome comparisons, the high data for each sample suggests that transcriptomics data between two groups (CK and YP) were well accurate and reliable ( Figure 4A ). The PCA showed significant sample separation between the YP group and the CK group, with greater variation within the CK group than the YP group ( Figure 4B ). This is consistent with the metabolomics data, suggesting that the accumulation of DEMs after bagging is regulated by the DEGs. The percentage of genomic regions occupied by different instances of the samples in this study is reasonable ( Figure 4C ). And the overall profile of gene expression levels of the samples in this study is presented as a violin plot ( Figure 4D ), with similar interquartile ranges and medians between the CK and YP groups, with the 95% confidence intervals for the latter being relatively more clustered.

Comparing the gene expression revealed that 24814 genes were co-expressed between the CK groups and YP groups, with the number of genes expressed in each of the two groups was 485 and 843, respectively ( Figure 4E ). In addition, volcano plots were generated to visualize the DEGs between two groups. A total of 1784 DEGs were identified from the two groups of samples, 821 genes were up-regulated and 923 down-regulated ( Figure 4F ).

Out of the 1784 DEGs, a total of 713 (40%) were annotated to different transcription factor families (Pkinase, p450 and Pkinase_Tyr etc.).

The 1784 DEGs obtained were further enriched and analyzed using the GO (Gene Ontology) database, and the DEGs were annotated on GO terms in three categories, namely Biological process, Cellular component and Molecular function. We chose the categories with a high percentage of GO secondary entries, in terms of Biological process, the main secondary classifications mainly include metabolic process (GO: 0008152), cellular process (GO: 0009987), localization (GO: 0051179), single-organism process (GO:0044699), and regulation of biological process (GO:0050789) in terms of Biological Process. In terms of Cellular components, they mainly include cell part (GO: 0044464), cell (GO: 0005623), organelle part (GO: 0044422), and membrane (GO: 0016020). In terms of Molecular function the main categories are catalytic activity (GO:0003824), binding (GO:0005488), and transporter activity (GO:0005215). Particularly, metabolic process genes are active ( Figure 5A ).

DEGs were categorized using the KEGG database, and they were annotated on primary (5), secondary (17), and tertiary pathways (40). And the DEGs were mostly enriched in metabolic pathways. Specifically, there was a focus on the carotenoid biosynthesis pathway (ko00906), which belongs to the metabolism of terpenoids and polyketides and is relevant to exocarp color ( Figure 5B ).

Totally, 11 structural genes were annotated for five enzymes in the carotenoid biosynthesis pathway. Subsequently, we conducted further analysis by combining the transcriptomics and metabolomics, correlation analysis was performed using quantitative values of genes and metabolites in all samples. The Pearson correlation coefficients for the genes and metabolites were calculated using the cor function in R, and from these, results greater than 0.80 and with a p-value less than 0.05 were selected. Considering the correlation coefficient, the value of |log2 (FoldChange)|, and the major metabolites, three candidate genes, namely Acat_3g005410, Acat_3g016820 and Acat_9g010750 from the DEGs were selected for further validation.

Based on the previous results and analyses, the main DEMs and candidate genes were simultaneously mapped onto KEGG pathway maps. This approach aimed to provide a clearer understanding of the relationship between metabolites and genes. For better illustration, we mapped the simplified carotenoid biosynthesis pathway ( Figure 6 ).

The correlation between the candidate genes and the main four DEMs (α-carotene, β-carotene, lutein, neoxanthin) was consistently greater than 0.8, with Acat_9g010750 having the highest correlation coefficient, approaching 1 ( Figure 7A ). In addition, we used the NCBI (National Center for Biotechnology Information) database to obtain genes from the palm family (Elaeis guineensis, Phoenix dactylifera) that showed high homology to the amino acid sequences of the candidate genes. We then constructed a phylogenetic tree containing the above genes using MEGA 11 software (Mega Limited, Auckland, New Zealand). Neighbor-joining method employed to infer evolutionary history (Bakker, 2005), and the bootstrap test was conducted with 1000 replicates. Markedly, the three candidate genes belonged to three different taxa ( Figure 7B ).

In our study, candidate genes were validated by qPCR and the results were consistent with phenotypic measurements ( Figure 7C ). This alignment suggests that the analysis of candidate genes was reliable.