Genome-wide identification and co-expression network analysis of Aux/IAA gene family in Salvia miltiorrhiza

Identification and characterization of the Aux/IAA gene family in S. miltiorrhiza

The Aux/IAA protein sequences of A. thaliana and O. sativa were downloaded from The Arabidopsis Information Resources (TAIR) database and Rice Genome Annotation Project (RGAP) database respectively. The genomic data of S. miltiorrhiza was downloaded from National Genomics Data Center (NGDC) database (project number GWHAOSJ00000000).

We employed two methods to identify potential members of the SmIAAs. First, we utilized the A. thaliana Aux/IAA protein sequence as a seed sequence for BLASTP homologous alignment in the S. miltiorrhiza local database. We then validated the domain structure of the resulting sequences using the Pfam database (El-Gebali et al., 2018) and retained only those sequences that contained the AUX_IAA domain as candidate members of the SmIAAs. Second, we utilized the iTAK database to comprehensively evaluate the entire S. miltiorrhiza protein sequence and identified additional members of the candidate SmIAAs (Zheng et al., 2016). By combining the candidate gene family members obtained by these two methods, we determined the final set of members comprising the SmIAAs.

To explore the evolutionary relationships among SmIAAs, we utilized the Neighbor-Joining (NJ) approach implemented in MEGA X software (Kumar et al., 2018), employing the following parameters: 50% partial deletion and 500 bootstrap replications. Additionally, we analyzed the basic physicochemical properties of the protein sequences using the Expasy online website (http://web.expasy.org/protparam/) and predicted their subcellular localization using the WoLF PSORT software (https://wolfpsort.hgc.jp/). To further investigate the chromosomal distribution of SmIAAs, we obtained relevant information from the S. miltiorrhiza genome annotation file and visualized it using the MG2C software (http://mg2c.iask.in/mg2c_v2.0/).

Genetic analysis of Aux/IAA gene sequence and promoter cis-acting elements in S. miltiorrhiza

The structure of SmIAA gene sequence was analyzed by TBtools software according to the GFF3 file and CDS file of S. miltiorrhiza (Chen et al., 2020). The protein sequence domain of SmIAA proteins were mined by CD-search online database (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi). Finally, Tbtools software was used to visually analyze the results (Chen et al., 2020). Furthermore, GO enrichment analysis was performed using R software. In addition, the promoter sequence of 2,000 bp upstream of the start codon of the SmIAAs was extracted from the S. miltiorrhiza genome. Subsequently, the promoter cis-acting elements were analyzed using the online website PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/), with the results of this analysis being visualized through TBtools software, as detailed in Table S1.

Comparative analysis of Aux/IAA gene family across different species

In order to explore the genetic and evolutionary relationships among members of the Aux/IAA gene family in various species, a multiple sequence alignment of 23 SmIAA proteins, 29 AtIAA proteins, and 31 OsIAA proteins was performed using the online program Clustal Omega, followed by MEGAX to create the phylogenetic tree (Table S2), with the parameters used in line with previous analyses (Kumar et al., 2018; Sievers & Higgins, 2014). Furthermore, we utilized the Maximum Likelihood (ML) method to construct an additional evolutionary tree to compare with the results obtained using the NJ method (Fig. S4). This was done with the implementation of specific parameters, including a 50% partial deletion criterion and 500 bootstrap replications. The phylogenetic tree was shown using the online software Evolview (https://www.evolgenius.info/evolview/) (Subramanian et al., 2019). MEME, an online program (https://meme-suite.org/meme/), was used to assess conserved motifs.

Collinearity analysis and selection pressure analysis

The GFF annotation file of the genome was utilized to acquire the chromosomal location data for SmIAAs. The MCScanX software, developed by Wang et al. (2012), was subsequently employed to compute the collinear block and tandem repeat gene information. The TBtools software (Chen et al., 2020) was utilized to visualize the chromosomal distribution and collinear relationships of the SmIAAs. For each identified homologous gene pair by MCScanX, the synonymous substitution rate (Ks) and non-synonymous substitution rate (Ka) were computed using KaKs_Calculator software (Wang et al., 2010) for each synonymous site. The Ka/Ks ratio was utilized to evaluate the rate of protein evolution and determine the selection pressure of the genes during evolution.

RNA-seq analysis of SmIAA expression in S. miltiorrhiza

The RNA-seq data of S. miltiorrhiza were obtained from the NCBI database (PRJNA744957& PRJNA798876). The received sequence read archive (SRA) files were first converted to the fastq format using the fastq-dump program. The downloaded raw data was then trimmed and preprocessed using TrimGalore (https://github.com/FelixKrueger/TrimGalore), and the resulting clean data was utilized for downstream analysis. To construct a genome index with default settings, the Hisat2-build command in Hisat2 (version 2.3.0) (Zhang et al., 2021) was used. The featureCounts software was then employed to calculate the raw counts of each gene in each sample. The standardized transcriptional abundance was determined by computing the transcripts per million (TPM) value. The raw counts of each sample were used to generate a gene expression matrix of 21 transcriptome samples (Table S3). Finally, the expression data for the SmIAAs was extracted and a corresponding heat map was generated using the TBtools software (Chen et al., 2020).

Gene co-expression network analysis in S. miltiorrhiza

Transcription factors were identified in the entire protein sequence of Salvia miltiorrhiza using the iTAK online software developed by Zheng et al. (2016). Gene co-expression network analysis was conducted using R software version 3.9.4 and the Weighted Gene Co-expression Network Analysis (WGCNA) package developed by Langfelder & Horvath (2008). First, the goodSamplesGenes function within the WGCNA package was utilized to identify and filter out any genes with missing values. The soft threshold power of 14 was then selected to construct the weighted gene network, with parameter settings of minModuleSize = 30, mergeCutHeight = 0.25, and maxBlockSize = 6,000. Additionally, to ensure the biological significance of each co-expression module, GO enrichment analysis was performed for each module using R software.

Gene expression analysis of S. miltiorrhiza using qPCR and transcriptomics data integration

Two-year-old S. miltiorrhiza seedlings were collected and frozen in a −80 °C refrigerator. The roots, stems, leaves, and flowers were used for specific expression detection. Fluorescent quantitative specific primers (Table S4) were designed using Primer 5.0, and SmUBQ was used as an internal reference gene. Total RNA was extracted using the RNAprep Pure Plant Plus Kit (Tiangen, Beijing) and reverse transcribed using the PrimeScriptTM reverse transcriptase kit (Takara Tokyo, Japan). Real-time PCR reactions were performed using the PerfectStart Green qPCR SuperMix Kit (Transgene, Beijing), and fluorescence signals were collected using a qTOWER 3G PCR instrument. Each sample was analyzed in triplicate, and relative expression was calculated using the 2^−ΔΔCT method (Table S5). The correlation between transcriptome data and qPCR data was calculated and visualized by Sangerbox website (http://sangerbox.com/) (Shen et al., 2022).

Identification of Aux/IAA gene family members in S. miltiorrhiza

In this study, two distinct methods were employed to identify the Aux/IAA gene family across the entire genome of S. miltiorrhiza. The initial approach involved using Arabidopsis Aux/IAA protein sequences to conduct a BLASTp search on the proteomic data of S. miltiorrhiza. Subsequently, the results were subjected to domain detection via submission to the Pfam database (El-Gebali et al., 2018), and the undesired outcomes were removed. This led to the identification of a total of 23 potential Aux/IAA gene family members. The second method involved employing the iTAK database to analyze the complete protein sequences of S. miltiorrhiza, which also resulted in the identification of 23 candidate genes. By comparing the outcomes of the two methods, it was established that the S. miltiorrhiza genome contained 23 members of the Aux/IAA gene family. In order to simplify future research, the nomenclature of Arabidopsis genes was adopted, and the genes were labeled as SmIAA1-SmIAA23 based on their distribution across the evolutionary tree (Fig. S1).

A comprehensive analysis was conducted to investigate the physicochemical properties of the SmIAA protein sequence. The results revealed that the number of amino acids in the SmIAA protein ranged from 87 to 371 aa, with an average of 211.61 aa. The molecular weight of the SmIAA protein was found to be between 10,364.81 and 40,197.45 Da, while the theoretical isoelectric point ranged from 4.82 to 9.81. Moreover, the subcellular localization analysis indicated that some proteins, namely SmIAA1, SmIAA4, SmIAA5, SmIAA8, SmIAA10, SmIAA11, SmIAA14, SmIAA16, SmIAA18, SmIAA22, and SmIAA23, were located in the cytoplasm, whereas the remaining proteins were all localized in the nucleus (Table 1).

Furthermore, we referred to the genome annotation information of S. miltiorrhiza to identify the distribution patterns of SmIAAs, with the exception of SmIAA23 on contig685. Our analysis revealed that these genes were situated on chromosomes other than chromosome 7, with most genes located at the upper and lower ends of the chromosome. The highest frequency of gene distribution was observed on chromosome 5 (eight genes), while the lowest was on chromosomes 6 and 8 (one gene). Moreover, the genes on chromosomes 2 and 5 were more densely distributed compared to those on other chromosomes, which were more sparsely and randomly distributed (Fig. S2).

Analysis of gene structure and protein sequence

Previous studies have demonstrated that gene structural diversity is a significant driver for the evolution of many gene families (Singh & Jain, 2015). To further explore the gene structure of SmIAAs (Fig. 1A), we conducted an analysis which revealed variable splicing in all genes. The number of introns ranged from 1 to 6. Notably, some genes did not have UTRs. SmIAA2, SmIAA3, SmIAA4, SmIAA6, SmIAA13, SmIAA17, and SmIAA23 contained only CDS and no UTR.

To gain more insights into the properties of SmIAA proteins, we performed domain analysis using the Pfam database. The results showed that all proteins contained Aux/IAA domains, although their placement and size varied (Fig. 1B). We also conducted multiple sequence alignment to evaluate the conserved amino acid sites in the SmIAA proteins. The domains of SmIAA proteins were further divided into Domain I (EAR-motif), Domain II (Degron motif), Domain III (Basic motif), and Domain IV (OPCA-like motif) (Figs. 2 and S3). Some SmIAA proteins were found to be atypical Aux/IAA proteins, lacking at least one conserved domain. This structural diversity may provide these proteins with diverse functions, although further research is required.

Phylogenetic analysis of SmIAA proteins

To investigate the potential relationship between typical and atypical AUX/IAA proteins, we identified and selected the Aux/IAA proteins from two model organisms, A. thaliana and O. sativa, as well as S. miltiorrhiza for further study. We utilized MEGA X to construct phylogenetic trees of the Aux/IAA proteins in three species. The results obtained from the ML method and NJ method were largely consistent, indicating the reliability of our findings. Additionally, we employed the MEME database to analyze the motifs present in these proteins. The resulting evolutionary tree was primarily divided into two groups, I and II, respectively. As illustrated in Fig. 3, the evolutionary tree and conserved motifs were consistent with previous studies on the Aux/IAA family of plants. The majority of the members were classified in Group I, with nearly all of them being canonical Aux/IAA proteins possessing four canonical domains (motif1, 2, 3, and 4). On the other hand, nearly all the members of Group II were found to be atypical Aux/IAA proteins, lacking at least one conventional domain. Atypical Aux/IAA proteins have been found to be prevalent in the Aux/IAA gene family of various plants, where they play a significant role in plant adaptability to changing environmental conditions, although their specific functions remain largely unknown (Luo, Zhou & Zhang, 2018).

Analysis of collinearity and selection pressure

We conducted a collinear analysis of SmIAAs to elucidate their relationships further. The results showed that nine lineal homologous gene pairs were identified among all SmIAAs, suggesting that chromosomal replication processes play a crucial role in the amplification of the S. miltiorrhiza Aux/IAA gene family members (Fig. 4).

To gain a better understanding of the evolutionary processes of genes and proteins, studying selective interaction patterns can be a useful approach. Selection has been shown to be a valuable technique for speculating about gene function (Wang et al., 2009). The Ka/Ks value, which is the ratio of the two protein-coding genes’ non-synonymous substitution rates (Ka) and synonymous substitution rates (Ks), can be used to identify whether the gene is under selection pressure. When Ka/Ks > 1, the gene is thought to be positively selected; when Ka/Ks = 1, the gene is thought to be purified and selected; and when Ka/Ks < 1, the gene is thought to be neutrally selected. Except for SmIAA8/SmIAA14 and SmIAA8/SmIAA19, the Ka/Ks value of the other S. miltiorrhiza lineal homologous gene pairs is less than 1, which means that these SmIAAs evolved under purification selection, indicating that the sequence of the Aux/IAA gene is relatively conservative. These findings suggest that the Aux/IAA gene plays a significant role in S. miltiorrhiza and attains an optimal state, which is consistent with the conservative function of several SmIAAs.

Analysis of GO enrichment and cis-acting element

The GO categorization system is a widely used method for characterizing genes and proteins in various organisms. The system is organized into three distinct categories, including molecular function, biological process, and cellular component, with each category corresponding to specific functions or properties. To investigate the potential functions of the SmIAAs, we utilized TBtools to examine the GO enrichment of SmIAA and identified the top 10 most significantly enriched GO terms across the three categories (if less than 10, all).

Figure 5A demonstrates that the SmIAAs are associated with 10, five, and three GO terms related to biological process, cellular component, and molecular function, respectively. The majority of SmIAAs are believed to play a role in various biological activities. Their molecular functions include transcription regulator activity, DNA-binding transcription factor activity, and identical protein binding. The cell composition mainly includes the nucleus, nuclear, intracellular membrane-bounded organelle, membrane-bounded organelle, and other cellular components. Furthermore, SmIAAs are primarily involved in biological processes related to the response to auxin, organic substance, endogenous stimulus, and hormone, among others.

To gain a more comprehensive understanding of the role of SmIAAs in the growth and development of S. miltiorrhiza, we analyzed the cis-acting elements in the promoter region of the related genes, which is 2,000 bp upstream of the initiation codon (Fig. 5B and Table S1). Our analysis revealed the presence of numerous cis-acting elements associated with auxin, salicylic acid, abscisic acid, methyl jasmonate, and gibberellin in the promoter region of SmIAAs, in addition to key cis-acting elements associated with transcription and light response. The presence of these cis-regulatory components suggests that hormones have the capacity to regulate SmIAAs. Overall, our analysis of cis-acting elements reveals that different family members contain a variety of response motifs, implying that the SmIAA gene family may participate in various physiological processes of S. miltiorrhiza by responding to various response elements such as hormone, low temperature, light, and so on.

Analysis of gene expression pattern

The temporal and spatial expression of a gene may provide insights into its functional features (Wang et al., 2021). In this study, transcriptome data retrieved from NCBI was used to extract the TPM values of root, stem, leaf, and flower tissues of S. miltiorrhiza, and subsequently analyze the expression pattern of the SmIAAs. The results of gene expression patterns based on transcriptome data indicated considerable variability in the expression of several SmIAAs across different tissues of S. miltiorrhiza, as well as in the expression trends of individual genes. Based on the results of expression cluster analysis, SmIAAs could be categorized into six groups (A~F) (Fig. 6). Group A consisted of two genes (SmIAA4 and SmIAA6) that were primarily expressed in the leaves of S. miltiorrhiza. Group B included four genes (SmIAA19, SmIAA16, SmIAA10, and SmIAA11) that showed similar expression levels in the leaves, stems, and flowers of S. miltiorrhiza, and relatively low expression levels in the roots. Group C comprised five genes (SmIAA14, SmIAA1, SmIAA5, SmIAA23, and SmIAA22) that were predominantly expressed in the flowers and roots of S. miltiorrhiza. Group D encompassed six genes (SmIAA12, SmIAA18, SmIAA21, SmIAA7, SmIAA17, and SmIAA20) that showed the highest expression levels in the stems, roots, and flowers of S. miltiorrhiza, but minimal expression in the leaves. Group E consisted of three genes (SmIAA13, SmIAA15, and SmIAA8) that were primarily expressed in the stems and roots of S. miltiorrhiza. However, the three genes in group E (SmIAA2, SmIAA9, and SmIAA3) were significantly expressed only in the stem of S. miltiorrhiza. In summary, SmIAAs were expressed in clusters in various tissues of S. miltiorrhiza, with the majority of SmIAAs exhibiting high expression levels in roots, stems, and particularly flowers.

Weighted correlation network analysis

Although the Aux/IAA gene family is crucial for auxin signal transduction, prior studies have suggested that its members may interact with other transcription factors (TFs) besides ARFs (Luo, Zhou & Zhang, 2018). To explore these interactions, we employed the iTAK tool to identify all TFs in the S. miltiorrhiza genome and utilized weighted correlation network analysis (WGCNA) to examine additional transcriptome data from S. miltiorrhiza flower tissue at different stages (withering, full blooming, and bud). Our analysis identified 68 TFs, including ARF, ERF, MYB, and others (Tables 2 and S6). We used a soft threshold of β = 14 to divide the 1814 TFs into 10 modules, with gene numbers ranging from 45 to 490 in each module (Table S7). The outcomes of the GO enrichment analysis of distinct modules revealed that each of these modules has distinct biological significance (Table S8). Furthermore, utilizing the Pearson correlation data, we observed a significant association between SmIAAs and the modules, suggesting that SmIAAs may be linked with the module hub genes (Fig. 7). We selected the red and cyan modules for further investigation, as they were highly associated with SmIAA14 and negatively related to SmIAA22, respectively. We identified the top 10 hub genes in each module and annotated their homologs in the iTAK database (Table 3). Surprisingly, only one ARF was among the 20 key genes identified across various modules (EVM0019866), suggesting that SmIAAs may have other important functions beyond their role in the auxin pathway. However, more research is needed to confirm this hypothesis.

The validation of qRT-PCR

To confirm the RNA-Seq results, we employed the quantitative reverse transcription polymerase chain reaction (qRT-PCR) technique to evaluate the expression levels of 12 genes, comprising of four SmIAAs, four hub genes selected by co-expression modules, and four transcription factors chosen randomly (Fig. 8). The findings demonstrated that the alteration in gene expression detected by qRT-PCR and RNA-Seq techniques was similar, and the outcomes of the two datasets were positively correlated, suggesting that the co-expression network established in this study is reliable.