Transcriptome-based analysis of the hormone regulation mechanism of gender differentiation in Juglans mandshurica Maxim

Plant material

Based on the previous research results of the phenology, morphology, physiology, and cytology of J. mandshurica natural forest in eastern Liaoning, China, the morphological development period, physiological differentiation period, and vegetative growth period of J. mandshurica flower bud differentiation process were determined. In 2019, the flower buds were selected on March 17 and 25, and April 2, 10, and 15 to be quick-frozen in liquid nitrogen and stored in a refrigerator at –80 °C for differential gene expression analysis and metabolomics. Simultaneously, the female and male flower buds collected on April 2 during the physiological differentiation period of the male and female precursors (Fig. 1) were selected for transcriptomic research. The female flower buds of the female precursor type were named T1, which were repeated three times, namely A1, A2, and A3, respectively. The male flower buds were named T2, which were repeated three times, namely A4, A5, and A6. The female flower buds of the male precursor type were named M1, which were repeated three times, namely B1, B2, and B3. The male flower buds were named M2, which were repeated three times, namely B4, B5, and B6. There were 12 samples in total.

Construction and sequencing of library

Oligo(dT) magnetic beads were used to enrich mRNA with polyA structure in total RNA, and the RNA was interrupted to 200–300 bp fragments by ion interruption. After the construction of the library was completed, PCR amplification was adopted to enrich the library fragments, and then library selection was carried out according to the fragment size, where the library size is 300–400 bp. Next, the library was inspected using an Agilent 2100 Bioanalyzer, and the total concentration and effective concentration of the library were detected. After RNA extraction, purification, and library building, these libraries were sequenced using the second next-generation sequencing technology based on an Illumina HiSeq sequencing platform. The sequencing raw read length was 150 bp, and the original sequence of mRNA in the physiological differentiation stage of male and female flower buds of J. mandshurica was obtained.

Sequencing data sorting and transcript assembly

We counted the Raw Data of each sample separately, spliced Clean Reads with Trinity software to obtain transcripts, and then carried out the subsequent analysis. Trinity was used for D. novo assembly software for transcriptome splicing, which splices high-quality sequences based on the De Bruijn Graph splicing principle. After splicing was completed, the transcript sequence file in the FASTA format was obtained. The longest transcript under each gene was extracted as the representative sequence of the gene, called unigene.

Unigene annotation classification and expression analysis

The database used for gene function annotation included NR (NCBI non-redundant protein sequences), GO (Gene Ontology), KEGG (Kyoto Encyclopedia of Genes and Genome), eggnog (evolutionary genealogy of genes: Non-supervised Orthologous Groups), and Swiss-Prot.

Differential expression and enrichment analysis

By comparing the RNA-seq data of the physiological differentiation stage of female and male flower buds, the gene expression was analyzed by DESeq. The conditions for screening differentially expressed genes were as follows: expression difference multiple |log2FoldChange| > 1, significant P-value < 0.05. Significant enrichment analyses were conducted of differentiated genes based on the number of differentially expressed genes contained at different levels of each KEGG pathway, after multiple tests and corrections, the pathway with a Q-value ≤ 0.05 was selected as the pathway significantly enriched in differentially expressed genes, and then the metabolic pathway and signal pathway in which differentially expressed genes were mainly involved were determined.

qRT-PCR verifies sequencing data

Twelve differential expression genes were randomly selected, and their expression profiles were detected via qRT-PCR to verify transcriptome data.

Endogenous hormone determination

Samples from different periods were tested for phytohormones, including IAA, GA, CTK, ABA, ETH, BR, JA, and SA, then analyzed and compared (Weiler, Jourdan & Conrad, 1981).

Transcriptome sequencing analysis of male and female flower buds

Raw transcriptome sequencing reads all flower bud samples with raw sequencing reads of more than 40 M, an average reading length of 150 bp, a total number of sequencing bases of over 6G nt, a Q20 content of over 95%, and a Q30 content of over 90%. Statistical analysis revealed that the sequencing data were of high quality and could meet the requirements of further analysis. The data filtering used Cutadapt to remove the 3′-end linker, and the removed part had at least 10 bp overlap (AGATCGGAAG) with the known linker, allowing a 20% base mismatch. Remove read with an average mass fraction lower than Q20. The number of clean reads obtained from female flower buds of female precursor type ranged from 43,894,726 to 46,469,762. The number of high-quality sequence reads was 6,610,627,884 to 6,991,929,010, the percentage of high-quality sequence reads in sequencing reads was 99.41–99.49%, and the percentage of high-quality sequence bases in sequencing bases was 99.13–99.16% (Table S1).

Using the reads assembly software Trinity to assemble from scratch, the assembly fragments contig were obtained according to the overlap of reads, and these contigs were linked together by overlap to obtain unigene. Transcript and unigene sequences were 344,218 and 138,138, respectively, and the transcript and unigene sequences had average lengths of 1,182 bp and 788 bp, respectively. The number of unigene in the physiological differentiation stage of male and female flower buds of J. mandshurica was approximately 80,000, while that of unigene larger than 500 bp was almost 60% of the total. This indicated that the sequencing results were of high splicing quality, which was conducive to further functional annotation analysis (Fig. 2).

All-unigene obtained by transcriptome sequencing was annotated with protein and nucleic acid functions (Table 1). Among the 108,828,259 bp unigene assembled in total, all-unigene sequences were sequentially aligned to the protein database through blastx, and unigene was annotated with gene functions. The databases used for gene function annotation included NR (NCBI Non-Redundant Proteins Sequences), GO (Gene Ontology), KEGG (Kyoto Encyclopedia of Genes and Genome), eggNOG (evolutionary genealogy of genes: Non-supervised Orthologous Groups), and Swiss-Prot. The number of successful annotations of unigene into the Nr protein database reached 60,415, accounting for 43.74% of the total number of all-unigene. The number of unigenes annotated in the GO database reached 29,402, accounting for 21.28% of the total unigenes. The number of unigenes annotated in the eggNOG database reached 57,250, accounting for 41.44% of the total unigenes. The number of unigenes annotated in the Swissprot database reached 45,347, accounting for 32.83% of the total unigenes, while the number of unigenes annotated in the KEGG database was at least 6,845, accounting for only 4.96% of the total. The number of unigenes annotated in all databases was 5,312, accounting for only 3.85% of the total unigenes.

According to the statistics of unigenes annotated in the Nr protein database by E-value distribution, e-values in 1e⁻⁵–1e⁻¹⁵ accounted for 26.72%, and 1e⁻¹⁵–1e⁻³⁰ accounted for 22.33% (Fig. S1A). According to the annotation similarity distribution statistics, 19.45% of the sequences were concentrated in the 80–90% range, with a total of 55.51% sequence similarity higher than 60% (Fig. S1B). Homology analysis among different species demonstrated that Vitis vinifera had the highest homology, reaching 19.26%, followed by Ziziphus jujuba (6.29%), and other species (58.27%). The distribution of species homology revealed high homology with grape (Fig. S1C).

The gene was annotated by eggNOG comparison, and the eggNOG number of the best comparison result was assigned to the corresponding gene. Further, using the correspondence between the eggNOG number and the eggNOG classification catalog, each gene was classified into the eggNOG classification catalog. A total of 77,173 unigenes were compared to the eggNOG database (Fig. 3). Among them, 21,187 were compared with the predicted general function prediction only, and 15,111 were compared with the predicted unknown function. Then, compared with the predicted number of aggregation classifications, the sequence was 4,728 signal transmission mechanisms, 4,354 transcriptions, 4,069 replications, recombination, and repair. Post-translational modification, protein transformation, 4,042 molecular chaperones, 2,896 translations, ribosome structure, and biogenesis. The above approach demonstrated that substance and energy metabolism were very vigorous in the process of J. mandshurica flower bud differentiation, cell growth, proliferation, and division speed were very fast, and hormone signal transduction participated in the differentiation process.

GO annotation was completed using BLAST2 GO software, and the annotation adopted the default parameters of BLAST2 GO. The results of the GO annotation were mapped to go term, and the number of genes annotated in the second-level classification was counted (Fig. S2). According to classification statistics, there were 29,402 (21.8%) unigene genes annotated in the J. mandshurica transcription group into the GO database. There were three major categories of annotations: 32,959 (24.5%) of molecular function, 51,672 (38.4%) of biological processes, and 50,032 (37.2%) cellular components. There were 67 functional groups. Among the 24 functional groups of biological processes, unigene included 14,613 metabolic processes (GO:0008152) and 14,514 cellular processes (GO:0009987). There were 10,152 cells (GO:005623) and 10,063 cells part (GO:0044464). Among the 19 functional groups of molecular functions, there were 14,843 bindings (GO:0005488) and 14,061 catalytic activities (GO:0003824). Through GO function annotation, a large number of unigenes related to metabolic process, biological regulation, cell process, cell, membrane, and membrane component organelles, and other unigenes were obtained in the three major types of functions. These genes can provide a theoretical basis for the study of sex differentiation and flowering related functional genes of J. mandshurica.

A total of 6,288 unigenes were mapped to 35 KEGG metabolic pathways, of which 606 genes were annotated by the signal transduction pathway in environmental information processing. A total of 504 genes were annotated by the translation pathway in genetic information processing, and 487 genes were annotated by carbohydrate metabolism in the metabolic pathway (Fig. S3).

Analysis of GO and KEGG pathway in male and female flower buds

According to the classification and annotation of the GO function of differentially expressed genes in male and female flower buds of female precursor type (T1 vs T2), we can observe that there were 837 functions, including 230 molecular functions, 517 biological processes, and 90 cell compositions. We analyzed the significance of GO function enrichment of female and male flower buds of female precursor type (when Bonfeeoni-corrected P-value ≤ 0.05). The results revealed that T1 vs. T2 significantly enriched 36 terms, 18 terms belonged to the biological process, eight terms belonged to molecular function, and 10 terms belonged to cellular component (Fig. 4). The enriched term included nucleic acid binding transcription factor activity (GO:0001071), regulation of RNA biosynthesis processes (GO:2001141), ribosomal subunit (GO:0044391), etc. The cluster frequencies are 6.1%, 6.7%, and 9.9%, respectively. The main annotated genes include GAI, SOC1, AP1, etc. The metabolism of endogenous substances changed during the physiological differentiation of female and male flower buds of J. mandshurica, and, simultaneously, the genes related to flower formation were pre-activated to a certain extent, laying a foundation for the differentiation of flower organs. There were 121 genes annotated on the pathway compared in T1 vs. T2. Differential Expression Genes (DEGs) significantly enriched 20 pathways in KEGG (Fig. 5), of which plant hormone signal transduction was the highest (ko04075). The second was carbon metabolism (ko01200). Additionally, glyoxylate and dicarboxylate metabolism (ko00630), AMPK signaling pathway (ko04152), and zeatin biosynthesis (ko00908) were also significantly enriched. These results demonstrate that hormones play a leading role in regulating the differentiation of female and male flower buds of the female precursor type, and carbohydrate was also an essential metabolic substance during the transformation from physiological differentiation of flower buds to morphological differentiation.

Male and female flower buds of male precursor type (M1 vs. M2) had a total of 1,013 functions annotated, including 649 biological processes, 94 cell compositions, and 270 molecular functions (Fig. 6). A total of 54 terms were significantly enriched, included 28 belonged to biological process, 21 belonged to molecular function and five belonged to cellular components. The enriched term included protein dimerization activity (GO:0046983), regulation of transcription, DNA-templated (GO:0006355), extracellular region (GO:0005576), etc. The cluster frequencies are 5.6%, 9.8%, and 4.0%, respectively. The main annotated genes include ARF5, GA2OX1, TEM1, SOC1, etc. A total of 125 differentially expressed genes DEGs of M1 vs. M2 were enriched in the KEGG pathway (Fig. 7). Among them, the ribosome (ko03010) was the most abundant, followed by plant hormone signal transmission (ko04075) and carbon metabolism (ko01200), which indicated that M1 vs. M2 was affected by gene regulation, hormone regulation, and physiological metabolites in the process of sex differentiation. The KEGG significant enrichment analysis (Bonferroni corrected P-value ≤ 0.05) revealed that the KEGG of M1 vs. M2 significantly enriched the 18 pathways. Among them, plant hormone signal transduction (ko04075) was the highest, followed by the neurotrophin signaling pathway (ko04722) and starch and sucrose metabolism (ko00500). This indicates that hormones play a leading role in the physiological differentiation of male and female flower buds of male precursor type and regulated the growth and differentiation of organs in all parts of the flower.

Analysis of DEGs

According to the screening of differentially expressed genes in the flower bud transcriptome of J. mandshurica, we obtained the DEGs comparison library of female and male flower buds: male pre-type female flower buds and male flower buds (M1 vs. M2, 1,357 genes), male pre-type female flower buds and female pre-type female flower buds (M1 vs. T1, 1,461 genes), male pre-male flower buds and female pre-type male flower buds (M2 vs. T2, 3,476 genes), female pre-type female flower buds, and male flower buds (T1 vs. T2, 1,822 genes) (Fig. 8A). This indicated that most of the DEGs appeared in flower buds of different genders, and 11 of them were differentially expressed in the four libraries. A total of 8,116 differentially expressed genes were obtained in the four cDNA libraries, including 2,840 upregulated genes and 5,276 downregulated genes (total 8,116 genes) in Figs. 8A and 8B. The downregulated genes in male and female flower buds of the two mating types were generally higher than those of the upregulated genes. In M1 vs. M2, 585 genes were upregulated and 772 genes were downregulated. In T1 vs. T2, 359 genes were upregulated and 1,463 genes were downregulated (Fig. 8B) (Data related to differential expressed genes have been submitted to NCBI database and TSA database, TSA SUB10193123).

In the signal transduction and synthesis pathways of plant hormones, we compared the differentially expressed genes in M1 vs. M2 and T1 vs. T2 libraries and found that the sex differentiation of flower buds was mainly related to ARF, SAUR, DELLA, TF, PR-1, etc (TSA SUB10193123, also uploaded as the Supplementary File, Data related to differential expressed genes. zip).

Differential expression of TFs in male and female flower buds

A total of 195 TFs were differentially expressed in M1 vs. M2, 168 TFs were differentially expressed in T1 vs. T2, and 83 TFs were differentially expressed in M1 vs. M2 and T1 vs. T2. The differentially expressed TFs included members of ARF, AP2, C2H2, ERF, bHLH, MADS, MYB, WRKY, NAC and other TF family members were included. MADS and bHLH families were dominant among the families. We have screened out many transcription factors related to plant hormones, and based on the results, we have identified transcription factors with differential expression in auxin and ethylene signaling pathways (Fig. 9). There were eight identical differentially expressed transcription factors in M1 vs. M2 and T1 vs. T2.

Verification of gene expression via real-time PCR

To verify the gene expression pattern identified by RNA-Seq data, the transcription levels of 12 randomly selected DEGs were detected by quantitative PCR. We used these 12 genes for qRT-PCR, with three biological replications in each group. The gene expression pattern detected by qRT-PCR was consistent with the RNA sequence data (Fig. 10). The qRT-PCR experiments confirmed that the DEGs obtained from the assembled transcriptome were accurate and the gene expression profile obtained from RNA-Seq data was reliable.

Determination of endogenous hormones in flower buds of J. mandshurica

To explore the relationship between the flower bud differentiation of J. mandshurica and plant endogenous hormones, we determined the content of eight plant hormones in flower buds at the early stage of physiological differentiation. The content of plant endogenous hormones affected flowering time. By measuring the plant hormones in the male and female flower buds of J. mandshurica at the early stage of physiological differentiation (Table 2), we found that the contents of auxin (IAA), brassinolide (BR), and salicylic acid (SA) demonstrated an overall upward trend. Gibberellin (GA) content first increased and then decreased, and a peak appeared on approximately April 10th. The content of cytokinin (CTK) fluctuated: up-down-up. The contents of ethylene (ETH) and abscisic acid (ABA) first decreased and then increased. Jasmonic acid (JA) content in flower buds of female precocious J. mandshurica was in a relatively stable state before April 10th and began to increase after this date, while JA content in flower buds of male precocious J. mandshurica revealed an overall upward trend. Interestingly, on April 2nd, the contents of IAA, GA, CTK, ETH, and SA were all found in the flower buds of J. mandshurica: T1 < T2 and M1 > M2. The female flowers of the female pre-type J. mandshurica open first, and the male flowers disperse pollen later. The male flowers of the male pre-type J. mandshurica are scattered first, and the female flowers open later. The contents of the above five hormones in the male and female flower buds of the two types of J. mandshurica were opposite, indicating that there were genes regulating gender differentiation in these five hormones, which was consistent with the differential expression of genes in the transcriptome data. We inferred that these hormones participate in gender differentiation to varying degrees and that there may be interactions. BR and JA contents were as follows: T1 < T2 and M1 < M2. The contents of BR and JA in female flower buds were lower than those in male flower buds, indicating that the regulation of these two hormones was more important in male flower buds. ABA content revealed that: T1 > T2 and M1 > M2, while the ABA content in female flower buds was higher than that in male flower buds, indicating that ABA regulation was more needed in female flower buds.