BACKGROUND
apple fruit acidity refers to the sensory intensity of tartness of fruit flesh tissues.
the stronger the tartness taste, the higher the fruit acidity levels.
chemically, fruit acidity can be quantified by measuring fruit juice ph and/or titratable acidity .
it has been shown that several organic acids are present in mature apple fruit, including malic acid, quinic acid, citric acid and others, but malic acid consists of more than 90 % of the total , thereby largely determining fruit acidity.
for dessert apples, fruit acidity could be grouped into three categories by ta: low , normal , and high , and categories high and low are not acceptable  <cit> .
this makes fruit acidity an essential quality component not only for determining the fate of existing varieties, but also for prompting routine evaluations of acidity levels in apple breeding as genotypes of fruit acidity beyond the acceptable normal range can make up 25–50 % in breeding populations.
as a result, fruit acidity has long been an important subject area of investigations in apple genetics.
apple fruit acidity is primarily controlled by the major gene or qtl on chromosome  <dig>  called malic acid  alongside a significant qtl on chromosome  <dig> and a few other qtls with relatively smaller effects .
with regard to the allelic interactions at the major qtl ma, normal/high acid allele  is largely dominant over the low acid allele .
genotypes mama and mama have fruit of normal to high acidity while genotype mama of low acidity with little or none commercial value.
to focus on the most important genetic factor, we  <cit>  and others  <cit>  have recently isolated the major qtl ma.
the common findings were that the ma locus harbors two new members of the aluminum-activated malate transporter <dig>  gene family, called ma <dig> and ma <dig>  the studies further found that ma <dig> was expressed in significantly positive correlation with fruit acidity levels while the expression of ma <dig> was barely detectable in both high and low acid fruit, suggesting that it was gene ma <dig> rather than ma <dig> that was the very gene underlying ma  <cit> .
in addtion, a detailed analysis of the allele specific dna sequences of ma <dig> indicated that a single base mutation that would stop the protein translation process prematurely was almost completely associated with low acidity in a diverse set of apple germplasm studied, suggesting that the low acidity is caused by the malfunction of the ma <dig> protein due to the deduced truncation at the c-terminus  <cit> .
these latest findings have markedly increased our understanding on fruit acidity, but many questions remain to be answered.
for example, how does the landscape of transcriptomes differ between genotypes ma_  and mama at maturity stage?
if fruit acidity is governed by a gene network in which ma <dig> is the genetic determinant, what would the other possible members in the network?
what are the potential biological process and/or regulatory mechanisms responsible for the contrast acidity levels between genotype groups ma_ and mama?
the development of mrna sequencing  technology that unlocks the power of high throughput of next generation sequencing has provided an ideal means to address these questions.
since its inception  <cit> , rna-seq has been rapidly adapted in transcriptomics studies in plants such as arabidopsis  <cit> , grape  <cit> , maize  <cit>  and rice  <cit> .
in apple, rna-seq based studies have recently been reported as well .
more importantly, to begin resolving the low coverage issue of the current version of apple reference transcriptome, we have improved it with rna-seq reads from fruit of golden delicious , the source of the reference genome  <cit> , which is available at the genome database for rosaceae   <cit> .
we have also used the improved reference transcriptome through rna-seq approach to construct a co-expression gene network associated with developmental regulation of malate levels varying from  <dig>  to  <dig>  mg/g fw  in developing fruit of ‘golden delicious’ of genotype mama  <cit> .
the objectives of this study are to address the questions aforementioned.
to do so, we first sequenced  <dig> rna-seq libraries representing transcriptomes of mature fruit from ten apple varieties of genotypes mama, and mama and mama, and then mapped the rna-seq reads against the improved apple reference transcriptome  <cit> .
using ma <dig> as a guide gene, a series of downstream analyses was conducted, leading to identification of weighted co-expression gene network modules significantly correlated with malate, most significant acidity genes, intramodular hub genes, regulator genes, enriched gene ontology  terms and mapman sub-bins, and others.
to the best of our knowledge, this is first report attempting to understand the ma1-mediated gene network regulating fruit acidity in apples of diverse genetic background.
methods
plant materials and fruit acidity quantification
ten apple varieties of known genotype at the ma locus were chosen, including four of mama–britegold , sweet delicious , novosibirski sweet  and pi <dig> , four of mama–fuji , rome beauty law , cox’s orange pippin  and jonathan , and two of mama–empire  and granny smith .
genotypes mama and mama were jointly designated ma_ to represent a group that had at least one functional allele of ma <dig>  the trees were budded onto rootstock p <dig> and grown in a research orchard of cornell university at geneva, new york.
fruit of three replicates per variety and 8– <dig> fruit per replicate were harvested at maturity from  <dig> to  <dig> trees in fall  <dig> .
each fruit was cross-sectioned into two halves: one half were used for maturity evaluation and fruit juice extraction, and the other half were sliced and immediately frozen in liquid nitrogen for rna isolation and for quantitation of malate and other metabolites.
the evaluation of fruit maturity and fruit acidity was conducted as previously described  <cit> .
fruit of cornell starch index  <dig> – <dig>   <cit>  were considered matured and only matured fruit  were used for analysis.
two indicators of fruit acidity were measured: ph by a ph meter  and titratable acidity  by an autotitrator .
fruit organic acids were extracted and measured using an agilent 7890a gc/5795c ms  with the same configurations and settings as described previously  <cit> .
metabolites were identified by comparing fragmentation patterns with those in mass spectral libraries and quantified based on standard curves for each metabolite and the internal standard ribitol.
rna-seq library construction and data analysis
fruit rna isolation and rna-seq library construction were performed as previously described  <cit> .
briefly, total rna was isolated from 3 g of ground fruit flesh tissues and then treated with dnase i .
for mrna isolation and strand specific rna-seq library construction, nebnext poly mrna magnetic isolation module and nebnext ultra directional rna library prep kit for illumina  were used.
the mrna was isolated from 5 μg of total rna and was fragmented at 94 °c for 10 min.
the first strand cdna was reverse transcribed from the fragmented mrna with datp mix and second strand cdna was synthesized from the first strand cdna with dutp mix.
the resulting double strand cdna was end-repaired, adaptor-ligated and size-selected, and then followed by the user enzyme digestion of the second strand cdna.
the intact modified first strand cdna was pcr-enriched for 14– <dig> cycles to obtain the individual sequencing ready rna-seq libraries.
these libraries were multiplexed in equal amount for single end 100-base sequencing in three lanes of hiseq  <dig>  at the cornell university biotechnology resource center .
the ten apple varieties were sequenced in three biological replicates with one replicate per lane.
illumina sequencing of the pooled rna-seq libraries generated  <dig> fastq files of sequences with a total of  <dig>  million raw reads .
the raw reads were fed into bowtie  <cit>  to remove bar codes and adapters and then aligned to the rrna reference sequences allowing up to three mismatches, which were downloaded from the silva rrna database   <cit> .
the rrna depleted reads were imported into clc genomics workbench v <dig>   to trim low quality reads and/or bases using the quality limit of  <dig>  and the ambiguous limit of  <dig>  the resultant clean and high quality reads were mapped by clc genomics workbench  against the improved apple reference transcriptome  <cit> , which is available at url http://www.rosaceae.org/species/malus/malus_x_domestica/cu_rna_seq_genes of gdr  <cit> .
in addition to  <dig>  of the  <dig>  originally predicted genes or mdps , the reference transcriptome includes  <dig>  novel transcripts  <cit> .
for convenience, these novel transcripts will be referred to ‘genes’ and named ‘g######s’ as in ‘g101234’, and ‘mdp0000’  in the original gene ids will be abbreviated to ‘m’  hereafter.
gene expression levels were calculated and represented by reads per kilobase of exon model per million mapped reads   <cit> .
genes of rpkm > <dig>  were considered expressed according to a previous study  <cit> .
identification of fruit acidity associated genes
fruit acidity associated genes  were defined for those that were differentially expressed between genotype groups ma_ and mama and/or those that were expressed similarly as the guide gene ma <dig>  to identify differentially expressed genes  between the two genotype groups, the gene expression data rpkms were subjected to baggerly’s test  <cit> .
to control the false discovery rate , the original p values in baggerly’s test were adjusted for multiple testing using benjamini-hochberg correction  <cit> .
the cutoff for declaring a gene expressed differentially was pfdr < <dig> .
however, any genes with an rpkm range < <dig> in the  <dig> rna seq samples were removed from the list.
to identify genes expressed similarly as ma <dig>  two filters were applied: 1) the mean difference in rpkm between the two genotype groups was greater than ten , and 2) the fold change in mean rpkm between the two genotype groups was greater than  <dig>  .
genes passed both filters were considered of similar expression as ma <dig>
inferring fruit acidity associated co-expression gene network modules
the co-expression gene network modules  were inferred from the fruit acidity associated genes  using weighted gene co-expression network analysis , an r package  <cit> .
before network inferring, the expression data  were normalized by square root transformation.
the automatic one-step network construction and module detection method with default settings were used, which include an unsigned type of topological overlap matrix , a power β of  <dig>  a minimal module size of  <dig>  and a branch merge cut height of  <dig> .
the module eigengene  value was calculated and used to test the association of modules with acidity in the ten genotypes of  <dig> samples.
gene significance , total network connectivity , and module membership  that is also known as eigengene-based intramodular connectivity , were calculated for each of the  <dig> faags in the ten modules.
the most significant module  was visualized using cytoscape  <dig>   <cit>  and were also analyzed using network analyzer  <cit> , a cytoscape plugin.
identification of regulator genes
detection of regulator genes was performed using lemon-tree software  <cit> , which is an improved and ‘one-stop shop’ version of the stochastic bayesian module network algorithm learning module networks , which requires run scripts on both matlab and perl  <cit> .
using expression data from the  <dig> faags, ten independent gibbs sampler runs  <cit>  were conducted with lemon-tree.
from these ten runs, we were able to generate a single set of tight independent clusters, and then select  <dig>  of at least ten genes per cluster for testing  <dig> candidate regulators.
the candidate regulators were genes annotated as transcription factors or signal transducers in go terms and/or mapman  bins in the  <dig> genes.
we considered candidate regulators as regulators when their probabilistic  scores associated with their assignment to any of the  <dig> tight clusters were within top 1 %.
functional annotation and enrichment analysis of wgcna modules
we previously used both mapman bins  <cit>  and gene ontology  terms to annotate the improved apple reference transcriptome  <cit> .
briefly, for annotation with mapman  bins, several databases were blast-searched with the apple reference transcriptome of  <dig>  genes using mercator, a web-based search tool  <cit> .
the searched databases include: tair-arabidopsis tair proteins , ppap-swissprot/uniprot plant proteins, chlamy-jgi chlamy release  <dig> augustus models, oryza-tigr <dig> rice proteins, kog-clusters of orthologous eucaryotic genes database, cdd- conserved domain database, and ipr-interpro scan.
for go term annotations, relevant information for apple genes  is available at gdr  <cit>  and was directly adapted.
however, go terms for the new transcripts  <cit>  were obtained by blast2go program  <cit>  together with blastx search against the ncbi reference proteins database .
enrichment analyses of the wgcna modules were all conducted by hypergeometric tests using clc gw against the same background of  <dig>  expressed genes.
the hypergeometric test is similar to the unconditional gostats test  <cit> , and the cutoff for significant enrichment is pfdr< <dig> .
depending upon actual test, the background was defined by  <dig>  annotated genes in mapman  bins, by  <dig>  in go biological process , or by  <dig>  in go molecular functions .
qrt-pcr analysis
qrt-pcr analyses were conducted with the same total rna samples  used for the rna-seq library construction.
two micrograms of total rna were reverse transcribed using the superscript iii rt module .
the resulting first strand cdna was diluted by 5-fold and used as templates for qrt-pcr analysis on a lightcycler  <dig> , where an apple actin gene  was served as reference.
the primer sequences for the reference gene and the eight target genes were listed .
the actual qrt-pcr reactions and subsequent expression analysis were conducted similarly as described previously  <cit> .
statistical analyses
statistical analyses, such as anova, student’s t test and tukey’s hsd  test, were performed using jmp pro <dig> .
RESULTS
fruit metabolite profiling and acidity evaluation
fruit metabolite profiling was conducted using gc-ms with three biological replicates in the ten apple varieties.
a total of  <dig> metabolites were quantified, including  <dig> soluble sugars and seven organic acids.
among the  <dig> sugars, only sorbitol showed a significantly  higher concentration in genotype group mama  than in ma_ .
other sugars, especially the highly abundant fructose, sucrose and glucose, had similar levels across these varieties.
among the seven organic acids , malate was the most abundant  while quinate  and maleate  were distantly followed as the second and third abundant, respectively.
malate concentrations varied significantly  in the ten apple varieties , and were more than 3-fold higher in genotype group ma_  than in mama  .
the concentrations of maleate, succinate, fumarate and citrate were also significantly higher in genotype group ma_ than in mama , but they were at much lower levels  .
however, quinate showed no significant difference between groups ma_ and mama.
the trend of fruit juice titratable acidity  was similar to that of malate in the ten varieties  while an expected reverse trend was observed for fruit juice ph .fig.
<dig> evaluation of fruit acidity and causal organic acids.
standard deviations were shown with the error bars.
a ta and malate concentrations measured in the ten apple varieties, which were abbreviated as the following: b britegold, n novosibirski sweet, p pi <dig>  s sweet delicious, c cox’s orange pippin, e empire, f fuji, g granny smith, j jonathan, r rome beauty law.
the first four are of genotype mama whereas the last six of genotype ma_.
b mean ta and malate concentrations in genotype groups ma_ and mama.
c concentrations of succinate, fumarate, citrate and maleate in genotype groups ma_ and mama.
d fruit ph readings in the ten apples
transcriptome analysis in genotype groups ma_ and mama
rna-seq data analysis was performed using the improved apple reference transcriptome  <cit> .
after removing reads derived from rrna and/or of low quality, the reads for rna-seq mapping were  <dig>  million in total .
overall the total mapped and uniquely mapped reads were  <dig>  million  and  <dig>  million , respectively.
the mean mapped reads per sample were  <dig>  ±  <dig>  million in total and  <dig>  ±  <dig>  million in unique.
however, the sample for replicate i of granny smith was excluded from the downstream analysis as  <dig>  % of its raw reads was removed by the rrna and quality filters .
the maximum number of expressed genes was  <dig>  in this study, a tally counting all genes of rpkm > <dig>  in at least one of the  <dig> samples.
however, measuring by the threshold in group means , the expressed genes were  <dig>  .
of these,  <dig>  were expressed in both genotype groups, and  <dig> and  <dig> were specifically in ma_ and mama, respectively .fig.
<dig> summary of rna-seq data analysis.
a venn diagram representation of the number of genes expressed  in genotype groups ma_  and mama .
the cutoff is group mean rpkm > <dig> .
b mean number of total and uniquely mapped reads per sample in genotype groups ma_ and mama.
standard deviations were shown with the error bars
the total mapped reads of the expressed genes were  <dig>  ±  <dig>  million per sample in the ma_ group and  <dig>  ±  <dig>  million per sample in mama .
the uniquely mapped reads per sample in the ma_ and mama groups were  <dig>  ±  <dig>  million and  <dig>  ±  <dig>  million, respectively .
identification of malate associated co-expression gene network modules using wgcna
transcriptome-wide comparison of the  <dig>  expressed genes revealed  <dig> fruit acidity associated genes , including  <dig>  of differentially  expressed genes  between groups ma_ and mama, and  <dig> expressed similarly as the guide gene ma <dig> .
using wgcna, a co-expression gene network was constructed from the  <dig> faags .
within the network, the total network connectivity of genes showed a significant correlation  with the gene significance for malate , suggesting that more connected genes would have higher gene significance for acidity in the  <dig> faags.fig.
<dig> weighted gene co-expression network analysis  of fruit acidity associated genes  in the ma_ and mama genotype groups.
a hierarchical cluster tree showing ten modules of co-expressed genes.
each of the  <dig> faags is represented by a leaf in the tree while each of the ten modules by a major tree branch.
the lower panel shows modules in designated colors.
b correlation between gene total network connectivity and absolute gene significance for malate in the whole co-expression gene network.
c module-fruit acidity correlations and corresponding p-values .
the left panel shows ten modules  and the number of their member genes.
the right panel is a color scale for module trait correlation from − <dig> to 1
further wgcna analysis identified ten network modules in the co-expression network, designated turquoise, black, brown, pink, magenta, red, green, blue, yellow, and grey  .
investigating the relationships between module eigengene and acidity uncovered that the correlation coefficients varied widely from − <dig>  to  <dig>  in malate and from − <dig>  to  <dig>  in ta .
the eigengenes  of five modules turquoise, black, brown, blue and yellow showed significant correlations  with malate, suggesting these five modules had greater relevance in fruit acidity although blue and yellow were negatively correlated .
the five modules comprised  <dig> genes with module turquoise being the largest of  <dig> genes.
notably, the guide gene ma <dig> was assigned to turquoise, which also had the highest modular correlation  with acidity .
inspecting the correlation between the module memberships  and gene significance  revealed that all modules were significant except blue .
again, module turquoise had the most significant correlation between mm and gs , indicating genes of higher mm values were more likely of greater gs in turquoise than in others.fig.
<dig> analysis of module turquoise.
a module eigengene values across the  <dig> samples, including  <dig> in ma_ on left and  <dig> in mama on right.
samples are represented by the combination of a letter  and a number .
b correlation between module membership  and gene significance  for malate
identification of the most significant acidity genes
the  <dig> faags had a range of absolute gene significance  for malate from  <dig>  to  <dig>  .
we define that the most significant acidity genes  are those of gs for malate ≥ <dig> , the observed gs of the guide gene ma <dig>  this allowed to identify  <dig> msags , including three transporters , two receptor kinases , one ethylene-insensitive <dig> -like transcription factor , one pyruvatedehydrogenase complex component e <dig> , one endo- <dig> -beta-xylanase/ hydrolase , one glutamine synthetase cytosolic isozyme , one photosystem ii subunit r , and seven  of unknown function.
the  <dig> msags were distributed in three modules with  <dig> in turquoise, three in black and another three in blue , further indicating module turquoise is of greater relevance in acidity.fig.
<dig> the most significant acidity genes .
a expression profile  of the  <dig> msags.
each row represents a gene as listed on right.
each column stands for a sample as specified by sample names  at lower panel.
clustering of genes and samples was shown with distance on left and top panels, respectively.
the expression of genes is color coded from low  through mean  and through high .
the bright yellow color across the samples for ma <dig> and m <dig> indicate the expression of these two genes were all close to or exceeded the high end of the scale at bottom.
b distribution of the  <dig> msags in modules.
modules with zero genes are not shown.
note that genes of ids beginning with letter ‘g’ are referred to the novel transcripts  <cit> , and the apple reference gene ids are abbreviated
identification of intramodular hub genes and evaluation of their total network connectivity changes in genotype groups ma_ and mama
to identify intramodular hub genes, we similarly used the guide gene ma <dig>  which had an mm of  <dig> , ranking 11th in module turquoise.
based on this information, we defined top  <dig> genes in mm values from each of the five modules as intramodular hub genes, leading to identification of  <dig> hub genes .
to evaluate how the  <dig> hub genes behaved in genotype groups ma_ and in mama, their respective total network connectivity  were calculated and compared with the overall group of both genotypes .
most of the  <dig> hub genes showed a large increase in ktotal in mama while a marked reduction in ma_, suggesting a major pattern of negative correlation between the hub genes’ ktotal changes and acidity levels .
this pattern was true in modules blue and yellow , largely true in brown and black , but largely untrue in turquoise , seemingly providing another means depicting the module-acidity relationship shown early .fig.
<dig> intramodular hub genes and their total network connectivity changes in genotype groups ma_ and mama when compared with the combined groups of both genotypes.
shown are modules black, blue, brown, turquoise and yellow
there were three additional patterns in ktotal changes : the first was observed only on m <dig> in turquoise, i.e.
a reverse trend by having an increased ktotal in ma_ and a decreased ktotal in mama.
m <dig> encodes a leucine-rich repeat receptor kinase .
the second was represented by ma <dig>  i.e.
the ktotal was decreased in both ma_ and mama, suggesting reduced ktotal of ma <dig> in mama might have important ramifications for low acidity.
this pattern was also observed on m <dig>  binding protein involved in calcium signaling), m <dig>  and m <dig>  in turquoise, m <dig>  m <dig> and g <dig> in black, and m <dig> and m <dig> in brown .
the third pattern was characterized by an increased ktotal in both ma_ and mama, which was observed on m <dig> , m <dig>  and g <dig>  in turquoise, and g <dig> and m <dig> in black .
to graphically view the intramodular connectivity of these hub genes, the  <dig> hub genes along with other co-expressed  genes in module turquoise was presented using cytoscape .
it showed that the  <dig> hub genes were of the most edges from  <dig> to  <dig>  and were in the module network core .fig.
<dig> graphic view of module turquoise.
a cytoscape representation of co-expressed genes with edge weight ≥ <dig> .
the edge number of the genes ranges from  <dig>  to  <dig> , which is color coded by the scale at top left from blue through red.
some important genes are noted, including  <dig> intramodular hub genes m <dig>  m <dig>  m <dig>  m <dig>  m <dig>  m <dig>  m <dig>  m <dig>  m <dig>  g <dig>  and ma <dig> , two regulator genes m <dig> and m <dig>  and one hub and regulator gene m <dig> .
b distributions of edges.
the  <dig> hub genes are among those of most edges  of weight ≥ <dig>
identification of regulator genes
regulator genes are not readily identifiable from the wgcna co-expression network modules as they are non-directed.
to identify regulator genes, we used lemon-tree software  <cit>  first to generated  <dig> tight clusters  covering  <dig> of the  <dig> faags through ten gibbs sampler runs  <cit> , and then to examine the  <dig> candidate regulators  with the clusters.
to statistically validate the results obtained from the  <dig> candidate regulators, another  <dig> randomly selected genes were used as control  and examined in parallel with the  <dig> candidate regulators.
as expected, the probabilistic  scores for genes assigned as a regulator for the  <dig> tight clusters showed significant difference  between the  <dig> candidate regulators  and ck .
at the top 1 % of p. scores ,  <dig> of the  <dig> candidate regulators were identified as regulators and assigned to  <dig> of  <dig> tight clusters .
of the  <dig> regulators, five  were from module turquoise , two  from brown , and five from modules green, red and pink which were irrelevant for acidity .
below are brief descriptions of three regulator genes, which were among the  <dig> msags or among the  <dig> most connected intramodular hub genes in module turquoise.fig.
<dig> regulator genes and their assigned tight clusters.
a regulator m <dig>  and cluster  <dig>  of  <dig> genes.
the expression of genes is color coded from low  through high .
each row stands for a gene  and each column for a sample named at the bottom .
the regulator is assigned based on the hierarchical tree on top, which indicates how samples were clustered together with the red vertical lines.
the short red vertical line in the sample list shows where the two primary clades diverge while the short green vertical line marks where the secondary clades depart within a primary clade.
the genotype names in red indicate they are not in agreement with one of the two primary clusters ma_ or mama.
b regulator m <dig> and cluster  <dig> of  <dig> genes.
c regulator m <dig> and cluster  <dig> of  <dig> genes.
elements in  and  and their contents, formats and messages are the same as those in
the first is m <dig>  encoding an ein3-like transcription factor.
m <dig> was assigned to six tight clusters  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> with p. scores  <dig> – <dig>  .
notably, lemon-tree clustered ma <dig> into cluster  <dig>  to which m <dig> was assigned as a regulator with the highest p. score  , suggesting m <dig> might potential regulate ma <dig>  m <dig> was one of the  <dig> msags  and of high mm value  in module turquoise.
in the six m <dig> regulated clusters, the fruit samples were all clustered into two primary clades, which were nearly perfectly alongside the line dividing the two genotypes ma_ and mama .
such clustering pattern of fruit samples were not observed in any of the remainder  <dig> clusters regulated by the other  <dig> regulator genes , further supporting the potential regulatory role of m <dig> in fruit acidity.
the second is m <dig>  encoding a receptor-like protein kinase.
m <dig> was assigned to cluster  <dig> with p. score  <dig>  .
similar to m <dig>  m <dig> is also one of the  <dig> msags  that had a high mm value  in module turquoise.
the third is m <dig>  one of the  <dig> intramodular hub genes in module turquoise .
m <dig> was assigned to cluster  <dig> with p. score  <dig>  .
as aforementioned, m <dig> encodes an atiqd2-like cam binding protein involved in calcium signaling and behaved similarly as ma <dig> in ktotal changes in genotype groups ma_ and mama.
detection of enriched go terms and mapman  bins in network modules
to detect enriched go terms and mapman  bins in the five significant modules, a series of hypergeometric tests were performed against the background of  <dig>  expressed genes.
overall, eight go terms and ten mapman sub-bins were enriched  in modules turquoise, brown, blue and yellow, but not in black, and no terms were repressed .
due to multifaceted roles of genes and due to both go and mapman systems were used, some of the enriched terms had the same or mostly overlapping member genes.
in total, the  <dig> enriched terms covered  <dig> non-redundant genes, including  <dig> in module turquoise,  <dig> in brown,  <dig> in blue and  <dig> in yellow .table  <dig> go terms and/or mapman  bines significantly enriched in wgcna modules significantly correlated with fruit aciditya
athe go terms and/or mapman  bins that have the same or mostly overlapping member genes in the same module are highlighted with same color
bbai et al.
<cit>
in turquoise, two go  terms  and two mapman sub-bins  were enriched, which are all related to photosynthesis .
however, go: <dig> and sub-bin  <dig> .
<dig>  had the same five member genes.
notably, m <dig>  a member in go: <dig>  is one of the  <dig> msags .
in brown, three go  terms and one mapman sub-bin  <dig> .
<dig>  were enriched .
the three go terms were highly related as the four member genes in go: <dig>  and go: <dig>  were identical and were also among the six member genes in go: <dig> .
m <dig>  a common member gene in these go terms, is one of the  <dig> hub genes in brown .
in blue, we detected three go terms and five mapman sub-bins, the most in this study .
but the ten member genes in go: <dig>  were fully covered by the  <dig> genes in sub-bin  <dig> .
<dig> .
similarly, the two genes in sub-bin  <dig> .
<dig>   and three genes in go: <dig>  were covered by the  <dig> genes in go:  <dig> .
sub-bins  <dig>  ,  <dig> .
<dig>  and  <dig>   were unique.
remarkably, four of the  <dig> hub genes in blue were present, i.e.
m <dig> and m <dig> in go:  <dig>  and m <dig> and m <dig> in sub-bin  <dig>  .
in yellow, two mapman sub-bins  <dig>   and  <dig> .
<dig>   were enriched, and m <dig> from sub-bin  <dig>  is one of the  <dig> hub genes in yellow .
qrt-pcr confirmation of gene expression
to evaluate if and how the rpkm values reflected the gene expression levels, a set of eight genes were analyzed using qrt-pcr .
the eight genes include ma <dig>  the ein3-like regulator m <dig>  another most significant gene for acidity m <dig> encoding a protein serine/threonine kinase, and five others.
the data not only confirmed that the relative expressions of the eight genes in qrt-pcr were significantly  correlated with their rpkm values in rna-seq, but also confirmed their differential expression between the two genotypes groups ma_ and mama.fig.
<dig> expression confirmation of eight selected genes using qrt-pcr.
a-h the normalized expression of target genes relative to a control gene  in qrt-pcr was shown in light grey, and their corresponding rpkm values from rna-seq were in black.
the correlation coefficient  and associated p value  were shown accordingly.
standard deviations were shown with the error bars
discussion
identification of individual fruit acidity associated genes
in this study, ten apple varieties were carefully chosen to represent the two genotype groups ma_  and mama of contrast fruit acidity levels .
these varieties were grafted on the same rootstock p <dig> and grown in the same orchard with the same management routines to curtail the environmental impacts on gene expression.
however, since the apples are of diverse genetic background and matured at varying dates, fruit transcriptome at maturity were inevitably influenced by these genetic and environmental variables.
to minimize such biases, three replicates per variety and at least four varieties per genotype group were sampled and used for rna-seq analysis.
this strategy appeared to be appropriate and adequate for the defined objectives in this study as we identified  <dig> fruit acidity associated genes .
in developing fruit of ‘golden delicious ’ of varying malate levels,  <dig> differentially expressed genes  were identified  <cit> .
in comparison of the two lists of genes,  <dig>  of the  <dig> faags were not present in the  <dig> degs, suggesting most genes relevant for fruit acidity variation across diverse genetic background were largely different from those important for malate variations in the developing fruit of gd.
overall, we found  <dig>  genes expressed  in the ma_ group and  <dig>  genes in the mama group.
as there were  <dig>  genes expressed in both groups, the total number of expressed genes was  <dig>  .
varying rpkm thresholds were employed in published studies  <cit> , we view rpkm > <dig>  used previously  <cit>  fit this set of rna-seq data well.
of the  <dig>  expressed genes,  <dig> were regarded faags.
given the significant correlation between rpkm in rna-seq and relative expression in qrt-pcr in the eight genes tested, we consider that the majority, if not all, of the genes relevant for fruit acidity variation at maturity in this study were included in the  <dig> genes.
co-expression gene network modules and go terms and mapman  bins relevant for fruit acidity
wgcna  <cit>  is one of the programs used widely for inferring co-expression network modules .
using wgcna, we identified five modules highly associated with fruit acidity .
hypergeometric tests identified eight go terms and ten mapman  bins that were enriched in four of the five modules .
in the most significant module turquoise that includes ma <dig>  the two enriched go terms and two mapman sub-bins are all related to photosynthesis .
in the developing fruit, the two go terms and the two mapman sub-bins were also enriched  <cit> , and  <dig>  of the  <dig> member genes in the four enriched terms were identified in both studies .
this suggested that fruit photosynthesis could likely be a part of the mechanism for the contrast malate levels between the two genotype groups at maturity.
it was proposed that the relative more active photosynthesis in young fruit may facilitate malate biosynthesis, thereby contributing to higher malate levels in young fruit than in maturing fruit  <cit> .
however, chloroplasts are only present in the hypodermal green tissues and the inner perivascular tissues in mature apple fruit, and are essentially absent in parenchyma cells that constitute the major portion of flesh  <cit> .
how enhanced photosynthesis activities in a fraction of cells would lead to increased overall malate levels in fruit remains to be understood.
in module brown, the enrichment of mapman sub-bin  <dig> .
<dig>  was unique to this study.
however, the three enriched go terms  were essentially a reflection of mapman sub-bin  <dig>   , which was enriched in negative malate dependent manner in developing fruit  <cit> .
in module blue, the enrichment of mapman sub-bins  <dig>  ,  <dig> .
<dig>  and  <dig>   were identified in this study only.
but go: <dig> , and mapman sub-bins  <dig> .
<dig>  and  <dig> .
<dig>   were also enriched in developing fruit  <cit> .
in module yellow, the enrichment of two mapman sub-bins  <dig>   and  <dig> .
<dig>   was unique in these diverse mature fruit.
overall,  <dig> enriched terms , including four from turquoise , three from brown , and five from blue , were also enriched in developing fruit of ‘golden delicious’  <cit> , whereas the remainder six terms  were enriched specifically in this study .
this suggested that common and unique processes relevant for malate variations in both the diverse mature fruit and the developing fruit of ‘golden delicious’ exist.
however,  <dig> of the  <dig> genes in blue and all of the six genes in brown were not found to be relevant in developing fruit although  <dig> of the  <dig> genes related to photosynthesis in turquoise were , suggesting the common processes may operate but most likely with their unique gene sets.
transcriptional regulation of fruit acidity mediated by ma1
in an ongoing survey of the usda malus repository at geneva, new york, we identified  <dig> diverse apple accessions of genotype ma1ma <dig> using marker caps <dig> that can detect the stop codon leading snp in ma <dig>  <cit>  and found that each of them had fruit ph > <dig>  , a typical low fruit acidity characteristic in genotype group mama.
these data, together with the previous studies  <cit> , suggested us that ma <dig>  a member of the almt <dig> gene family  <cit> , is the genetic determinant underlying the ma locus controlling acidity levels in mature apple fruit although the functional genetic complementation data are not available.
ma <dig> is thought to exert its effect on fruit acidity by both the truncation leading mutation and expression levels, where the mutation plays greater role  <cit> .
in this study, ma <dig> showed not only high gene significance for malate , but also high intramodular connectivity in module turquoise , further supporting this notion.
however, in our previous effort to understand the developmental regulation of malate levels in the developing fruit of ‘golden delicious’, ma <dig> was not present in the major co-expression gene network due to its less significant  correlation with malate although the involvement of ma <dig> was clear  <cit> .
this is likely caused by the fact that the study focused on developing fruit of a single variety  of genotype mama that had malate levels equivalent to the normal to high acidity range during fruit development, i.e.
the truncation mutation effect of ma <dig> for low acidity in genotype group mama was absent.
the current models accounting for the almt <dig> mediated plant tolerance to soil aluminum toxicity in arabidopsis and rice comprise several common elements  <cit> .
these elements include environmental stimuli , a receptor  on the plasma membrane that interacts with al3+/h+ to initiate a signal transduction pathway , transcription factors , protein kinase  activating the tfs by phosphorylation, and responsive genes almt <dig>  abc transporters  and others  <cit> .
in the developing fruit of gd, several groups of genes functionally similar to these elements were co-enriched or co-suppressed in a malate dependent manner, including  <dig> c2h <dig> transcriptional factors,  <dig> protein kinases, and  <dig> receptor kinases for signaling  <cit> .
in addition, g <dig>  one of  <dig> c2h <dig> transcription factors, was annotated as a stop1-like protein by mercator  <cit> .
these had led to a speculation that a similar transcriptional regulation model for malate in developing fruit may exist although almt <dig> targets plasma membrane  <cit> .
in this study, some of these functionally similar components were also identified.
for example, within the  <dig> msags, we see m <dig> encoding a leucine-rich repeat  transmembrane receptor protein kinase for signaling, m <dig> a protein serine/threonine kinase, m <dig> an abc transporter, and of cause ma <dig> an almt <dig> like malate transporter .
interestingly, lemon-tree based regulator analysis identified the llr receptor kinase encoding gene m <dig> as a regulator of cluster  <dig> of  <dig> genes .
moreover, m <dig>  the most intramodular connected gene in module turquoise , encodes a protein kinase.
these data suggested that the almt <dig> mediated plant tolerance models somewhat be plausible for the transcriptional regulation of fruit acidity.
however, this study identified only one c2h <dig> zinc finger tf  in the  <dig> faags assigned to module red not highly correlated with acidity.
there is no evidence for c2h <dig> zinc finger tfs being part of the transcriptional regulation of ma <dig> mediated fruit acidity.
identification of regulator m <dig> by lemon-tree may have suggested an emerging mechanism for fruit acidity regulatory network mediated by ma <dig>  again, m <dig> is one of the msags  as well as a highly connected gene in module turquoise .
lemon-tree assigned m <dig> as a regulator to six clusters of  <dig> genes in total, including cluster  <dig> containing ma <dig> , cluster  <dig> containing g <dig> and m <dig> , and cluster  <dig> containing m <dig> , which are all among the  <dig> msags .
although g <dig> and m <dig> are of unknown function, m <dig> was speculated as a possible protein kinase component in the almt <dig> model aforementioned.
another relevant observation is that samples in these six clusters were all clustered alongside the line between the two genotype groups ma_ and mama, a pattern not observed for any other regulators, suggesting m <dig> is likely a critical regulator for fruit acidity.
m <dig> shows some sequence similarity to ein3-like genes, which are transcription factors important in ethylene response and signal transduction in plant  <cit> .
previous studies documented that the root growth inhibition by aluminum stress is mediated by ethylene in bean  <cit> , lotus japonicus  <cit> , and arabidopsis is  <cit> .
a more recent study reported that aluminum-induced malate efflux from wheat roots and tobacco cells transformed with taalmt <dig> was negatively regulated by ethylene although the molecular mechanisms remain unknown  <cit> .
these reports hint that it is possible for m <dig> to be a regulator of malate in apple fruit.
m <dig>  one of the  <dig> hub genes in module turquoise , is a regulator assigned to cluster  <dig> containing m <dig>  one of the  <dig> msags .
m <dig> encodes a cam binding protein similar to atiqd <dig> involved in calcium  signaling  <cit> .
m <dig>  another hub gene of module turquoise , encodes a protein similar to cam-like <dig>  of role in calcium signaling  <cit> .
among the  <dig> msags , m <dig> encodes a calcium regulated channel similar to atcngc  <dig>  it has been demonstrated that atcngc <dig> is localized to tonoplast  <cit>  as well as plasma membrane  <cit> , and to interact with all cam isoforms through its iq domain  <cit> .
in a study investigating carbohydrate metabolism in two apple genotypes that sharply differ in malate accumulation in a way similar to those between the two genotype groups ma_ and mama found that calcium contents were significant higher in high acid fruit than in low acid fruit  <cit> .
these observations suggested ca2+ signaling might also be an important direction for better understanding the regulation of malate levels in mature apple fruit.
interestingly, ca2+ signals have also been considered one of the possible signal transduction pathways for aluminum tolerance lately  <cit> .
CONCLUSIONS
the most relevant finding of this study is the identification of a weighted gene co-expression network analysis  module  of  <dig> genes that not only encompasses ma <dig> but also shows the highest modular correlation with acidity in mature fruit of diverse genetic background.
based on lemon-tree algorisms, mdp <dig>  which putatively encodes an ein3-like transcriptional factor, is likely the most important regulator of ma <dig> and its mediated gene network.
the two gene ontology biological process terms  and two mapmap sub-bins  that were significantly enriched in module turquoise implicated that photosynthesis related pathways are likely essential for acidity.
overall this study for the first time provides not only important insight into the ma1-mediated gene network controlling acidity in mature apple fruit of diverse genetic background, but also relevant clues for future biological validation, including the three putative regulators, mdp <dig>  mdp <dig> and mdp <dig>
supporting data
the data sets supporting the results of this article are available in the ncbi sequence read archive  repository  under accession srx <dig>
additional files
additional file 1: table s <dig>  fruit harvest date and weight.
additional file 2: table s <dig>  overview of rna-seq reads mapping.
additional file 3: table s <dig>  list of primers used in qrt-pcr.
additional file 4: table s <dig>  list of the  <dig> faags and their rpkm values and assignment in wgcna modules.
additional file 5: figure s <dig>  analysis of modules black, brown, blue and yellow.
module eigengene values across the  <dig> samples, including  <dig> in ma_ on left and  <dig> in mama on right.
samples are represented by the combination of a letter  and a number  .
correlation between module membership  and gene significance  for malate.
additional file 6: table s <dig>  list of the most significant acidity genes .
additional file 7: table s <dig>  list of intramodular hub genes in modules highly associated with acidity.
additional file 8: table s <dig>  regulator genes identified by lemon-tree.
additional file 9: figure s <dig>  regulator m <dig> and its other five assigned tight clusters.
elements and their contents, formats and messages are same as those noted in fig.
8a.
cluster  <dig> of  <dig> genes.
cluster  <dig> of  <dig> genes.
cluster  <dig> of  <dig> genes.
cluster  <dig> of  <dig> genes.
cluster  <dig> of  <dig> genes.
additional file 10: figure s <dig>  other regulators from modules turquoise and brown and their assigned tight clusters.
elements and their contents, formats and messages are same as those noted in fig.
8a.
regulator m <dig> and cluster  <dig> of  <dig> genes.
regulator m <dig> and cluster  <dig> of  <dig> genes.
regulator m <dig> and cluster  <dig> of  <dig> genes.
regulator m <dig> and cluster  <dig> of  <dig> genes.
regulator m <dig> and cluster  <dig> of  <dig> genes.
regulator m <dig> and cluster  <dig> of  <dig> genes.
regulator m <dig> and cluster  <dig> of  <dig> genes.
additional file 11: table s <dig>  enriched gene ontology  terms and mapman -bins in wgcna modules.
abbreviations
almt1aluminum-activated malate transporter1
degsdifferentially expressed genes
faagsfruit acidity associated genes
fdrfalse discovery rate
gdgolden delicious
ma1malic acid1
msagsmost significant acidity genes
qtlquantities trait locus
rna-seqmrna sequencing
rpkmreads per kilobase of exon model per million mapped reads
tatitratable acidity
competing interests
the authors declare that they have no competing interests.
authors’ contributions
yb and ld carried out fruit sampling, acidity and metabolite evaluation, rna isolation and rna-seq library construction and qpcr validation of rna-seq.
lc and gz participated in the design and coordination of the study.
kx conceived the study.
yb and kx analyzed the data and wrote the manuscript.
all authors read and approved the final manuscript.

 apple fruit acidity refers to the sensory intensity of tartness of fruit flesh tissues. the stronger the tartness taste, the higher the fruit acidity levels. chemically, fruit acidity can be quantified by measuring fruit juice ph and/or titratable acidity. it has been shown that several organic acids are present in mature apple fruit, including malic acid, quinic acid, citric acid and others, but malic acid consists of more than 90 % of the total, thereby largely determining fruit acidity. for dessert apples, fruit acidity could be grouped into three categories by ta: low, normal, and high, and categories high and low are not acceptable. this makes fruit acidity an essential quality component not only for determining the fate of existing varieties, but also for prompting routine evaluations of acidity levels in apple breeding as genotypes of fruit acidity beyond the acceptable normal range can make up 25–50 % in breeding populations. as a result, fruit acidity has long been an important subject area of investigations in apple genetics. apple fruit acidity is primarily controlled by the major gene or qtl on chromosome called malic acid alongside a significant qtl on chromosome and a few other qtls with relatively smaller effects. with regard to the allelic interactions at the major qtl ma, normal/high acid allele is largely dominant over the low acid allele. genotypes mama and mama have fruit of normal to high acidity while genotype mama of low acidity with little or none commercial value. to focus on the most important genetic factor, we and others have recently isolated the major qtl ma. the common findings were that the ma locus harbors two new members of the aluminum-activated malate transporter gene family, called ma and ma the studies further found that ma was expressed in significantly positive correlation with fruit acidity levels while the expression of ma was barely detectable in both high and low acid fruit, suggesting that it was gene ma rather than ma that was the very gene underlying ma. in addtion, a detailed analysis of the allele specific dna sequences of ma indicated that a single base mutation that would stop the protein translation process prematurely was almost completely associated with low acidity in a diverse set of apple germplasm studied, suggesting that the low acidity is caused by the malfunction of the ma protein due to the deduced truncation at the c-terminus. these latest findings have markedly increased our understanding on fruit acidity, but many questions remain to be answered. for example, how does the landscape of transcriptomes differ between genotypes ma_ and mama at maturity stage? if fruit acidity is governed by a gene network in which ma is the genetic determinant, what would the other possible members in the network? what are the potential biological process and/or regulatory mechanisms responsible for the contrast acidity levels between genotype groups ma_ and mama? the development of mrna sequencing technology that unlocks the power of high throughput of next generation sequencing has provided an ideal means to address these questions. since its inception, rna-seq has been rapidly adapted in transcriptomics studies in plants such as arabidopsis, grape, maize and rice. in apple, rna-seq based studies have recently been reported as well. more importantly, to begin resolving the low coverage issue of the current version of apple reference transcriptome, we have improved it with rna-seq reads from fruit of golden delicious, the source of the reference genome, which is available at the genome database for rosaceae. we have also used the improved reference transcriptome through rna-seq approach to construct a co-expression gene network associated with developmental regulation of malate levels varying from to  mg/g fw in developing fruit of ‘golden delicious’ of genotype mama. the objectives of this study are to address the questions aforementioned. to do so, we first sequenced rna-seq libraries representing transcriptomes of mature fruit from ten apple varieties of genotypes mama, and mama and mama, and then mapped the rna-seq reads against the improved apple reference transcriptome. using ma as a guide gene, a series of downstream analyses was conducted, leading to identification of weighted co-expression gene network modules significantly correlated with malate, most significant acidity genes, intramodular hub genes, regulator genes, enriched gene ontology terms and mapman sub-bins, and others. to the best of our knowledge, this is first report attempting to understand the ma1-mediated gene network regulating fruit acidity in apples of diverse genetic background. methods plant materials and fruit acidity quantification ten apple varieties of known genotype at the ma locus were chosen, including four of mama–britegold, sweet delicious, novosibirski sweet and pi, four of mama–fuji, rome beauty law, cox’s orange pippin and jonathan, and two of mama–empire and granny smith. genotypes mama and mama were jointly designated ma_ to represent a group that had at least one functional allele of ma the trees were budded onto rootstock p and grown in a research orchard of cornell university at geneva, new york. fruit of three replicates per variety and 8– fruit per replicate were harvested at maturity from to trees in fall. each fruit was cross-sectioned into two halves: one half were used for maturity evaluation and fruit juice extraction, and the other half were sliced and immediately frozen in liquid nitrogen for rna isolation and for quantitation of malate and other metabolites. the evaluation of fruit maturity and fruit acidity was conducted as previously described. fruit of cornell starch index – were considered matured and only matured fruit were used for analysis. two indicators of fruit acidity were measured: ph by a ph meter and titratable acidity by an autotitrator. fruit organic acids were extracted and measured using an agilent 7890a gc/5795c ms with the same configurations and settings as described previously. metabolites were identified by comparing fragmentation patterns with those in mass spectral libraries and quantified based on standard curves for each metabolite and the internal standard ribitol. rna-seq library construction and data analysis fruit rna isolation and rna-seq library construction were performed as previously described. briefly, total rna was isolated from 3 g of ground fruit flesh tissues and then treated with dnase i. for mrna isolation and strand specific rna-seq library construction, nebnext poly mrna magnetic isolation module and nebnext ultra directional rna library prep kit for illumina were used. the mrna was isolated from 5 μg of total rna and was fragmented at 94 °c for 10 min. the first strand cdna was reverse transcribed from the fragmented mrna with datp mix and second strand cdna was synthesized from the first strand cdna with dutp mix. the resulting double strand cdna was end-repaired, adaptor-ligated and size-selected, and then followed by the user enzyme digestion of the second strand cdna. the intact modified first strand cdna was pcr-enriched for 14– cycles to obtain the individual sequencing ready rna-seq libraries. these libraries were multiplexed in equal amount for single end 100-base sequencing in three lanes of hiseq at the cornell university biotechnology resource center. the ten apple varieties were sequenced in three biological replicates with one replicate per lane. illumina sequencing of the pooled rna-seq libraries generated fastq files of sequences with a total of million raw reads. the raw reads were fed into bowtie to remove bar codes and adapters and then aligned to the rrna reference sequences allowing up to three mismatches, which were downloaded from the silva rrna database. the rrna depleted reads were imported into clc genomics workbench v to trim low quality reads and/or bases using the quality limit of and the ambiguous limit of the resultant clean and high quality reads were mapped by clc genomics workbench against the improved apple reference transcriptome, which is available at url http://www.rosaceae.org/species/malus/malus_x_domestica/cu_rna_seq_genes of gdr. in addition to of the originally predicted genes or mdps, the reference transcriptome includes novel transcripts. for convenience, these novel transcripts will be referred to ‘genes’ and named ‘g######s’ as in ‘g101234’, and ‘mdp0000’ in the original gene ids will be abbreviated to ‘m’ hereafter. gene expression levels were calculated and represented by reads per kilobase of exon model per million mapped reads. genes of rpkm > were considered expressed according to a previous study. identification of fruit acidity associated genes fruit acidity associated genes were defined for those that were differentially expressed between genotype groups ma_ and mama and/or those that were expressed similarly as the guide gene ma to identify differentially expressed genes between the two genotype groups, the gene expression data rpkms were subjected to baggerly’s test. to control the false discovery rate, the original p values in baggerly’s test were adjusted for multiple testing using benjamini-hochberg correction. the cutoff for declaring a gene expressed differentially was pfdr <. however, any genes with an rpkm range < in the rna seq samples were removed from the list. to identify genes expressed similarly as ma two filters were applied: 1) the mean difference in rpkm between the two genotype groups was greater than ten, and 2) the fold change in mean rpkm between the two genotype groups was greater than. genes passed both filters were considered of similar expression as ma inferring fruit acidity associated co-expression gene network modules the co-expression gene network modules were inferred from the fruit acidity associated genes using weighted gene co-expression network analysis, an r package. before network inferring, the expression data were normalized by square root transformation. the automatic one-step network construction and module detection method with default settings were used, which include an unsigned type of topological overlap matrix, a power β of a minimal module size of and a branch merge cut height of. the module eigengene value was calculated and used to test the association of modules with acidity in the ten genotypes of samples. gene significance, total network connectivity, and module membership that is also known as eigengene-based intramodular connectivity, were calculated for each of the faags in the ten modules. the most significant module was visualized using cytoscape and were also analyzed using network analyzer, a cytoscape plugin. identification of regulator genes detection of regulator genes was performed using lemon-tree software, which is an improved and ‘one-stop shop’ version of the stochastic bayesian module network algorithm learning module networks, which requires run scripts on both matlab and perl. using expression data from the faags, ten independent gibbs sampler runs were conducted with lemon-tree. from these ten runs, we were able to generate a single set of tight independent clusters, and then select of at least ten genes per cluster for testing candidate regulators. the candidate regulators were genes annotated as transcription factors or signal transducers in go terms and/or mapman bins in the genes. we considered candidate regulators as regulators when their probabilistic scores associated with their assignment to any of the tight clusters were within top 1 %. functional annotation and enrichment analysis of wgcna modules we previously used both mapman bins and gene ontology terms to annotate the improved apple reference transcriptome. briefly, for annotation with mapman bins, several databases were blast-searched with the apple reference transcriptome of genes using mercator, a web-based search tool. the searched databases include: tair-arabidopsis tair proteins, ppap-swissprot/uniprot plant proteins, chlamy-jgi chlamy release augustus models, oryza-tigr rice proteins, kog-clusters of orthologous eucaryotic genes database, cdd- conserved domain database, and ipr-interpro scan. for go term annotations, relevant information for apple genes is available at gdr and was directly adapted. however, go terms for the new transcripts were obtained by blast2go program together with blastx search against the ncbi reference proteins database. enrichment analyses of the wgcna modules were all conducted by hypergeometric tests using clc gw against the same background of expressed genes. the hypergeometric test is similar to the unconditional gostats test, and the cutoff for significant enrichment is pfdr<. depending upon actual test, the background was defined by annotated genes in mapman bins, by in go biological process, or by in go molecular functions. qrt-pcr analysis qrt-pcr analyses were conducted with the same total rna samples used for the rna-seq library construction. two micrograms of total rna were reverse transcribed using the superscript iii rt module. the resulting first strand cdna was diluted by 5-fold and used as templates for qrt-pcr analysis on a lightcycler, where an apple actin gene was served as reference. the primer sequences for the reference gene and the eight target genes were listed. the actual qrt-pcr reactions and subsequent expression analysis were conducted similarly as described previously. statistical analyses statistical analyses, such as anova, student’s t test and tukey’s hsd test, were performed using jmp pro.  fruit metabolite profiling and acidity evaluation fruit metabolite profiling was conducted using gc-ms with three biological replicates in the ten apple varieties. a total of metabolites were quantified, including soluble sugars and seven organic acids. among the sugars, only sorbitol showed a significantly higher concentration in genotype group mama than in ma_. other sugars, especially the highly abundant fructose, sucrose and glucose, had similar levels across these varieties. among the seven organic acids, malate was the most abundant while quinate and maleate were distantly followed as the second and third abundant, respectively. malate concentrations varied significantly in the ten apple varieties, and were more than 3-fold higher in genotype group ma_ than in mama. the concentrations of maleate, succinate, fumarate and citrate were also significantly higher in genotype group ma_ than in mama, but they were at much lower levels. however, quinate showed no significant difference between groups ma_ and mama. the trend of fruit juice titratable acidity was similar to that of malate in the ten varieties while an expected reverse trend was observed for fruit juice ph.fig. evaluation of fruit acidity and causal organic acids. standard deviations were shown with the error bars. a ta and malate concentrations measured in the ten apple varieties, which were abbreviated as the following: b britegold, n novosibirski sweet, p pi s sweet delicious, c cox’s orange pippin, e empire, f fuji, g granny smith, j jonathan, r rome beauty law. the first four are of genotype mama whereas the last six of genotype ma_. b mean ta and malate concentrations in genotype groups ma_ and mama. c concentrations of succinate, fumarate, citrate and maleate in genotype groups ma_ and mama. d fruit ph readings in the ten apples transcriptome analysis in genotype groups ma_ and mama rna-seq data analysis was performed using the improved apple reference transcriptome. after removing reads derived from rrna and/or of low quality, the reads for rna-seq mapping were million in total. overall the total mapped and uniquely mapped reads were million and million, respectively. the mean mapped reads per sample were  ±  million in total and  ±  million in unique. however, the sample for replicate i of granny smith was excluded from the downstream analysis as  % of its raw reads was removed by the rrna and quality filters. the maximum number of expressed genes was in this study, a tally counting all genes of rpkm > in at least one of the samples. however, measuring by the threshold in group means, the expressed genes were. of these, were expressed in both genotype groups, and and were specifically in ma_ and mama, respectively.fig. summary of rna-seq data analysis. a venn diagram representation of the number of genes expressed in genotype groups ma_ and mama. the cutoff is group mean rpkm >. b mean number of total and uniquely mapped reads per sample in genotype groups ma_ and mama. standard deviations were shown with the error bars the total mapped reads of the expressed genes were  ±  million per sample in the ma_ group and  ±  million per sample in mama. the uniquely mapped reads per sample in the ma_ and mama groups were  ±  million and  ±  million, respectively. identification of malate associated co-expression gene network modules using wgcna transcriptome-wide comparison of the expressed genes revealed fruit acidity associated genes, including of differentially expressed genes between groups ma_ and mama, and expressed similarly as the guide gene ma. using wgcna, a co-expression gene network was constructed from the faags. within the network, the total network connectivity of genes showed a significant correlation with the gene significance for malate, suggesting that more connected genes would have higher gene significance for acidity in the faags.fig. weighted gene co-expression network analysis of fruit acidity associated genes in the ma_ and mama genotype groups. a hierarchical cluster tree showing ten modules of co-expressed genes. each of the faags is represented by a leaf in the tree while each of the ten modules by a major tree branch. the lower panel shows modules in designated colors. b correlation between gene total network connectivity and absolute gene significance for malate in the whole co-expression gene network. c module-fruit acidity correlations and corresponding p-values. the left panel shows ten modules and the number of their member genes. the right panel is a color scale for module trait correlation from − to 1 further wgcna analysis identified ten network modules in the co-expression network, designated turquoise, black, brown, pink, magenta, red, green, blue, yellow, and grey. investigating the relationships between module eigengene and acidity uncovered that the correlation coefficients varied widely from − to in malate and from − to in ta. the eigengenes of five modules turquoise, black, brown, blue and yellow showed significant correlations with malate, suggesting these five modules had greater relevance in fruit acidity although blue and yellow were negatively correlated. the five modules comprised genes with module turquoise being the largest of genes. notably, the guide gene ma was assigned to turquoise, which also had the highest modular correlation with acidity. inspecting the correlation between the module memberships and gene significance revealed that all modules were significant except blue. again, module turquoise had the most significant correlation between mm and gs, indicating genes of higher mm values were more likely of greater gs in turquoise than in others.fig. analysis of module turquoise. a module eigengene values across the samples, including in ma_ on left and in mama on right. samples are represented by the combination of a letter and a number. b correlation between module membership and gene significance for malate identification of the most significant acidity genes the faags had a range of absolute gene significance for malate from to. we define that the most significant acidity genes are those of gs for malate ≥, the observed gs of the guide gene ma this allowed to identify msags, including three transporters, two receptor kinases, one ethylene-insensitive -like transcription factor, one pyruvatedehydrogenase complex component e, one endo- -beta-xylanase/ hydrolase, one glutamine synthetase cytosolic isozyme, one photosystem ii subunit r, and seven of unknown function. the msags were distributed in three modules with in turquoise, three in black and another three in blue, further indicating module turquoise is of greater relevance in acidity.fig. the most significant acidity genes. a expression profile of the msags. each row represents a gene as listed on right. each column stands for a sample as specified by sample names at lower panel. clustering of genes and samples was shown with distance on left and top panels, respectively. the expression of genes is color coded from low through mean and through high. the bright yellow color across the samples for ma and m indicate the expression of these two genes were all close to or exceeded the high end of the scale at bottom. b distribution of the msags in modules. modules with zero genes are not shown. note that genes of ids beginning with letter ‘g’ are referred to the novel transcripts, and the apple reference gene ids are abbreviated identification of intramodular hub genes and evaluation of their total network connectivity changes in genotype groups ma_ and mama to identify intramodular hub genes, we similarly used the guide gene ma which had an mm of, ranking 11th in module turquoise. based on this information, we defined top genes in mm values from each of the five modules as intramodular hub genes, leading to identification of hub genes. to evaluate how the hub genes behaved in genotype groups ma_ and in mama, their respective total network connectivity were calculated and compared with the overall group of both genotypes. most of the hub genes showed a large increase in ktotal in mama while a marked reduction in ma_, suggesting a major pattern of negative correlation between the hub genes’ ktotal changes and acidity levels. this pattern was true in modules blue and yellow, largely true in brown and black, but largely untrue in turquoise, seemingly providing another means depicting the module-acidity relationship shown early.fig. intramodular hub genes and their total network connectivity changes in genotype groups ma_ and mama when compared with the combined groups of both genotypes. shown are modules black, blue, brown, turquoise and yellow there were three additional patterns in ktotal changes: the first was observed only on m in turquoise, i.e. a reverse trend by having an increased ktotal in ma_ and a decreased ktotal in mama. m encodes a leucine-rich repeat receptor kinase. the second was represented by ma i.e. the ktotal was decreased in both ma_ and mama, suggesting reduced ktotal of ma in mama might have important ramifications for low acidity. this pattern was also observed on m binding protein involved in calcium signaling), m and m in turquoise, m m and g in black, and m and m in brown. the third pattern was characterized by an increased ktotal in both ma_ and mama, which was observed on m, m and g in turquoise, and g and m in black. to graphically view the intramodular connectivity of these hub genes, the hub genes along with other co-expressed genes in module turquoise was presented using cytoscape. it showed that the hub genes were of the most edges from to and were in the module network core.fig. graphic view of module turquoise. a cytoscape representation of co-expressed genes with edge weight ≥. the edge number of the genes ranges from to, which is color coded by the scale at top left from blue through red. some important genes are noted, including intramodular hub genes m m m m m m m m m g and ma, two regulator genes m and m and one hub and regulator gene m. b distributions of edges. the hub genes are among those of most edges of weight ≥ identification of regulator genes regulator genes are not readily identifiable from the wgcna co-expression network modules as they are non-directed. to identify regulator genes, we used lemon-tree software first to generated tight clusters covering of the faags through ten gibbs sampler runs, and then to examine the candidate regulators with the clusters. to statistically validate the results obtained from the candidate regulators, another randomly selected genes were used as control and examined in parallel with the candidate regulators. as expected, the probabilistic scores for genes assigned as a regulator for the tight clusters showed significant difference between the candidate regulators and ck. at the top 1 % of p. scores, of the candidate regulators were identified as regulators and assigned to of tight clusters. of the regulators, five were from module turquoise, two from brown, and five from modules green, red and pink which were irrelevant for acidity. below are brief descriptions of three regulator genes, which were among the msags or among the most connected intramodular hub genes in module turquoise.fig. regulator genes and their assigned tight clusters. a regulator m and cluster of genes. the expression of genes is color coded from low through high. each row stands for a gene and each column for a sample named at the bottom. the regulator is assigned based on the hierarchical tree on top, which indicates how samples were clustered together with the red vertical lines. the short red vertical line in the sample list shows where the two primary clades diverge while the short green vertical line marks where the secondary clades depart within a primary clade. the genotype names in red indicate they are not in agreement with one of the two primary clusters ma_ or mama. b regulator m and cluster of genes. c regulator m and cluster of genes. elements in and and their contents, formats and messages are the same as those in the first is m encoding an ein3-like transcription factor. m was assigned to six tight clusters and with p. scores –. notably, lemon-tree clustered ma into cluster to which m was assigned as a regulator with the highest p. score, suggesting m might potential regulate ma m was one of the msags and of high mm value in module turquoise. in the six m regulated clusters, the fruit samples were all clustered into two primary clades, which were nearly perfectly alongside the line dividing the two genotypes ma_ and mama. such clustering pattern of fruit samples were not observed in any of the remainder clusters regulated by the other regulator genes, further supporting the potential regulatory role of m in fruit acidity. the second is m encoding a receptor-like protein kinase. m was assigned to cluster with p. score. similar to m m is also one of the msags that had a high mm value in module turquoise. the third is m one of the intramodular hub genes in module turquoise. m was assigned to cluster with p. score. as aforementioned, m encodes an atiqd2-like cam binding protein involved in calcium signaling and behaved similarly as ma in ktotal changes in genotype groups ma_ and mama. detection of enriched go terms and mapman bins in network modules to detect enriched go terms and mapman bins in the five significant modules, a series of hypergeometric tests were performed against the background of expressed genes. overall, eight go terms and ten mapman sub-bins were enriched in modules turquoise, brown, blue and yellow, but not in black, and no terms were repressed. due to multifaceted roles of genes and due to both go and mapman systems were used, some of the enriched terms had the same or mostly overlapping member genes. in total, the enriched terms covered non-redundant genes, including in module turquoise, in brown, in blue and in yellow.table go terms and/or mapman bines significantly enriched in wgcna modules significantly correlated with fruit aciditya athe go terms and/or mapman bins that have the same or mostly overlapping member genes in the same module are highlighted with same color bbai et al. in turquoise, two go terms and two mapman sub-bins were enriched, which are all related to photosynthesis. however, go: and sub-bin. had the same five member genes. notably, m a member in go: is one of the msags. in brown, three go terms and one mapman sub-bin. were enriched. the three go terms were highly related as the four member genes in go: and go: were identical and were also among the six member genes in go:. m a common member gene in these go terms, is one of the hub genes in brown. in blue, we detected three go terms and five mapman sub-bins, the most in this study. but the ten member genes in go: were fully covered by the genes in sub-bin.. similarly, the two genes in sub-bin. and three genes in go: were covered by the genes in go:. sub-bins,. and were unique. remarkably, four of the hub genes in blue were present, i.e. m and m in go: and m and m in sub-bin. in yellow, two mapman sub-bins and. were enriched, and m from sub-bin is one of the hub genes in yellow. qrt-pcr confirmation of gene expression to evaluate if and how the rpkm values reflected the gene expression levels, a set of eight genes were analyzed using qrt-pcr. the eight genes include ma the ein3-like regulator m another most significant gene for acidity m encoding a protein serine/threonine kinase, and five others. the data not only confirmed that the relative expressions of the eight genes in qrt-pcr were significantly correlated with their rpkm values in rna-seq, but also confirmed their differential expression between the two genotypes groups ma_ and mama.fig. expression confirmation of eight selected genes using qrt-pcr. a-h the normalized expression of target genes relative to a control gene in qrt-pcr was shown in light grey, and their corresponding rpkm values from rna-seq were in black. the correlation coefficient and associated p value were shown accordingly. standard deviations were shown with the error bars discussion identification of individual fruit acidity associated genes in this study, ten apple varieties were carefully chosen to represent the two genotype groups ma_ and mama of contrast fruit acidity levels. these varieties were grafted on the same rootstock p and grown in the same orchard with the same management routines to curtail the environmental impacts on gene expression. however, since the apples are of diverse genetic background and matured at varying dates, fruit transcriptome at maturity were inevitably influenced by these genetic and environmental variables. to minimize such biases, three replicates per variety and at least four varieties per genotype group were sampled and used for rna-seq analysis. this strategy appeared to be appropriate and adequate for the defined objectives in this study as we identified fruit acidity associated genes. in developing fruit of ‘golden delicious ’ of varying malate levels, differentially expressed genes were identified. in comparison of the two lists of genes, of the faags were not present in the degs, suggesting most genes relevant for fruit acidity variation across diverse genetic background were largely different from those important for malate variations in the developing fruit of gd. overall, we found genes expressed in the ma_ group and genes in the mama group. as there were genes expressed in both groups, the total number of expressed genes was. varying rpkm thresholds were employed in published studies, we view rpkm > used previously fit this set of rna-seq data well. of the expressed genes, were regarded faags. given the significant correlation between rpkm in rna-seq and relative expression in qrt-pcr in the eight genes tested, we consider that the majority, if not all, of the genes relevant for fruit acidity variation at maturity in this study were included in the genes. co-expression gene network modules and go terms and mapman bins relevant for fruit acidity wgcna is one of the programs used widely for inferring co-expression network modules. using wgcna, we identified five modules highly associated with fruit acidity. hypergeometric tests identified eight go terms and ten mapman bins that were enriched in four of the five modules. in the most significant module turquoise that includes ma the two enriched go terms and two mapman sub-bins are all related to photosynthesis. in the developing fruit, the two go terms and the two mapman sub-bins were also enriched, and of the member genes in the four enriched terms were identified in both studies. this suggested that fruit photosynthesis could likely be a part of the mechanism for the contrast malate levels between the two genotype groups at maturity. it was proposed that the relative more active photosynthesis in young fruit may facilitate malate biosynthesis, thereby contributing to higher malate levels in young fruit than in maturing fruit. however, chloroplasts are only present in the hypodermal green tissues and the inner perivascular tissues in mature apple fruit, and are essentially absent in parenchyma cells that constitute the major portion of flesh. how enhanced photosynthesis activities in a fraction of cells would lead to increased overall malate levels in fruit remains to be understood. in module brown, the enrichment of mapman sub-bin. was unique to this study. however, the three enriched go terms were essentially a reflection of mapman sub-bin, which was enriched in negative malate dependent manner in developing fruit. in module blue, the enrichment of mapman sub-bins,. and were identified in this study only. but go:, and mapman sub-bins. and. were also enriched in developing fruit. in module yellow, the enrichment of two mapman sub-bins and. was unique in these diverse mature fruit. overall, enriched terms, including four from turquoise, three from brown, and five from blue, were also enriched in developing fruit of ‘golden delicious’, whereas the remainder six terms were enriched specifically in this study. this suggested that common and unique processes relevant for malate variations in both the diverse mature fruit and the developing fruit of ‘golden delicious’ exist. however, of the genes in blue and all of the six genes in brown were not found to be relevant in developing fruit although of the genes related to photosynthesis in turquoise were, suggesting the common processes may operate but most likely with their unique gene sets. transcriptional regulation of fruit acidity mediated by ma1 in an ongoing survey of the usda malus repository at geneva, new york, we identified diverse apple accessions of genotype ma1ma using marker caps that can detect the stop codon leading snp in ma and found that each of them had fruit ph >, a typical low fruit acidity characteristic in genotype group mama. these data, together with the previous studies, suggested us that ma a member of the almt gene family, is the genetic determinant underlying the ma locus controlling acidity levels in mature apple fruit although the functional genetic complementation data are not available. ma is thought to exert its effect on fruit acidity by both the truncation leading mutation and expression levels, where the mutation plays greater role. in this study, ma showed not only high gene significance for malate, but also high intramodular connectivity in module turquoise, further supporting this notion. however, in our previous effort to understand the developmental regulation of malate levels in the developing fruit of ‘golden delicious’, ma was not present in the major co-expression gene network due to its less significant correlation with malate although the involvement of ma was clear. this is likely caused by the fact that the study focused on developing fruit of a single variety of genotype mama that had malate levels equivalent to the normal to high acidity range during fruit development, i.e. the truncation mutation effect of ma for low acidity in genotype group mama was absent. the current models accounting for the almt mediated plant tolerance to soil aluminum toxicity in arabidopsis and rice comprise several common elements. these elements include environmental stimuli, a receptor on the plasma membrane that interacts with al3+/h+ to initiate a signal transduction pathway, transcription factors, protein kinase activating the tfs by phosphorylation, and responsive genes almt abc transporters and others. in the developing fruit of gd, several groups of genes functionally similar to these elements were co-enriched or co-suppressed in a malate dependent manner, including c2h transcriptional factors, protein kinases, and receptor kinases for signaling. in addition, g one of c2h transcription factors, was annotated as a stop1-like protein by mercator. these had led to a speculation that a similar transcriptional regulation model for malate in developing fruit may exist although almt targets plasma membrane. in this study, some of these functionally similar components were also identified. for example, within the msags, we see m encoding a leucine-rich repeat transmembrane receptor protein kinase for signaling, m a protein serine/threonine kinase, m an abc transporter, and of cause ma an almt like malate transporter. interestingly, lemon-tree based regulator analysis identified the llr receptor kinase encoding gene m as a regulator of cluster of genes. moreover, m the most intramodular connected gene in module turquoise, encodes a protein kinase. these data suggested that the almt mediated plant tolerance models somewhat be plausible for the transcriptional regulation of fruit acidity. however, this study identified only one c2h zinc finger tf in the faags assigned to module red not highly correlated with acidity. there is no evidence for c2h zinc finger tfs being part of the transcriptional regulation of ma mediated fruit acidity. identification of regulator m by lemon-tree may have suggested an emerging mechanism for fruit acidity regulatory network mediated by ma again, m is one of the msags as well as a highly connected gene in module turquoise. lemon-tree assigned m as a regulator to six clusters of genes in total, including cluster containing ma, cluster containing g and m, and cluster containing m, which are all among the msags. although g and m are of unknown function, m was speculated as a possible protein kinase component in the almt model aforementioned. another relevant observation is that samples in these six clusters were all clustered alongside the line between the two genotype groups ma_ and mama, a pattern not observed for any other regulators, suggesting m is likely a critical regulator for fruit acidity. m shows some sequence similarity to ein3-like genes, which are transcription factors important in ethylene response and signal transduction in plant. previous studies documented that the root growth inhibition by aluminum stress is mediated by ethylene in bean, lotus japonicus, and arabidopsis is. a more recent study reported that aluminum-induced malate efflux from wheat roots and tobacco cells transformed with taalmt was negatively regulated by ethylene although the molecular mechanisms remain unknown. these reports hint that it is possible for m to be a regulator of malate in apple fruit. m one of the hub genes in module turquoise, is a regulator assigned to cluster containing m one of the msags. m encodes a cam binding protein similar to atiqd involved in calcium signaling. m another hub gene of module turquoise, encodes a protein similar to cam-like of role in calcium signaling. among the msags, m encodes a calcium regulated channel similar to atcngc it has been demonstrated that atcngc is localized to tonoplast as well as plasma membrane, and to interact with all cam isoforms through its iq domain. in a study investigating carbohydrate metabolism in two apple genotypes that sharply differ in malate accumulation in a way similar to those between the two genotype groups ma_ and mama found that calcium contents were significant higher in high acid fruit than in low acid fruit. these observations suggested ca2+ signaling might also be an important direction for better understanding the regulation of malate levels in mature apple fruit. interestingly, ca2+ signals have also been considered one of the possible signal transduction pathways for aluminum tolerance lately.  the most relevant finding of this study is the identification of a weighted gene co-expression network analysis module of genes that not only encompasses ma but also shows the highest modular correlation with acidity in mature fruit of diverse genetic background. based on lemon-tree algorisms, mdp which putatively encodes an ein3-like transcriptional factor, is likely the most important regulator of ma and its mediated gene network. the two gene ontology biological process terms and two mapmap sub-bins that were significantly enriched in module turquoise implicated that photosynthesis related pathways are likely essential for acidity. overall this study for the first time provides not only important insight into the ma1-mediated gene network controlling acidity in mature apple fruit of diverse genetic background, but also relevant clues for future biological validation, including the three putative regulators, mdp mdp and mdp supporting data the data sets supporting the results of this article are available in the ncbi sequence read archive repository under accession srx additional files additional file 1: table s fruit harvest date and weight. additional file 2: table s overview of rna-seq reads mapping. additional file 3: table s list of primers used in qrt-pcr. additional file 4: table s list of the faags and their rpkm values and assignment in wgcna modules. additional file 5: figure s analysis of modules black, brown, blue and yellow. module eigengene values across the samples, including in ma_ on left and in mama on right. samples are represented by the combination of a letter and a number. correlation between module membership and gene significance for malate. additional file 6: table s list of the most significant acidity genes. additional file 7: table s list of intramodular hub genes in modules highly associated with acidity. additional file 8: table s regulator genes identified by lemon-tree. additional file 9: figure s regulator m and its other five assigned tight clusters. elements and their contents, formats and messages are same as those noted in fig. 8a. cluster of genes. cluster of genes. cluster of genes. cluster of genes. cluster of genes. additional file 10: figure s other regulators from modules turquoise and brown and their assigned tight clusters. elements and their contents, formats and messages are same as those noted in fig. 8a. regulator m and cluster of genes. regulator m and cluster of genes. regulator m and cluster of genes. regulator m and cluster of genes. regulator m and cluster of genes. regulator m and cluster of genes. regulator m and cluster of genes. additional file 11: table s enriched gene ontology terms and mapman -bins in wgcna modules. abbreviations almt1aluminum-activated malate transporter1 degsdifferentially expressed genes faagsfruit acidity associated genes fdrfalse discovery rate gdgolden delicious ma1malic acid1 msagsmost significant acidity genes qtlquantities trait locus rna-seqmrna sequencing rpkmreads per kilobase of exon model per million mapped reads tatitratable acidity competing interests the authors declare that they have no competing interests. authors’ contributions yb and ld carried out fruit sampling, acidity and metabolite evaluation, rna isolation and rna-seq library construction and qpcr validation of rna-seq. lc and gz participated in the design and coordination of the study. kx conceived the study. yb and kx analyzed the data and wrote the manuscript. all authors read and approved the final manuscript. 
guided by ma we analyzed the transcriptomes of mature fruit of contrasting acidity from six apple accessions of genotype ma_ and four of mama using rna-seq and identified fruit acidity associated genes, among which were most significant acidity genes. we also identified intramodular hub genes from each of the five modules and enriched gene ontology terms and mapman sub-bines, including two go terms and two mapmap sub-bins related to photosynthesis in module turquoise. the most relevant finding of this study is the identification of the msags, intramodular hub genes, enriched photosynthesis related processes, and regulator genes in a wgcna module turquoise that not only encompasses ma but also shows the highest modular correlation with acidity. several organic acids are present in apple fruit, but malic acid is predominant and determines fruit acidity. we hypothesize that fruit acidity is governed by a gene network in which ma is key member. malusfruit acidityma1aluminum-activated malate transporter gene networkissue-copyright-statement© the author 2015 using lemon-tree algorithms, we identified regulator genes of probabilistic scores –, including mdp, mdp and mdp of greater interest for being one of the msags or one of the intramodular hub genes in turquoise, and/or a regulator to the cluster containing ma network inferring using weighted gene co-expression network analysis revealed five co-expression gene network modules of significant correlation with malate. overall, this study provides important insight into the ma1-mediated gene network controlling acidity in mature apple fruit of diverse genetic background. the trait is largely controlled by the malic acid locus, underpinning which ma that putatively encodes a vacuolar aluminum-activated malate transporter -like protein is a strong candidate gene. the goal of this study is to identify the gene network and the potential mechanisms through which the network operates.
 acidity is a major contributor to fruit quality. several organic acids are present in apple fruit, but malic acid is predominant and determines fruit acidity. the trait is largely controlled by the malic acid locus, underpinning which ma that putatively encodes a vacuolar aluminum-activated malate transporter -like protein is a strong candidate gene. we hypothesize that fruit acidity is governed by a gene network in which ma is key member. the goal of this study is to identify the gene network and the potential mechanisms through which the network operates.  guided by ma we analyzed the transcriptomes of mature fruit of contrasting acidity from six apple accessions of genotype ma_ and four of mama using rna-seq and identified fruit acidity associated genes, among which were most significant acidity genes. network inferring using weighted gene co-expression network analysis revealed five co-expression gene network modules of significant correlation with malate. of these, the ma containing module of genes showed the highest correlation. we also identified intramodular hub genes from each of the five modules and enriched gene ontology terms and mapman sub-bines, including two go terms and two mapmap sub-bins related to photosynthesis in module turquoise. using lemon-tree algorithms, we identified regulator genes of probabilistic scores –, including mdp, mdp and mdp of greater interest for being one of the msags or one of the intramodular hub genes in turquoise, and/or a regulator to the cluster containing ma  the most relevant finding of this study is the identification of the msags, intramodular hub genes, enriched photosynthesis related processes, and regulator genes in a wgcna module turquoise that not only encompasses ma but also shows the highest modular correlation with acidity. overall, this study provides important insight into the ma1-mediated gene network controlling acidity in mature apple fruit of diverse genetic background. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords malusfruit acidityma1aluminum-activated malate transporter gene networkissue-copyright-statement© the author 2015 