Genome-wide quantification of homeolog expression ratio revealed nonstochastic gene regulation in synthetic allopolyploid Arabidopsis.

5.5 MB·English

by  Akama, Satoru

Published online 13 January 2014 Nucleic Acids Research, 2014, Vol. 42, No. 6 e46
doi:10.1093/nar/gkt1376
Genome-wide quantification of homeolog
expression ratio revealed nonstochastic gene
regulation in synthetic allopolyploid Arabidopsis
Satoru Akama1,y, Rie Shimizu-Inatsugi2,y, Kentaro K. Shimizu2,* and Jun Sese1,*
1Department of Computer Science, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo
152-8550, Japan and 2Institute of Evolutionary Biology and Environmental Studies and Institute of Plant
Biology, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland
Received June 5, 2013; Revised November 26, 2013; Accepted December 14, 2013
ABSTRACT INTRODUCTION
Genome duplication is a common phenomenon in
Genome duplication with hybridization, or allopoly-
Eukaryota including animals, fungi and plants. It is
ploidization, occurs commonly in plants, and is
estimated that all flowering plant species have
considered to be a strong force for generating new
experienced one or more rounds of genome duplication
species. However, genome-wide quantification
(polyploidization) in their history (1). Moreover, 34.5%
of homeolog expression ratios was technically
of land plants are inferred to be polyploid relative to
hindered because of the high homology between their base chromosome number of the genus, suggesting
homeologous gene pairs. To quantify the homeolog their recent origins (2). The phenomenon is classified into
expression ratio using RNA-seq obtained from poly- autopolyploidization (derived from intraspecific genome
ploids, a new method named HomeoRoq was de- duplication)andallopolyploidization(orgenomeduplica-
veloped, in which the genomic origin of sequencing tion with hybridization). Allopolyploidization in particu-
readswasestimatedusingmismatchesbetweenthe lar is considered to have contributed to the adaptation to
broader and novel environmental niches (3–8). Polyploid
read and each parental genome. To verify this
plants often live in extreme habitats, such as cold arctic
method, we first assembled the two diploid parental
regions, or dry habitats such as in the Mediterranean
genomes of Arabidopsis halleri subsp. gemmifera
climate. Theoretical research performed during recent
and Arabidopsis lyrata subsp. petraea (Arabidopsis
decades has revealed that the duplication of some func-
petraea subsp. umbrosa), then generated a syn-
tionalgeneswouldbeadvantageous(9).Inaddition,many
thetic allotetraploid, mimicking the natural of high-yield agricultural species are polyploids, which
allopolyploidArabidopsiskamchatica.Thequantified suggests that genome duplication may have played a fun-
ratios corresponded well to those obtained damental role in the evolutionary history of speciation or
by Pyrosequencing. We found that the ratios of adaptation to new environments (10).
homeologs before and after cold stress treatment To identify the molecular mechanisms underlying po-
were highly correlated (r=0.870). This highlights the tential adaptation and diversification, the gene expression
presenceofnonstochasticpolyploidgeneregulation patterns of allopolyploids have been extensively studied
(11–13). Because the genome of an allopolyploid consists
despitepreviousresearchidentifyingstochasticvari-
of those of two (or more) closely related species, the two
ation in expression. Moreover, our new statistical
gene copies that originated from the parents are, in
test incorporating overdispersion identified 226
general, highly similar. These duplicated gene copies are
homeologs (1.11% of 20369 expressed homeologs)
called homeologs. Many studies of synthetic and natural
with significant ratio changes, many of which were
allopolyploid plants showed a nonadditive and unstable
related to stress responses. HomeoRoq would con- gene expression pattern compared with that of the
tributetothestudyofthegenesresponsibleforpoly- parental species, i.e. the sum of the gene expression
ploid-specific environmental responses. levels of homeologs is often different from the simple
*To whom correspondence should be addressed. Tel:[phone]; Fax:[phone]; Email: [email protected]
CorrespondencemayalsobeaddressedtoKentaroK.Shimizu.Tel:[phone];Fax:[phone];Email:[email protected]
yThese authors contributed equally to the paper as first authors.
(cid:2)TheAuthor(s)2014.PublishedbyOxfordUniversityPress.
ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommonsAttributionLicense(http://creativecommons.org/licenses/by/3.0/),which
permitsunrestrictedreuse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited.
e46 NucleicAcidsResearch,2014,Vol.42,No.6 PAGE2OF15
sum of those of the two parents. Synthetic and natural diploid has been studied in medical and evolutionary
allotetraploids of Arabidopsis suecica derived from research (24–33). These studies focused on model species
Arabidopsis thaliana and Arabidopsis arenosa has been with high-quality genomes. For example, in a study
used to study the expression pattern in polyploids of parent-of-origin expression in F1 mice, reads were
(13,14). In cotton species, the expression level dominance mapped on a reference that contained coding sequences
has long been discussed, with studies of these plants of both parental strains (22,27,31,33). However, methods
showing that the expression ratios of two different for allele-specific gene expression is not readily applicable
genomes are not even, and that gene expression could to polyploid studies because high-quality reference
show any possible pattern (11). Such novel expression genome of both of the parental species are usually not
patterns are proposed to contribute to the adaptation of available. In addition, methods in allele-specific gene ex-
polyploid species (5). pression are usually designed to focus on cis-trans regula-
However,thedifficultyindistinguishinghomeologshas tionorimprintingandthusratiochangesofhomeologsin
beenamajor limitationof expression studiesofpolyploid polyploid may not be analyzed.
species. Because of the high sequence similarity between When analyzing count data of RNA-seq reads, there
two homeologs, they may be detected similarly in micro- is a recognized problem of overdispersion where the
array or in polymerase chain reaction (PCR); thus, the variance is greater than expected from a simple Poisson
design of homeolog-specific probes or primers for each distribution (34–36). Overdispersion is attributed to un-
gene is typically necessary (15,16). Moreover, reference
controllable differences between replicates such as in
genome sequence data for polyploid species are not
microenvironment or experimental conditions. When
readily available in most cases. Geneticists have avoided
overdispersion is not taken into account, the number of
using polyploids as model species for two main reasons:
differentially regulated genes would be overestimated.
first,theirgenomesizeis,bydefinition,largerthanthatof
However, to our knowledge, little is known about
the parental diploid species and second, the genome
overdipersionintheexpressionratioofhomologsinpoly-
assembly of polyploid species is extremely difficult, if
ploid species.
not impossible (17) because genome-wide homeologous
Here we focus on the allopolyploid species Arabidopsis
sequences hinder correct assembly.
kamchatica(FischexDC.)K.Shimizu&Kudoh,whichis
Another obstacle to the separate detection of
derived from Far East populations of Arabidopsis lyrata
homeologs is the lack of large-scale observation tech-
and Arabidopsis halleri (37–40). A biogeographic study of
niques. Laborious molecular techniques such as real-time
Arabidopsis showed that the allotetraploid A. kamchatica
PCRorPyrosequencing (PyroMark)can oftenbe applied
(previously known as A. lyrata subsp. kamchatica) grows
to detect homeologs even with minor differences in single
in a broader climatic niche of temperature and precipita-
nucleotide polymorphisms (SNPs). However, these tech-
tion compared with diploid A. lyrata and A. halleri (41).
niquesallowtheidentificationofonlyalimitednumberof
Arabidopsis lyratatends tolive incolder regionsathigher
genes (18). Single-strand conformation polymorphism
latitudes, and A. halleri in warmer regions at lower lati-
assay was applied to study the expression ratio changes
tudes, while the niche of A. kamchatica ranges from cold
of 60 duplicated pairs in natural allopolyploid and 20
towarmhabitats.Thissuggeststhattheallopolyploidhas
pairs in synthetic allopolyploid in response to abiotic
adapted to a broad range of environments by combining
stresses (19). The use of standard microarray approaches
the two parental genomes. Moreover, a wealth of func-
excluding the cross-hybridization between two
tional, genomic and evolutionary data are available on
homeologous sequences is difficult when the number of
the model genus Arabidopsis including A. thaliana,
SNPs per probe is limited. Recent studies have applied
A. halleri and A. lyrata (42–45).
homeolog-specific microarray to cotton, but only a few
In the present study, we analyzed genome-wide
hundreds of homeologs were identified (20,21).
homeologous gene expression of allotetraploid plants
Recenttechnicaladvancesinnext-generationsequencing
using RNA-seq. We produced an artificial allotetraploid
(NGS) have allowed whole-genome transcriptomic
analysis, coined RNA-seq, for a variety of species of between A. lyrata and A. halleri. We developed a new
model and nonmodel organisms. Even though NGS has bioinformatics method, termed Homeolog Ratio and
been applied to the transcriptomic analysis of many quantification (HomeoRoq), to quantify the contribution
diploid and polyploid species, only a few studies have of gene expression of each homeolog. In this method, we
undertaken the challenge of distinguishing the expression assigned RNA-seq reads to each parental genome based
levels of homeologous genes. The most intensive analysis on the number of mismatches between the reads and
performedtodateusedcotton(22,23),inwhichESTdata the genome sequence of each parental species. Counting
and SNPs between partial genomic sequences were used. the number of the reads on each homeolog provided the
However, the experimental verification of the quantifica- expression level of the homeolog. The expression ratios
tion of the ratio of homeologs was not presented. between homeolog pairs, which were computed using
Furthermore, the number of analyzed homeologous this bioinformatics method, corresponded well to those
genes was less than half the total number of genes, measuredusingPyrosequencing.Wealsodevelopedastat-
which may have been due to the difficulty of computa- istical method to identify genes with the change of
tional analysis and limited parental genome information. homeolog expression ratio under different conditions.
Similar to homeolog-specific expression, allele-specific Although the ratio of homeologs before and after cold
or parent-of-origin specific gene expression analysis in stress treatment were highly correlated, the test identified
PAGE3OF15 NucleicAcidsResearch,2014,Vol.42,No.6 e46
226 homeologs with significant ratio changes, which these autopolyploidized F1 hybrid individuals were
includes many stress response genes. collected, germinated and grown for tissue sampling.
Plant incubation
MATERIALS AND METHODS All plants were grown in a climate chamber (22(cid:2)C, 60%
humidity, 150mE light intensity, 16h light/8h dark). For
Artificial allotetraploid plant production
cold treatment, plants were moved into another climate
The individual of A. halleri subsp. gemmifera chamber (4(cid:2)C, 80mE light intensity, continuous light)
(Matsumura) O’Kane & Al-Shehbaz used in this article and incubated for 7 days.
was collected from a population located in Tadaginzan,
Inagawa, Osaka, Japan. Arabidopsis halleri is self-incom- DNA extraction and Illumina sequencing
patibleandharborshighlevelsofheterozygosity,whichis
GenomicDNAfromA.halleriandA.lyratawasextracted
not ideal for genome assembly. We therefore conducted
fromleaf tissues using DNeasyplantkit (Qiagen) andthe
five rounds of self-fertilization by bud pollination in the
totalDNAwassequencedusingIlluminaHiSeq2000with
lab to reduce heterozygosity [clone of the individual
insert size of 200, 500 and 800bp. The library synthesis
halgem2 used in (39), and W302 with S-haplogroup A in
wasconductedattheFunctionalGenomicsCenterZurich
(46)]. Arabidopsis lyrata subsp. petraea (L.) O’Kane &
(http://www.fgcz.ch/)withTruSeqDNASamplePrepKit
Al-Shehbaz was collected by Valentin Yakubov from a
v2 (Illumina).
population at the banks of the Suharnaya river,
alluvium of Kolyma, Yakutia (Sakha Republic), Far RNA extraction and Illumina sequencing
East Russia [corresponds to lyrpet4 used in (39)]. The
Leaf tissues for RNA extraction were collected from
taxon has been treated as a species or as a subspecies,
the mature vegetative rosette of three individual
including Arabis umbrosa Turcz. ex Steud. (1840) (47),
allotetraploid plants. The RNA was extracted with
Arabis media N. Busch (1922) (48) and A. petraea subsp.
TRIzol (Invitrogen) and purified further with an RNeasy
umbrosaElven&D.E.Murray(2008)(49).O’Kane&Al-
Mini Kit (Qiagen). RNA concentration was measured
Shehbaz (1997) and Al-Shehbaz & O’Kane (2002) (50)
using Qubit (Invitrogen). The RNA samples (1mg each)
considered them as a synonym of A. lyrata subsp.
were submitted to the library synthesis by TruSeq RNA
petraea (48,51). Shimizu-Inatsugi et al. (39) did not find
Sample Prep Kit v2 (Illumina) at the Functional
heterozygosity in a few loci, which is consistent with high
Genomics Center Zurich and the library was sequenced
homozygosity due to auto-pollination. Following germin-
using Illumina HiSeq 2000.
ationfromasingleseed,asingleroundofself-fertilization
was conducted.
Genome assembly and gene annotation
Arabidopsis kamchatica (Fisch. ex DC.) K. Shimizu &
Kudoh is an allopolyploid that originated from A. halleri De novo assemblies were performed on A. lyrata and
subsp.gemmiferaandA.lyratasubsp.petraeainFarEast A. halleri independently. For assembly, we used
Asia(37,39).UsingthemutationrategivenbyKochetal. SOAPdenovo 1.05 (57) with the option ‘-d -D’. We
(2006)(52),thespeciationbetweenA.lyrataandA.halleri varied K-mer’s K size in SOAPdenovo, and determined
isestimatedtohaveoccurred2533980yearsago(95%CI: K=73forA.halleriandK=83forA.lyratabychecking
1307952–5166833 years ago) (53) or 337000 years ago N50length,thenumberofpredictedgenesandthecorres-
(95% CI: 272800–438200 years ago) (54) depending on pondence of the predicted genes to A. thaliana genes.
the analysis, while the origin of A. kamchatica by Assembly and raw read statistics of the assembled
allopolyploidization is estimated to have occurred 20417 genomes are shown in Supplementary Table S1. To
years ago (95% CI: 0–75460 years ago) (40). Thus, the predict gene regions from the genome sequences, we
allopolyploidization events occurred much more recently used AUGUSTUS 2.6 (58) with the A. thaliana model.
(one or two orders of magnitude more recently, in years) Because many short genes were predicted, we selected
than the split of the parental species. It should be noted genes whose predicted coding sequences were >200bp.
that the absolute value of the time estimation depends We detected orthologous genes between the predicted
largely on the mutation rate, using the mutation rate genes and A. thaliana genes in TAIR 10 (59) using a re-
given in (55), A. kamchatica is estimated to have ciprocal best hit (best-to-best) manner based on BLAST
originated 245070 years ago (95% CI: 37385–532953 E-values of 10(cid:3)15. The numbers of associated genes in
H- and L-genome were 21263 and 21166, respectively.
years ago) (40), although this will not affect the relative
We adopted the annotations of A. thaliana orthologous
time estimate.
genes for the gene annotation of the genes in A. halleri
Toobtainasynthetichybrid,wegeneratedanF1hybrid
and A. lyrata.
bycrossingA.halleri(female)andA.lyrata(male),which
isinthesamedirectionasinthefourindependentnatural
Mismatch ratio on coding regions among A. halleri and
hybridization events estimated from chloroplast (plastid)
two subspecies of A. lyrata
sequences (39). The obtained F1 seeds were germinated
and grown on sand until cotyledons developed fully, and To estimate the mismatch rate in coding sequences
the shoot apical meristem was treated with colchicine between our assemblies of East Asian A. halleri, that of
solution (0.5%) overnight (56). The seeds from one of East Asian A. lyrata subsp. petraea and the previously
e46 NucleicAcidsResearch,2014,Vol.42,No.6 PAGE4OF15
released A. lyrata subsp. lyrata sequenced by Joint from a homeolog, it indicates that HomeoRoq can
Genome Institute (JGI) (45), we performed BLAST quantify the homeolog-specific expression from the gene.
search between all gene pairs over the target species, and
took the average value of the identity scores calculated. Mapping of RNA-seq reads
The value between A. halleri and our A. lyrata subsp.
RNA-seq reads were obtained from synthetic
petraea is 2.44%, from which we expected that most of
A. kamchatica in two conditions (control and cold
Illumina paired-end reads could be mapped to both
stress) with three replicates. Supplementary Table S3
genomes and have an average of 4.93 mutations per
shows the statistics from the RNA-seq experiments. The
paired-end read (202bp length) used to identify genomic
reads were mapped onto the L- and H-genome independ-
originofthereads.Additionally,aconsiderablemismatch
ently.AlthoughthemappingsoftwareTopHat(60)isused
rate (0.85%) between our A. lyrata subsp. petraea from
frequently for RNA-seq, its upper limit of the number of
East Asia and JGI’s A. lyrata subsp. lyrata was found,
mismatches is small and allows the mapping of only
which suggests divergence at the subspecies or species
<50% of reads. We therefore used STAR-2.2.0 for
level (see ‘Discussion’ section).
mapping. We ran STAR without any special option and
it allowed the maximum 10 mismatches for the alignment
Homeolog identification
becausethealignmentaccuracywashighevenunderthese
Toidentifyhomeologs,weassociatedA.hallerigeneswith parameters (Supplementary Table S2). If a read was
A.lyratagenesusingBLASTwithanE-valuethresholdof mapped on multiple positions, the aligned position with
10(cid:3)15 and regarded the associations whose conserved a ‘primary’ flag (position with the best alignment score)
regions are >200bp as homeologs. We identified 31749 was adopted to quantify expression levels to avoid am-
homeologs.Resultsfromthismethodwerenotcompletely biguous result caused from multiple hits.
symmetrical:whenweselectedA.lyratagenesasthequery
and A. halleri as the database to perform BLAST, the Pyrosequencing
number of homeologs identified becomes 31560. The
Homeologous gene expression ratios were verified using
reason is that single H-homeologs may be associated
the Pyrosequencing of five genes. This technique has
with multiple L-homeologs caused from the division of
beenusedforallelefrequencydetection(18),DNAmethy-
single L-homeolog gene into two or more scaffolds,
lation analysis in combination with bisulfite conversion
and vice versa. As the scaffold N50 of the H-genome is
(61) and others during the past decade. Pyrosequencing
longerthanthatofL-genome,andthusH-genomeismore
can provide relative abundance of SNPs by comparing
reliable than L-genome, we selected A. halleri as a
the intensity of light emission generated by the incorpor-
standard in all analyses.
ation of allele-specific bases during the sequencing
reaction by polymerase. We applied this technique
Evaluation of HomeoRoq with simulated RNA-seq reads
to detect the frequency of homeologous cDNA that
To assess the measurement accuracy of our strategy, we originated from the H- and L-genomes.
generated simulated RNA-seq reads and estimated their The set of gene-specific (PCR amplification) primers
genomic origin. Even though each gene should exhibit and sequencing primer(s) was designed by PyroMark
variable expression in actual tissues, we generated the Assay Design v2.0 software (Qiagen). The gene-specific
simulated reads on the assumption that all genes had the amplification primers were designed at conserved regions
same expression level. We also supposed that the total between H- and L-homeologs, so that they include SNP
expression levels of H-homeologs were equal to those of positions inside the amplified region. The sequencing
L-homeologs. The average insert length of the simulated primer was designed at adjacent conserved region to a
reads was 200bp, and we introduced random mutation targeted SNP position(s) between H- and L-homeologs.
at 0.1% rate, mimicking the RNA-seq error rate by Two sequencing primers were designed for each set
Illumina. of amplification primers, and one or more SNPs were
We computed the accuracy of mapped reads by designated as the target per sequencing primer
checking the correctness of the gene where the reads are (Supplementary Table S4).
mapped. For the H-origin read, we checked whether the The RNA samples that were used for Illumina RNA
mapped position is in the gene from which the read was sequencing were converted into cDNA using the High
generated. For the L-origin read, we checked whether Capacity RNA-to-cDNA kit (Invitrogen). DNA frag-
the L-homeolog from which the read is generated is a ments were amplified with amplification primers by
homeolog of the H-homeolog on which the read is ExTaq polymerase (Takara). The amplified DNA frag-
mapped. The accuracies of H-origin and L-origin reads ments were sequenced using PyroMark ID system
were 99.9 and 89.4%, and their average was 94.6% (Qiagen) with the sequencing primers at the Genetic
(Supplementary Table S2). DiversityCenter,Zurich.Theobtainedpeakwasautomat-
To count homeologs whose expression ratios could be ically analyzed using the PyroMark system using AQ
quantifiedby HomeoRoq, we measured thepercentageof (allele quantification) mode to determine the SNP ratio.
simulated reads generated from H-homeolog identified as The multiple SNP ratios obtained by multiple sequencing
H-origin on each H-homeolog. We assume that when the primers for each gene fragment were averaged to decide
percentage of the identifiable reads is >50% of the reads the H- and L-homeolog ratios.
PAGE5OF15 NucleicAcidsResearch,2014,Vol.42,No.6 e46
Differential expression between H- and L-homeologs distribution was likely not to fit the distribution of ratios
especiallywhentheratiosbetweenthetwoconditionswere
To identify H- and L-homeologs whose expressions
not close, resulting in false-negative detection of ratio
changed statistically significantly under different condi-
changes.Wethereforeusedthetruncatednormaldistribu-
tions, the RNA-seq reads sorted either as H- or
tionwherethemeanandvarianceofthedistributionwere
L-homeolog using HomeoRoq were subjected to DESeq
computed from real data h and l .
(34). We set False Discovery Rate (FDR) threshold to ij ij.
Oncewefixedr,theprobabilityP(r)=P(R =R )was
0.05. C O
found from the distribution. Also P(h , l j r) was
Cj Cj
calculated from the normal distribution whose mean
Statistical method for the detection of ratio changes
is m and variance of f(m, r)+d. Similarly, P(h , l j r)
incorporating overdispersion i i i Oj Oj
was calculated, and hence P(h , l j R =R =r) P(h ,
Cj Cj C O Oj
Todetectthechangesintheratioofhomeologexpression l j R =R =r) was calculated. By selecting r from a
Oj C O
inresponsetothestresscondition,wedevelopedamethod truncated normal distribution whose mean and variance
tocheckwhethertheratioofH-homeologtoL-homeolog arecalculatedfromR,P(h ,l jr)P(h ,l jr)P(r)can
i Cj Cj Oj Oj
expression was changed between the conditions. In this be calculated. Because P(h , l ) is independent from the
ij ij
method, we calculated the probability that the homeolog sampling,calculationofP(h ,l jr)P(h ,l jr)P(r)one
Cj Cj Oj Oj
expression ratiosbetween thetwo conditions areidentical million times with random samplings and comparison of
from given RNA-seq data. Then, when the probability them with the value calculated from observed values
was less than significance level, the homeolog could be gave the value of P(R =R j h , l ) of a homeolog. We
C O ij ij
regarded as having different expression ratio under the computed the P-value for each homeolog. In total, we
conditions. observed 31749 raw P-values, and then we performed
LetR beahomeologexpressionratioinconditioni.Let Benjamini-Hochberg (BH) correction to them to
i
h and l be the numbers of H- and L-origin reads on compute corrected P-values for multiple tests.
ij ij
sample j in condition i. The probability that the ratios
between the two conditions C and O are equal can be Gene ontology analysis
described as P(R =R j h , l ). As we could not
C O ij ij Gene Ontology (GO) analyses were performed with the
compute this probability directly, we transformed it by
GOTermFinder Web site (62) on 27 February 2013. We
using Bayes theorem:
firstassociatedhomeologswithA.thalianaTAIRgeneID,
PðR ¼R jh ,l Þ¼Pðh ,l jR ¼R ÞPðR ¼R Þ=Pðh ,l Þ and then used GO terms related to the TAIR gene ID. In
C O ij ij ij ij C O C O ij ij
all, 20814 genes were associated with at least one GO
Supposethatthetwoconditionsareobservedindepend-
term. Among 226 homeologs with significant ratio
ently, it can be described as follows:
changes, 153 of these homeologs were associated with at
Pðh ,l jR ¼R ÞPðh ,l jR ¼R ÞPðR ¼R Þ=Pðh ,l Þ least one GO term.
Cj Cj C O Oj Oj C O C O ij ij
As we had only three samples from each condition, it
wasdifficulttodirectlycalculateP(h ,l jR =R )and RESULTS
Cj Cj C O
P(h ,l jR =R ).Toestimatetheprobability,weused
Oj Oj C O Generating synthetic allotetraploids
an expression-level model incorporated in the RNA-seq
read analysis (34,35). In the RNA-seq model, expression Inthepresentstudy,wegeneratedsyntheticallotetraploid
level obeyed a normal distribution, whose mean was plants from two parental taxa of A. kamchatica, i.e.
average expression level of each gene, and variance A. halleri subsp. gemmifera and A. lyrata subsp. petraea
depended on the mean expression level of each gene and (sensu O’Kane & Al-Shehbaz, 1997, which is also called
ashotnoise,whichwasuniquelydeterminedoneachcon- A. petraea subsp. umbrosa) (48,49,51), both collected in
ditionandeachgene.Inhomeologexpression,inaddition Far East Asia (Figure 1A). The genome size of each
tothetwovariables, thevariancedependedonexpression taxa was ca. 220Mb measured by flow cytometry
ratio of each homeolog pair. We therefore assumed that (CyFlow(cid:3) Space, Partec) with 40,6-diamidino-2-
the variance of homeolog expression levels was described phenylindole (DAPI) (CyStain UV precise P, Partec).
asf(m,R)+d,wherem wasthemeanexpressionlevelof The F1 hybrid individual of these two species was
i i i i
condition i, and d was the shot noise. f was computed by treated with the classic polyploid inducer colchicine, and
i
usingalocalregressionmethod(weusedthelocfitpackage allotetraploid seeds were obtained from the treated indi-
in R). Among the three parameters, if R was fixed, mean vidual. The three individuals that germinated from these
i
expression m was calculated from the total expression seeds were used for RNA-seq. The genome size of this
i
h +l multiplied by R. From the two values, regression synthetic allotetraploid was ca. 450Mb, which corres-
ij ij i
value f(m, R) was calculated. The shot noise ponds to the sum of those of the two parents.
i i
was calculated as the residue from the real variance to We denote the genomes of A. halleri and A. lyrata
f(m, R). as H-genome and L-genome, respectively, the origin
i i
Inrealcalculation,R wasunknown.Wethereforeused from A. halleri and from A. lyrata as H-origin and
i
a truncatednormal distribution ranged from0 to 1as the L-origin, respectively, and halleri-derived and lyrata-
priordistributiononP(R =R ).Althoughbetadistribu- derived homeolog as H-homeolog and L-homeolog,
C O
tion is frequently used as the distribution of ratios, this respectively.
e46 NucleicAcidsResearch,2014,Vol.42,No.6 PAGE6OF15
A B
Figure 1. Overview of the quantification of the expression level ratio of homeologs. (A) The generation of synthetic allopolyploid mimicking
A. kamchatica. HomeoRoq was applied to this species. (B) Homeolog discriminative RNA-Seq pipeline. Genome/gene sequences of parental
speciesareindividuallyassembled.RNA-seqreadsobservedfromthetargetallopolyploidaremappedontobothH-andL-genomesindependently.
Based on the number of mismatches, genomic origin of the reads is classified.
Method for quantifying the ratio of the expression A. lyrata subsp. petraea from Far East Asia. Moreover,
of homeologs we aimed to demonstrate that HomeoRoq works without
theassembledgenomereportedpreviously.Wesequenced
We developed a new algorithm HomeoRoq to detect
the genomic DNA of both parents of the synthetic poly-
changes in the ratios of the expression of homeologs
ploid using Illumina sequencing technology. The total
under different conditions using RNA-seq data obtained
sequence amounts were 36.5 Gb (166-fold coverage of
fromthesyntheticallopolyploid species.Ouranalysiswas
the estimated genome size) for A. halleri and 67.8 Gb
composed of two parts: the first part aimed to quantify
(308-fold coverage of the estimated genome size) for
homeolog-specific expression level (Figure 1B), and the
A. lyrata. The total genome lengths were 202.97 and
second part aimed to test statistically the ratio of change
221.14Mb, respectively, which were consistent with
in the expression of each homeologous gene pair under
the genome sizes estimated using flow cytometer.
different conditions.
The scaffold N50 sizes were 17686 and 7848bp,
During the first step of expression quantification,
respectively.
we assembled the genome of two parental taxa of
Next,wepredictedthegeneregionsofeachspecieswith
the allopolyploid A. kamchatica, i.e. A. halleri subsp.
gemmifera and A. lyrata subsp. petraea (Figure 1B, A. thaliana model parameters and identified 36737 genes
Table 1 and Supplementary Table S1). Because our for A. halleri and 35392 genes for A. lyrata (Table 1).
previous study (39) had shown that A. lyrata subsp. Homeologous gene pairs between them were identified
lyrata from North America, whose genome has been using a BLAST search (63) from all A. halleri genes to
sequenced and assembled by JGI (45), was not directly all A. lyrata genes. Based on the BLAST results with an
involved in the origin of A. kamchatica, we sequenced E-value threshold of 10(cid:3)15, a one-way hit strategy
PAGE7OF15 NucleicAcidsResearch,2014,Vol.42,No.6 e46
associated A. halleri genes with A. lyrata genes. This pro- allopolyploid accurately, we calculated the classification
cedure yielded 31749 homeologous gene pairs (Table 1). accuracies from simulated RNA-seq reads generated
Last, we identified the origin of RNA-seq reads from the L- and H-genomes. We generated 10 million
obtained from synthetic A. kamchatica. Here, we define artificial reads from the genomes, added the random
the genomic origin as whether the read was derived from sequence error and estimated the genomic origin of
the H- or L-genome. Each read was mapped separately reads using HomeoRoq. When the sequencing error rate
onto each of the H- and L-genomes and the number of was set to 0.1%, 72.2% of the reads were classified as
mismatches between the read and each genome was having either an H- or L-origin; in addition, 94.6% of
counted. The genomic origin was decided according to the classified reads were used to quantify correct
the lower number of mismatches (read 1 and read 2, in homeologs (see ‘Materials and Methods’ section). Even
Figure 1B). In the case of a perfect match of the read to if the error rate was varied, it had little effect on
both genomes, or the presence of an identical number of the accuracies of read classification (Supplementary
mismatches, the reads were considered ‘unclassified’ (read Table S2). Moreover, our method quantified homeolog-
3,inFigure1B)becausetheseconditionsdidnotallowthe specific expression from 27874 homeologs (75.9% of the
genomic origin of the read to be determined. If the read total genes, 87.8% of the homeologs) as far as they are
was mapped to only one of the parental genomes (read 4, expressed(see‘MaterialsandMethods’sectionfordetail).
inFigure 1B), we did notuse it,as using thistypeof read These simulations validated our algorithm as an effective
toquantifyhomeologexpressionlevelsmaycauseabiasin method for identifying the genomic origin of RNA-seq
theestimationofexpressionlevelsattributedtothediffer- reads.
ence in quality of the assembled genomes of two parental
species, even though part of them would truly represent RNA sequencing and homeolog-specific expression
species-specificfragments.Wecountedthenumbersofthe
RNA samples were extracted from leaf tissues of three
H- and L-origin reads on each gene on H-genome and
synthetic allotetraploid individuals grown at 22(cid:2)C
took the ratio of the reads on each gene, which is
(referred to as the control condition). These individuals
regarded as the expression level ratio of the homeologs.
were further incubated at 4(cid:2)C for 7 days, and leaf
To assess whether HomeoRoq can be used to estimate
tissues were sampled again (referred to as the stress con-
thegenomicoriginoftheRNA-seqreadsobservedforthe
dition) to extract RNA. In total, six RNA samples (three
individuals in two conditions) were sequenced using
Illumina (Table 2).
Table 1. Statistics of parental species’ genes
WeestimatedthegenomicoriginofeachRNA-seqread
Parental species A. halleri A. lyrata by mapping them onto H- and L-genomes as describe
above. In three control samples, 92.6 and 92.7% of the
Number of predicted genes 36737 35392 reads were mapped onto the H- and L-genome, respect-
Number of genes related 21263 (57.9%) 21166 (59.8%) ively, and 88.1% of them were mapped onto both
to A. thaliana
genomes, which was used to measure homeolog expres-
Number of homeologs 31749 (86.4% based on A. halleri)
Number of expressed homeologs 18928 (59.6%) 19186 (60.4%) sion. Of the reads that were mapped onto both genomes,
46.0, 43.0 and 11.0% were categorized as H-origin,
Between A. halleri and A. lyrata, similar number of genes was found L-originandunclassified,respectively(Table2).Thevisu-
and expressed. The relationships between predicted genes and genes in
alization of the mapping results confirmed that our
A. thaliana were computed using reciprocal best-hit strategy based on
BLAST alignment score. The expressed homeologs are counted from method was able to distinguish the genomic origin of the
the six RNA samples, and we regarded a gene/homeolog expressed RNA-seq reads (Figure 2). Similar results were obtained
when its RPKM is >0.1. from the samples with cold stress treatment (Table 2).
Table 2. Statistics of mapped results (total of three samples)
Control (three samples) Stress (three samples)
Sample condition Number of reads Rate in all reads Number of reads Rate in all reads
Number of all reads 119123927 102602463
Mapped on H-genome 110262849 92.60% 95775209 93.40%
Mapped on L-genome 110369840 92.70% 95026179 92.60%
Mapped on both genomes 104910672 88.10% 90702938 88.40%
Discarded reads (Read 4 in Figure 1) 14213255 11.90% 11899525 11.60%
Rate among reads Rate among reads
‘‘Mapped on both genomes’ ‘‘Mapped on both genomes’
H-origin reads (Read 1) 48278507 46.00% 40828935 45.00%
L-origin reads (Read 2) 45059363 43.00% 39073036 43.10%
Unclassified reads (Read 3) 11572802 11.00% 10800967 11.90%
e46 NucleicAcidsResearch,2014,Vol.42,No.6 PAGE8OF15
Figure 2. AvisualizationexampleofmappedreadsonH-genome.H-originreadsareperfectlymatchedtotheH-genome,whileL-originreadshave
a few mismatches to the H-genome. This confirms that the proposed method estimate genomic origin successfully.
The mapping results showed that 18928 (59.6%) consistent with that obtained using RNA-seq quantifica-
and 19186 (60.4%) genes were expressed [Reads Per tion. A known technical limitation of Pyrosequencing is
Kilobase per Million mapped reads (RPKM)>0.1] from that PCR amplification may be biased to a particular
H-homeologs and L-homeologs, respectively, and 17745 homeolog. Considering the potential technical bias
genes were expressed from both homeologs of the Pyrosequencing method (18), our estimation
(Supplementary Table S3). These statistics confirmed of homeolog ratios using HomeoRoq with Illumina
that both homeologs of the majority of genes were RNA-seq data works reasonably well.
expressed under these conditions.
Expression ratio change of homeologous pairs
Experimental validation using Pyrosequencing
Coldtreatmentisknowntoaffecttheexpressionlevelofa
To verify homeolog expression quantified using RNA-seq large number of plant genes (64,65). To find differentially
reads, we conducted a Pyrosequencing analysis using expressedH-andL-homeologs,thesequencereadssorted
PyroMark for five genes known to be involved in stress by HomeoRoq were subjected to standard RNA-seq
response[COLDREGULATED6.6(COR6.6),COR15B, analysis software DESeq (FDR<0.05). We found that
ETHYLENE RESPONSIVE FACTOR-2 (ERF2), the expression of 2897 and 4364 (9.12 and 13.7%) of
EARLY RESPONSE TO DEHYDRATION SIX-LIKE H- and L-homeologs, respectively, were changed signifi-
1 (ESL1) and AT1G16850] (Supplementary Figure S1). cantly by cold stress (Supplementary Figure S2A and B).
We amplified by PCR the target-gene fragments from the These results encouraged us to proceed to examine the
cDNAs that had been reverse-transcribed from the same changes in homeolog expression ratio.
RNA samples used for the Illumina RNA-seq and subse- Using the homeolog-specific RNA-seq data sorted by
quent HomeoRoq analysis. We measured the homeolog HomeoRoq at control and stress conditions, we checked
ratio at multiple SNP positions (see details in thedistributionofthehomeologexpressionratioforeach
Supplementary Table S4) and averaged the ratio of condition (Figure 4A and B). The overall histogram
H-homeolog (Figure 3). In most of the samples, the shapes of the homeolog expression ratio in both control-
results of Pyrosequencing were consistent with those of and stress-condition results were fairly symmetric, which
the HomeoRoq result. In particular, the RNA-seq confirmedthat bothhomeologswereexpressedinthetwo
results obtained for the COR15B gene were mostly conditions. The distributions created the impression that
inside the error bar of the Pyrosequencing results. Even the ratio may be determined randomly, as both distribu-
whentheRNA-seqquantificationsofafewindividualslay tions were similar to a normal distribution.
outside the error-bar range of the Pyrosequencing results, The ratios of homeologs, however, were highly
aswasthecaseforCOR6.6orESL1,theoveralltendency correlated between the two conditions (Figure 4C.
of homeolog expression ratio for whole individuals was Pearson’s r=0.870), suggesting that most homeologous
PAGE9OF15 NucleicAcidsResearch,2014,Vol.42,No.6 e46
COR6.6
COR15B
Control 1 Control 2 Control 3 Stress 1 Stress 2 Stress 3 Control 1 Control 2 Control 3 Stress 1 Stress 2 Stress 3
100 100
)
%
(
n 75 75
gi
ri
o
c
mi 50 50
o
n
e
G
25 25
0 0
HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM
ERF2 ESL1
100 100
)
%
(
n 75 75
gi
ri
o
c
mi 50 50
o
n
e
G
25 25
0 0
HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM
AT1G16850 (salt response)
100
HRQ=HomeoRoq
)
%
( PM=Pyroseqencing (PyroMark)
n 75
gi
ri
o
c H-origin ratio in HomeoRoq
mi 50
no H-origin ratio in PyroMark
e
G
25 L-origin ratio in HomeoRoq
L-origin ratio in PyroMark
0
HRQPM HRQPM HRQPM HRQPM HRQPM HRQPM
Figure 3. Comparison of homeolog expression ratios between the bioinformatics analysis and Pyrosequencing. The homeolog expression ratios of
five genes estimated by HomeoRoq and Pyrosequencing (PyroMark) are compared. The data are shown for three individuals at control (ctrl) and
stress conditions.
pairsmaintain theirexpression ratioevenaftercoldstress Fisher’sexacttest,thechi-squaredtestandthebinomial
treatment. Moreover, Figure 4C showed the existence of testaretypicallyusedtodetectexpressionratiodifferences
specifichomeologswithexpressionratiosthatweregreatly (22). However, by checking the variance of homeolog ex-
changed by cold treatment. pression levels, we confirmed that the distribution
e46 NucleicAcidsResearch,2014,Vol.42,No.6 PAGE10OF15
Figure 4. Theratioofhomeologexpressionlevel.(A)Histogramofexpressionratioincontrolcondition.(B)Histogramofexpressionratioinstress
condition. (C) Scatterplot of the ratio of homeolog expression level between two conditions. Each dot corresponds to one homeologous pair. Red
dots indicate homeologs whose expression ratio was significantly changed over the conditions. (D) Scatterplot of the ratio of homeolog expression
levels. Red dots indicate homeolog whose expression ratio was regarded as significantly changed over conditions with Fisher’s exact test.
was overdispersed when compared with a binomial red points). Supplementary Table S5 includes the list of
test (Supplementary Figure S3). In the presence of the homeologs.
overdispersion, highly expressed homeologs tend to be These 226 genes with significant change in homeolog
deemed significant even if the actual change is small. ratio were subjected to Gene Ontology analysis. Twenty-
TheuseofFisher’sexacttesttoidentifyhomeologexpres- three categories showed enrichment (adjusted P<0.01
sionchanges led tothe identification of11442 homeologs with Bonferroni correction), most of which include
(36.0%) that were deemed significant (adjusted P<0.05 similar sets of genes related to stress response (see
with BH multiple test correction). Figure 4D (red ‘Discussion’ section; Supplementary Table S6). In
points) indicates that homeologs with a ratio that was Supplementary Figure S4, the expression changes of 49
not greatly changed were regarded as significant in genes in the category of response to stress were depicted
Fisher’s exact test. This value appears too large with an arrow.
compared with the numbers of differentially expressed
H- and L-homeologs extracted by DESeq analysis (9.12
and13.7%), and thusFisher’s exact testdoes notseem to DISCUSSION
be an appropriate method.
A novel method to quantify the expression level ratio
To overcome this problem, we developed a statis-
of homeologs
tical test that can handle overdispersion. This method
computes the appearance rate of the observed counts Thelarge-scaleapproachesavailablecurrentlytomeasure
under the assumption that the two conditions have iden- homeolog-specific expression require deep knowledge
tical ratios (see ‘Materials and Methods’ section for about the genome/gene sequences of both parental
detail). This test led to the estimation of 226 homeologs species because they use known SNP positions to
(1.11% of 20369 expressed homeologs) as being signifi- measure homeolog expression. Yoo et al. (22) used
cant (adjusted P<0.05 with BH correction) (Figure 4C, RNA-seq and counted the number of reads on SNPs.