Understanding the molecular and systemic levels of hormone interactions will help us predict the effects of
disruption or over activation of certain parts of the network on the entire plant response. However, there is
largely unknown about the crosstalk of different phytohormones in allotetraploid oil crops Brassica napus.
In this
work, we treated 14-day-old B. napus seedlings by seven plant hormones at four time points, and used
RNA-seq for
examining the overlap in transcriptional effects among the various hormones. Differentially Expressed Genes (DEGs)
analysis was employed to reveal genes that co-regulated by seven hormones, and the regulatory networks among seven
hormones were constructed, which uncovered the new network between gibberellin and cytokinin. These studies reveal
that the long-term effects of all hormone treatments represent a “domino effect”, resetting many systems within
the plant, and provides a versatile resource of genome-wide gene expression for future hormone studies in plant
kingdom.
Plant materials and growth condition
ZS11 is a conventional cultivars of semi-winter Brassica napus. The gain-of-function mutant of BnaA6.RGA
(Bnaa6.rga-D)
and BnaRGA quadruple mutant (Bnarga) were generated in Brassica napus (Westar) using CRISPR/Cas9
technology
(Yang
et
al., 2017). The BnaA9.CXK2-OE and TCSn-GUS report line were obtained by transferred p1300-35S-BnaA9.CXK2-mCherry
and
TCSn-GUS vector into Brassica napus (Westar), respectively. All the materials were soaked in water for 7
days,
then
placed into black plastic box flooded Hoagland nutrient solution and grown for 7 days. The whole growth process
was
presented in a growth chamber with 16-h light/ 8-h dark photoperiod at 22 °C, with a light intensity of 100 μmol
m-2s-1.
Total RNA extraction and RNA-Seq library preparation
Total RNA was extracted by Trizol Reagent (Invitrogen Life Technologies, USA) according to the manufacturer's
instructions. A total amount of 1.5 μg RNA per sample was used to generate RNA-seq libraries by the TruSeq RNA
Sample
Preparation Kit (Illumina, San Diego, CA, USA). To select cDNA fragments of the preferred 200 bp in length, the
library
fragments were purified by AMPure XP system (Beckman Coulter, Beverly, CA, USA). DNA fragments with ligated
adaptor
molecules on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 15 cycle PCR reaction.
Products
were then purified (AMPure XP system) and quantified by the Agilent high sensitivity DNA assay on a Bioanalyzer
2100
system (Agilent). The library preparations were sequenced on a Hiseq platform (Illumina) by Shanghai Personal
Biotechnology Co., Ltd.
Identification and analysis of differentially expressed genes
The quality of the RNA sequencing reads was examined by FastQC (v0.11.9)
(https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Barcode adaptors and
low-quality reads (read quality
< 80
for paired-end reads) were removed by Trimmomatic (v0.38) (Bolger et al., 2014). Then, the filtered reads were
aligned
to the B. napus reference genome (ZS11) (Song et al., 2020) using Hisat2 (v2.1.0) (Kim et al., 2015) with
default
parameters. Bam files containing aligned reads were inputted into StringTie (v1.3.3b) (Pertea et al., 2016) to
measure
the expression level of genes. Gene-level raw count data files were generated using featureCounts (v1.6.4) (Liao
et
al., 2014). The raw count data were imported into Bioconductor package DESeq2 (Love et al., 2014) in the R
language to
identify the differentially expressed genes. Genes had a log2-converted fold change ≥ 1 or ≤ -1 with an FDR (False
Discovery Rate) ≤ 0.05 were considered as DEGs. The cluster analysis for DEGs under 7 hormone treatments was
performed
using Short Time-series Expression Miner (STEM) software with the maximum number of model profile set to 40 (Ernst
and
Bar-Joseph, 2006). The gene ontology (GO) of ZS11 was retrieved from tair (https://arabidopsis.org/index.jsp)
according to the homologous genes in BnTIR (Liu et al., 2021). Fisher's exact test was used to examine whether the
functional categories were overrepresented for the DEGs under 7 hormone treatments. Resulting P values were
adjusted
to Q values by the Benjamini-Hochberg correction and the FDR of 5% was applied.
Hormone crosstalk constructed
The genes involved in hormones biosynthetic, catabolism, and signal pathway were collected from previous reported
(Li et
al., 2017). Furthermore, the keywords of 'abscisic acid', 'ethylene', 'brassinosteroid', 'gibberellin', 'auxin',
'jasmonic acid', and 'cytokinin' were used to search hormones genes from GO file of ZS11. The highly
interconnected web
of genes involved in hormones predicts that treatment with one hormone should lead to changes in the levels of
multiple
hormones. Hormone crosstalk network were constructed based on the changes of hormones biosynthetic, catabolism,
and
signal genes. We defined hormone-responsive genes as those with the absolute of log2-converted fold change ≥ 1 in
at
least one time point or tissues. The upregulation effect under hormone treatment was defined as the number of
up-regulated hormone-responsive genes greater than 1.5 times that of down-regulated hormone-responsive genes. When
the
difference between up-regulated and down-regulated hormone-responsive genes was less than 1.5 times, the effect
was
ambiguous. Then, the statistical tests under hormone treatment were performed using a Fisher's exact test based on
the
following 2*2 table (Goda et al., 2008). We used the expression changes of ABA biosynthetic genes under IAA
treatment as
an example. Where the number of DEGs of 7 hormone biosynthetic gene under IAA treatment is designated
N1
(Up-regulated)
and N2 (Down-regulated); the number of DEGs of ABA biosynthetic gene under IAA treatment is
designated X1, and X2,
respectively. The significant effect of ABA biosynthesis genes under IAA treatment was generated based on the
following
2*2 table:
|
Up-regulated |
Down-regulated |
ABA Biosynthesis genes |
X1 |
X2 |
Not ABA Biosynthesis genes |
N1 - X1 |
N2 - X2 |
Co-expression analysis and gene network visualization
The hormone hub genes were defined as the differentially expressed genes under 7 hormone treatment,
simultaneously.
Hormone hub genes were assembled in a 'guide-gene set', which were used to build the co-expression network. All
hormone
genes were defined as a 'candidate-gene set'. The connectivity between guide genes and candidate genes were
retrieved.
The gene pairs were retained only if the absolute value of Pearson correlation coefficient (PCC) of gene pairs
were
larger than 0.5. The co-expression networks of hormone biosynthetic, catabolism, and signal hub genes were
constructed
based on the gene pairs. Cytoscape v3.6 software 'yFiles Organic Layout' was used to visualize the co-expression
network
when the PCC was larger than 0.5 (Shannon et al., 2003).
System architecture and software for database construction
As a module in the BnTIR, all the web interface of hormone transcriptome platform were build using HTML5, CSS3,
JQuery
(http://jquery.com), and Bootstrap (https://getbootstrap.com). The gene expression profile
heatmap, and the fold
change
of DEGs were visualized by Highcharts (https://www.highcharts.com/), ECharts (https://echarts.apache.org), and
plotly.js
(https://plotly.com/). The statistics information of genes was
managed by Handsontable (https://handsontable.com/)
and
DataTables (https://datatables.net). An 'Electronic
Fluorescent Pictograph' browser was used to exploring and
analyzing
the expression or the fold change of DEGs under hormone treatment (Winter et al., 2007). An interactive Scalable
Vector
Graphics (SVG) were designed for exploring genes involved in hormone pathway and crosstalk network.