We downloaded data containing sequence IDs, protein sequences, genomic sequences, and conserved domain database (CDD) information from the A. annua genome database (the NCBI accession number: PKPP00000000). The characterized AtHAT and OsHAT sequences were obtained from the Phytozome database (https://phytozome.jgi.doe.gov/pz/portal.html V12.1) to identify all HATs in the A. annua genome. Eleven putative AaHATs were characterized via BLASTP searches (E< 10^-10) against the NCBI database. Next, the hidden Markov model (HMM) profiles of HAG (PF00583), HAC (PF08214), HAF (PF09247), and HAM (PF01853) were downloaded from the Pfam database. We applied the same strategy to search for conserved domains in the 11 AaHATs and ensure that they belonged to the HAT family, and TBtools software was used for visualization. Detailed information on the AaHATs was downloaded in batches from the NCBI. Information on the AaHAT gene family, which contains gene IDs, protein IDs, location, CDS and protein length, protein molecular weight (MW), isoelectric point (pI), and numbers of exon and groups is provided in detail ( Supplementary Table S1 ).

Multiple sequences were aligned to AtHAT, OsHAT, and AaHAT, and the protein sequences were analyzed using MEGA7.1 (https://www.megasoftware.net/) to construct an unrooted phylogenetic tree via the neighbor joining method with 1000 bootstrap replicates.

The genome annotation of A. annua was obtained from the NCBI database. Gene structure analysis was conducted using TBtools-II (https://github.com/CJ-Chen/TBtools/releases) with the input of the A. annua GTF file. We used SWISS-MODEL to predict the structures of the AaHATs. The online tool multiple Em for motif elicitation was applied to analyze protein motifs.

We defined upstream 2Kb sequences counted from the ATG start codon as the promoter regions with the help of TBtools-II. These promoter sequences were then used to identify cis-elements in the predicted promoter in the website of PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/). The cis-elements were defined and grouped with their various functions. Cis-elements of abscise acid responsiveness (ABRE), gibberellin-responsive elements (P-box elements), and auxin responsive elements (TGA-element) were identified in the promoter regions of the 11 AaHAT.

Data of the transcription pattern of AaHAT genes in six tissues of A.annua were downloaded from a published NCBI database (https://www.ncbi.nlm.nih.gov/) (SRA accession number: SRP129502). These tissues included young and old leaves, stem, seeds, roots, and flowers. All RNA-seq library preparation was conducted as described in previous study. Heatmap Illustrator HemI 2.0 (http://hemi.biocuckoo.org/) was used to plot the heatmap of gene transcription.

Plantlets were cultivated in growth chambers under optimal conditions of temperature (25 ± 2°C) and light provided by cool white fluorescent tubes emitting a light intensity of 40 μmol m⁻² s⁻¹. All plantlets were subjected to a 16-hour photoperiod. Mass propagation of plantlets was achieved through sub-culture in MS medium every 15 days. To investigate the responses of AaHATs to various abiotic stresses, 2-week-old plants grown on plates with moist gauze were treated with a solution containing 200mM NaCl and 100 μM ABA. For cold treatments, the plants were exposed to 4°C. Leaves from all samples for RNA extraction were collected at 0 h, 1 h, 2 h, 5 h, and 24 h post-treatments.

Total RNA for quantitative real-time PCR (qRT-PCR) analysis was extracted using Trizol reagent (Invitrogen) following the manufacturer’s protocol. Subsequently, cDNA was generated from 2 μg of total RNA utilizing the EasyScript One-Step gDNA Removal and cDNA Synthesis Supermix (Transgene Biotech, China). The qRT-PCR analysis was conducted with TranStart qPCR SupMix chip (Transgene Biotech, China) on the ABI StepOnePlus real-time PCR system with a minimum of three biological replicates. Data analysis was performed using the comparative Ct (ΔΔCt) method in Excel. Statistical significance between treatment and control groups was determined using Student’s t-test, with p values below 0.05 or 0.01 considered significant.

We identified and characterized 11 HATs in A. annua with A. thaliana and rice HATs as queries ( Supplementary Table S2 ). Consistent with their evolutionarily conserved sequences and the identified HATs in A. thaliana, the 11 AaHATs were classified into five groups with short names: HAC, HAG1, HAG3, HAF, and HAM. Every group had individual conserved domains, which reinforced the relevance of the classification ( Figure 1 ). Approximately 60% of the AaHAC CDSs were 800–2000 aa long, while PWA56962.1 (385 aa), PWA64079.1 (561 aa), PWA77834.1 (521 aa), and PWA80289.1 (561 aa) were unique. The protein sequences of each group displayed substantial similarity. The molecular weights of the AaHATs ranged from 45.31 to 207.58 kDa. The isoelectric points varied from 5.188 to 8.443. PWA89578.1 encoded the largest protein with a molecular weight of 207.58 kDa, but PWA56962.1 encoded the smallest protein with a molecular weight of 45.31 kDa ( Supplementary Table S1 ). The protein nature of the AaHATs is highly similar to that of HATs from other plants (Pandey et al., 2002; Liu et al., 2012), which indicates that the functions of these AaHATs are evolutionarily conserved.

To probe the phylogenetic relationships of plant AaHATs, we generated twelve A. thaliana, eight O. sativa, and eleven A. annua HAT amino acid sequences which were used to construct a neighbor joining tree. There were no assumptions about ancestral representation based only on the relationships of the leaf nodes in unrooted trees. As a consequence, HATs from the three plant species, A. thaliana, rice, and A. annua, were divided into six groups, as expected ( Figure 2 ). The AaHAT proteins exhibited strong similarity with HATs from module plants. These genes were classified into clades with AtHATs and OsHATs, which had better bootstrap values. One AtHAT, one OsHAT, and two AaHATs were grouped into the HAG1 clade, while for the HAG3 group, each species had one HAT ( Figure 2 ). Intriguingly, there was no AaHAT in the HAG2 group ( Figure 2 ). These data indicated that the HAT in the A. annua genome has unique biological evolution compared with that of these module plants. Regardless of the plant species, the HAC group had the largest group, with five AtHATs, three OsHATs, and six AaHATs ( Figure 2 ). Two AtHATs and OsHATs were clustered into HAM and HAF ( Figure 2 ). Overall, we found that HATs in the A. annua genome are highly conserved in terms of evolution and are rich in highly variable regions; thus, HATs have become comparatively ideal objects of medical plant identification and systematic study.

To visualize the various protein structures of histone acetylation transferases, we randomly selected HATs from each group of three plant species for analysis ( Supplementary Figure S1 ). HAM proteins have similar structures among diverse plant species ( Supplementary Figure S1 ). In contrast, the HAC, HAF, and HAG protein structures appeared to diverge from each other in this study ( Supplementary Figure S1 ). In addition, we found that the conserved sequence of the HAT protein was complete, but only partial folding regions differed slightly ( Supplementary Figure S1 ). AaHACs had one special gene (PWA64079.1) that evolutionarily lost regional introns. Consequently, the 3D structures of these genes are similar to those of other family groupers of AaHAC ( Figure 1 and Supplementary Figure S2 ). In total, for the same groups of HATs, we found that these proteins from different plant species possess similar structures, and vice versa. Given that protein structure and character determine the properties and functions of organisms, these data suggest that HATs from these three plant species may have homologous biological functions.

To better understand gene family evolution, we compared the structures of the AaHATs ( Figure 3A ). Analysis of the genomic DNA sequences revealed that the number of introns varied from 8 to 20 ( Figure 3A ). AaHATs mostly shared highly similar structures, which were grouped together in three branches of the NJ tree ( Figure 3A ). The majority of the AaHACs had similar numbers of introns and exons, with the exception of one AaHAC gene (PWA64079.1). AaHAF and AaHAG3 had the greatest number of introns (20) among all the groups. Moreover, to identify putative motifs in the HAT family in A. annua, we predicted the protein sequences of the 11 AaHATs using the MEME website. Consequently, 20 motifs were identified from these proteins ( Figure 3B ; Supplementary Table S3 ). The group contained similar motifs, which indicated that the proteins might have some common functions. HAC genes are evolutionarily conserved in other plants and possessed the greatest number of motifs. There was at least one motif in each AaHAT, and they were arranged in the same order in each group except for AaHAM (PWA56962.1), which had no predicted motif. Overall, this analysis supported the presumption that motifs commonly present in histone acetyltransferases in the A. annua genome may be related to the conserved functions of AaHATs, whereas the AaHATs seem to be associated with specific functional rules.

Given that subcellular location details can provide some clues for the prediction of protein function, we found that all of the AaHATs were predicted to be localized in the nucleus with high reliability (RI > 6), except for HAM (PWA56962.1), which might be localized in the chloroplast, mitochondria, and nucleus ( Supplementary Table S4 ). In particular, HAC (PWA69804.1) was located only in the nucleus with high reliability (RI = 14), and other AaHATs were predicted to be located in at least two subcellular organelles ( Supplementary Table S4 ). Similarly, Uniprot identified possible nuclear signals for AaHATs, except for HAM (PWA56962.1) and HAG3 (PWA79021.1), which were not predicted. The identified AaHAT proteins in the A. annua genome exhibited different subcellular distributions and may be involved in medical plant developmental growth.

To determine the transcriptional levels of the AaHAT genes in different organs, we generated a heatmap with the RNA-seq data from the A. annua database ( Supplementary Table S4 ). We detected complex, specific, and overlapping AaHATs in different tissues. The tissues/organs were divided into six types: flowers, young leaves, old leaves, seeds, stems, and roots. According to the results, 11 AaHAT genes were distinctly expressed at different developmental stages in the tissues. The transcriptional level of the genes was dynamic among these tissues and organs. For example, PWA92522.1 (HAC) was barely detectable in whole plants; in contrast, transcriptional signals of the PWA53508.1 (HAC) gene were highly detected in the overall plantlet ( Figure 4 ). Additionally, the transcript levels varied within the same organs. In the six A. annua tissues, the expression of the members of the AaHAC gene family was relatively low ( Figure 4 ). Some genes were exclusively transcribed in single tissues or organs; for example, PWA69804.1 (HAC) had a high transcriptional level in old leaves. The transcriptional levels of the AaHATs varied among tissues and were involved in the development of A. annua.

To better understand the potential roles of AaHATs, we used a plant promoter database to identify the region 2 kb upstream of the transcription initiation site of AaHATs. Overall, we identified 619 cis-acting elements within the promoter regions of the AaHATs. The associated elements included stress-related, hormone-responsive, light-responsive, developmental, promoter and enhancer, site-binding-related, and other elements ( Figure 5A ). Compared with those of other elements, the abundance of light-responsive and promoter and enhancer elements increased ( Figure 5A ). In addition, cis-elements in the promoter regions of the AaHATs were found to be relevant to plant hormone classifications ( Figure 5B ). As a consequence, all 11 AaHATs contained abscisic acid responsive-elements (ABREs), while gibberellin-responsive (P-box elements) and auxin responsive (TGA-element) elements were identified in practically all A. annua HAT promoters ( Figure 5B ). Interestingly, obvious cis-elements between the paired genes, including gibberellin and MeJA response elements, were identified in the promoter of AaHAG1 (PWA77834.1) ( Figure 5B ). Notably, no cytokine-responsive elements were found in these promoters ( Figure 5B ), which was consistent with previous results (Xing et al., 2022). These results showed that AaHATs may affect phytohormone responsiveness, adaptability to stress, and developmental growth.

It is important to study gene function to construct gene family interaction networks. Here, we constructed a regulatory network controlled by A. annua HAT genes with STRING. Among the 11 AaHATs, 9 were identified to interact with proteins with high confidence (score > 0.8). AaHAC, which has four groups, namely PWA88299.1, PWA53508.1, PWA88299.1, and PWA64079.1; AaHAG1, which includes PWA80289.1 and PWA77834.1; PWA89578.1 (HAF); PWA56962.1 (HAM); and PWA79021.1 (HAG3) constitute an independent regulatory network with 34 proteins involved in crosstalk ( Supplementary Table S6 ).

In detail, we determined the PCC values corresponding to the A. annua transcriptome of dynamic organs in the network. In our analysis, the red rectangles represent AaHATs, and the red edges represent the PCC values of AaHATs and interaction scores > 0.5 ( Supplementary Figure S3 ). Among the nine AaHATs, nine expressed transcripts had PCCs > 0.5 for PWA88299.1 (AaHAC), which was highly expressed in young flowers and young leaves; eight genes had PCC values greater than 0.5 for PWA53508.1 (AaHAC), which was most highly expressed in young/old leaves ( Supplementary Figure S3 ). Five genes had PCC values > 0.5 for PWA64079.1 (AaHAC), which was preferentially expressed in young leaves, and 2 genes had PCCs > 0.5 for PWA88299.1 (AaHAC), which also exhibited high transcription in young leaves ( Supplementary Figure S3 ). Four genes had PCC values greater than 0.5 for PWA77834.1 (AaHAG1), and these genes were highly expressed in young leaves ( Supplementary Figure S2 ). The 3 genes had PCCs >0.5 for PWA80289.1 (AaHAG1), which was most highly expressed in old leaves, stems, and roots, and 2 genes had PCC values > 0.5 for PWA79021.1 (HAG3), which was preferentially expressed in seeds ( Supplementary Figure S3 ). Four genes had PCC values > 0.5 for PWA89578.1 (AaHAF) and were highly expressed in the seeds ( Supplementary Figure S3 ). The four genes had PCCs > 0.5 for PWA56962.1 (AaHAM) and were highly expressed in young leaves and roots ( Supplementary Figure S3 ). Interestingly, PWA94397.1, which is similar to the Arabidopsis and rice ARBE (abscisic acid responsive element) binding factors, was one of the genes expressed with AaHAG1 (PWA80289.1), with a PCC greater than 0.5 ( Supplementary Figure S3 ; Supplementary Table S6 ). Furthermore, the entire promoter of the predicted AaHATs contained ABRE cis-elements ( Figure 5 ). Thus, we speculate that AaHAG1 may interact with other binding factors of the AaHAT gene, which indicates that HATs in A. annua may have self-regulating systems.

HATs may be involved in plant responses to environmental stresses. In this study, we investigated the expression levels of AaHATs under three treatments, cold, NaCl, and ABA, via qRT-PCR. We found that the expression of the AaHAT genes increased at most cold treatment time points. In addition, AaHACs were induced quickly with NaCl treatment and decreased after 5 hours of NaCl exposure in comparison with those in the control group ( Figure 6 ). Interestingly, the expression of orthologs of these genes in rice increased after NaCl treatment, which implies the evolutionary conservation of HAT genes in the plant kingdom. Interestingly, the expression profiles of the AaHAC proteins exhibited obvious positive correlations with all PCCs > 0.6 in both our tested organs and stress treatments, suggesting functional redundancy of the AaHAC proteins in A. annua ( Supplementary Figure S3 ). In addition, we also found that some associated genes in the artemisinin biosynthetic pathway transcriptionally activate under stress. In detail, in terms of cold stress, the expression levels of the four artemisinin biosynthetic genes were also significantly increased, such as AaADS, AaCYP71AV1, and AaDBR2 ( Supplementary Figure S4A ). As for NaCl treatment, the expression of AaADS and AaDBR2 were induced ( Supplementary Figure S4B ). After ABA application, we found that only AaCYP71AV1 was active ( Supplementary Figure S4C ). The results reinforced the assumption that the HAC proteins in monocot plants have shared ancestral genes, which may activate artemisinin biosynthesis related genes by ectopic acetylation under cold and NaCl treatment. However, it was not obvious that AaHATs were induced under ABA treatment, including the downstream genes ( Figure 6 ). In summary, it was indicated various AaHATs had special roles under various environment stimuli, which may activate artemisinin biosynthesis associated genes to promote the plant-derived output.