BACKGROUND
performance is related to the subtle variation in gene expression and this relationship differs among individuals  <cit> .
in well-defined families, the first level of variation comes from genetic segregation and recombination of chromosomes.
as a result of segregation and chromosomal recombination, each individual has different genetic makeup.
upon a given genetic background, genetic potential carried on dna can only be realized when the genes are expressed.
at the whole genome level, expression of each gene is affected by its genetic regulatory element as well as trans-acting factors including the impact of environment.
a composite of genes, transcriptional regulation, post-transcriptional modification and regulation, translational regulation and post-translational modification and regulation, along with environmental impact and genotype-environment interactions eventually determines the phenotypes of individuals.
when considered at the whole genome level, expression of tens of thousands of genes and combination of these genes make the variation of performance traits extremely complex with huge variability.
the task of modern agricultural genomics is to gain understanding of such variations and their relationship in determination of production and performance traits.
traditionally, genetic and molecular biological studies are conducted to dissect these variables at different levels.
for instance, the effect of alleles can be dissected through genetic and qtl mapping analysis  <cit> .
gene expression can be analyzed using high throughput methodologies such as microarrays and rna-seq analysis  <cit> .
various epigenetic regulations have also been studied to understand the differences in gene expression with similar genetic background.
such analyses have been very powerful in determination of genetic and epigenetic factors affecting performance and production traits  <cit> .
however, performance and production traits are often highly complex and the outcome of agricultural operations is affected by variations at all levels.
for example, genetic background is very important because disease resistance genes allow the organism to survive the serious infections  <cit> .
in most cases, disease resistance genes have been studied through genetic linkage and qtl analyses that allow the identification of genomic regions containing disease resistance genes to be identified.
even with the most powerful molecular approaches, analysis of complex traits such as disease resistance can be extremely challenging.
in  <dig>  michelmore et al.
developed a method called bulked segregant analysis  to study disease resistance in plants  <cit> .
the basic idea of bsa was that phenotypic extremes should have drastic differences in genotypes.
when samples are selected from phenotypic extremes, say the best and the worst performers, and their genotypes are analyzed in bulk, a correlation of genotypes with phenotypes can be attained.
in other word, the variation among individuals may be quite subtle and difficult to detect; however, the pooled samples  of the phenotypic extremes should pose a strong contrast in their genotypes at the genomic location linked to the trait.
bsa has been used in numerous studies to associate phenotypes with related genomic locations.
bsa has been evolving along with various types of molecular markers including restriction fragment length polymorphisms   <cit> , random amplified polymorphic dnas   <cit> , simple sequence repeats   <cit> , amplified fragment length polymorphisms   <cit>  and single nucleotide polymorphisms   <cit> .
with the development of next-generation sequencing  technologies, bsa was first enhanced by the application of sequence-based markers such as restriction-site associated dna  markers  <cit>  and whole genome sequencing  <cit> .
in recent years, the application of rna-seq  <cit>  has allowed rapid and comprehensive understanding of transcriptome level of variations.
coupling of bsa with rna-seq has led to the development of bulked segregant rna-seq , and it has been successfully applied in plants  <cit> , but not yet demonstrated in animals.
apparently, bsr-seq possesses the advantages of both bsa and rna-seq, with the high throughput for deep coverage of the transcriptome as well as the strong ability to detect genetic differences underlining the traits.
such a technique is best suited to organisms with high fecundity such as many species of fish.
catfish is the major aquaculture species in the united states, accounting for over 60% of all us aquaculture production.
the two major cultivated catfish species are channel catfish  and blue catfish .
an inter-specific hybrid  has been popular for aquaculture because of strong heterosis  <cit> .
not only is the interspecific hybrid is popular for aquaculture, it is also a superior system to study disease resistance because of their strong phenotypic difference.
blue catfish is extremely resistant against esc disease while channel catfish is relatively susceptible.
genetic linkage analysis of f <dig> generation of the interspecific hybrids  <cit>  would allow identification of disease resistance/susceptibility genes.
in this study, we take advantage of bsr-seq for the analysis of disease resistant genes using the f <dig> generation backcross progenies  of the interspecific hybrids.
here we demonstrate that bsr-seq is capable of 1) revealing differentially expressed genes; 2) revealing positional candidates containing genes related to disease resistance after mapping snps on the whole genome; and 3) revealing allele-specific expression after bacterial infection.
RESULTS
sequence assembly and analysis
rna-seq was conducted using illumina sequencing with three pooled samples of resistant fish, susceptible fish, and control fish.
each pooled sample contained equal amount of rna collected from  <dig> individuals.
each pooled sample was barcoded such that reads from each of the three samples can be traced and analyzed separately.
a total of over  <dig> million reads were generated with  <dig> million reads from the resistant fish,  <dig> million from the susceptible fish, and  <dig> million from the control fish.
after quality trimming, a total of  <dig> million reads were carried forward for analysis, with an average read length of  <dig>  bp .
raw read data are archived at the ncbi sequence read archive  under accession srp <dig>
the reads after quality trimming were assembled de novo into  <dig>  non-redundant contigs  with a n <dig> contig length of  <dig>  bp and an average contig length of 825 bp.
of the assembled contigs, over  <dig>  had a length of over  <dig>  bp.
the assembled contigs were analyzed by blast searches to determine the gene identities of the contigs.
of the  <dig>  contigs,  <dig>  had hits to  <dig>  unigenes, of which  <dig>  were known genes and  <dig>  were hypothetical genes .
de novo
assembly of the catfish liver transcriptome with infection of esc generated by illumina sequencing and assembled with trinity
contigs were used as queries for blast searches against zebrafish and nr protein databases.
the cut-off value for gene identity was set at e-value ≤ 1e- <dig>
differentially expressed genes after infection
differentially expressed genes after infection in resistant fish and susceptible fish were determined by comparing their expression levels in rpkm  with that of the control group.
as summarized in table  <dig>  a total of  <dig> genes were differentially expressed in the resistant group as compared with the control group.
a significant fraction of the differentially expressed genes in resistant fish were highly regulated, with  <dig> genes  being up- or down-regulated 10-fold or more, and  <dig> additional genes up- or down-regulated 5– <dig> fold .
edwardsiella ictaluri
although a relatively small number of genes were differentially expressed in resistant fish, a large number of genes were differentially expressed in susceptible fish, with a total of  <dig>  genes being differentially expressed .
not only the number of differentially expressed genes was drastically more in susceptible fish, the number of highly regulated genes was also much greater in susceptible fish, with  <dig> genes being up- or down-regulated  <dig> folds or more, and additional  <dig> genes were up- or down-regulated 5– <dig> folds .
comparison of gene expression in resistant fish and susceptible fish after infection
although a small number of genes were differentially expressed in resistant fish when being compared with the control fish, a large number of genes exhibited differential expression in resistant fish when being compared with the susceptible fish.
a total of  <dig>  genes were differentially expressed in resistant and susceptible fish, with  <dig> genes expressed significantly higher in resistant fish than in susceptible fish and  <dig> genes expressed significantly lower in resistant fish than in susceptible fish .
of the genes expressed significantly higher in resistant fish,  <dig> genes were expressed over  <dig> times more than in susceptible fish;  <dig> were expressed 50– <dig> times more than in susceptible fish;  <dig> were expressed 10– <dig> times more than in susceptible fish .
of the genes expressed significantly lower in resistant fish,  <dig> genes were expressed over  <dig> times less in resistant fish than in susceptible fish;  <dig> were expressed 50– <dig> times less in resistant fish than in susceptible fish;  <dig> were expressed 10– <dig> times less in resistant fish than in susceptible fish .
edwardsiella ictaluri
these genes are listed in additional file 1: table s <dig>
identification of snps and significant snps
snps were identified by alignment of short reads to the reference assembly of the rna-seq.
in order to be qualified for snps, at least  <dig> reads were required for each group  and a total of minor allele reads count must be greater than  <dig> among three groups.
using software popoolation  <dig>  a total of  <dig>  snps were identified.
in order to be sure the identified snps were properly identified, a second software package, varscan  <dig>  <cit>  was also used.
a total of  <dig>  snps were identified.
the difference was caused by differences in cutoff of quality scores of sequence reads by the two programs.
although the total numbers of snps identified by the two software differed by approximately 6%, the vast majority  of snps identified by the two software were identical.
these snps were located within  <dig>  contigs.
in order to determine snps with significant difference in allele frequencies between the resistant and susceptible fish , fisher’s exact test was initially performed.
as shown in table  <dig>   <dig>  snps were identified as significant snps.
these significant snps were located in  <dig>  contigs.
of the  <dig>  contigs,  <dig>  had significant hits to known genes, and the remaining probably represented contigs assembled from 5′- and 3′-untranslated regions or from long non-coding rnas.
the  <dig>  contigs with hits to known genes represented  <dig>  unique genes harboring  <dig>  snps .
bulk frequency ratios
although significant snps identified through rna-seq analysis reflect the final ratios of different alleles at the rna level in the two bulked samples, statistical analysis for significant snps using fisher’s exact test was only the first step for screening snps that may be significantly associated with the trait.
in order to compare the snp allele frequencies more directly, bulk frequency ratios  were generated from the rna-seq data between the two bulks, the resistant fish and the susceptible fish.
the bfr of genes were determined by the maximum bfr of the significant snps generated from the fisher’s exact test located in this gene.
as shown in figure  <dig>  large proportion of genes containing the significant snps had a bfr of at least  <dig>  a total of  <dig> genes had a bfr equal or greater than  <dig>  among these genes,  <dig>  genes had a bfr of 4–16;  <dig> genes had a bfr over 16; and  <dig> genes had a bfr over  <dig> .
the four genes had the highest bfr are multidrug resistance-associated protein  <dig>  tumor suppressor candidate  <dig> homolog , uncharacterized protein loc <dig>  and dnaj subfamily a member  <dig>  the additional  <dig> genes with the largest bfr of over  <dig> are myoglobin, cytosolic phospholipase a <dig>  metalloreductase steap <dig>  suppression of tumorigenicity  <dig> protein, protein g7c-like, beta-tubulin, purine nucleoside phosphorylase 4b, ef-hand domain-containing protein d <dig>  serum amyloid p-component precursor, plasminogen activator inhibitor  <dig> precursor, si:dkey-269d <dig> , transmembrane protein c9orf <dig>  mhc class ii beta, bone morphogenetic protein 1a, interferon-induced very large gtpase 1-like, aldehyde dehydrogenase family  <dig> member a1-a, cc chemokine scya <dig>  rna-binding protein 8a and 40s ribosomal protein s3a .
a list of the genes with various bfr values were listed in additional file 3: table s <dig>
genes with large bfr caused by genetic segregation
as rna-seq data is analyzed in terms of rpkm at the rna level, the allele ratios obtained by rna-seq are compounded by two factors: the genotype allele frequencies at the dna level and the relative expression levels of the two alleles at the rna level.
for instance, the two alleles may have very different genotype allele frequencies in the two bulked samples, and in these cases, even if the expression is not regulated at the transcriptional level, the final ratio of the two alleles between the bulked samples are expected to be different.
however, if one of the two alleles is differentially regulated, the final allele ratio at the rna level would be different from the allele ratio at the dna level.
in order to differentiate snps with large bfr caused by genotype allele frequency difference from those with large bfr caused by allele-specific expression, the ratio of the two alleles was analyzed with combined bulk of both resistance and susceptible bulks.
theoretically, now the combined bulk should include alleles at expected segregation ratios without any connection with the phenotypes.
if the large bfr is caused by different allele frequencies at the dna level in the two bulks but not by allele-specific expression, the ratio of the two alleles in the combined bulk  should be relatively small and predictable with the family structure.
for instance, at an aa x ag snp site in a single family, the progenies should have a 3a:1g allele frequency at the genomic level.
however, the situation is more complex if more than one family is used.
nevertheless, the largest allele ratio at the dna level can still be predicted.
for instance, if two families are used as in this study, at an snp site, the largest possible allele ratio at the dna level is 7: <dig>  i.e., aa x ag in one family, and aa x aa in the second family, where the largest possible allele ratio at the dna level is 7a:1g.
any other combinations would result in a smaller combined allele ratio at the dna level .
therefore, we differentiated the large bfr caused by genetic segregation and those caused by allele-specific expression based on the combined allele ratio of the resistant and susceptible bulks.
while those large bfrs with very large combined allele ratios are likely caused, at least in part, by allele-specific expression, and those large bfrs with small combined allele ratios are likely caused by genetic segregation, the bfrs in the transitional zone could be caused by both genetic segregation and allele-specific expression .
of the  <dig> genes with large bfr ,  <dig> had a combined bulk allele ratio of  <dig> or less.
the vast majority of these had a combined bulk allele ratio of 1– <dig>  suggesting that the large bfr of these genes were not caused by allele-specific expression, and likely caused by genetic segregation.
genes with large bfr caused by allele-specific expression
as shown in figure  <dig>  a large number of genes harboring significant snps had a significantly higher combined allele ratio, with  <dig> genes had an allele ratio of greater than 9: <dig> .
considering that most of the genes with bfr > <dig> had a combined allele ratio of 1– <dig> , many of the genes with combined allele ratio of greater than  <dig> could be caused by allele-specific expression.
on the cautious side, even in the extreme cases of 7: <dig> allele ratios at the dna levels, twice the largest possible allele ratio at the dna level should be 14: <dig>  as shown in figure  <dig>  after the fisher’s exact test on the significant different level between the two alleles on the genes with allele ratio ≥  <dig>   <dig> genes had a combined allele ratio of greater than  <dig> with fdr p-value smaller than  <dig> , indicating that these large allele ratios are caused, at least in part, by allele-specific expression.
a list of these allele-specific expressed genes was provided in additional file 4: table s <dig>  of the  <dig> genes with high combined allele ratios,  <dig> genes were with bfr higher than  <dig> and allele ratio higher than  <dig>  they are plasminogen activator inhibitor  <dig>  interferon-induced very large gtpase 1-like, uncharacterized protein loc <dig> and cc chemokine scya <dig>  apparently, these large combined allele ratios were caused both by genetic segregation and by allele-specific expression.
location of genes with high bulk frequency ratio
to determine the genomic location of snps with high levels of bfr, genes containing snps with high bfr  were initially used as query for blast searches against the draft catfish genome sequence scaffolds in relation to the linkage map.
a total of  <dig> genes were identified to contain significant snps with high bfr .
blastn searches were conducted to determine the locations of the  <dig> genes on the scaffolds of the catfish genome draft sequence , and they were found to be within  <dig> genomic sequence scaffolds.
of the  <dig> scaffolds,  <dig> can be located in linkage groups, and the remaining cannot be located on linkage groups because no markers from these scaffolds were mapped.
in order to further analyze the linkage disequilibrium  and to reduce the possibility of false positives, we placed additional stringent parameters including snps with bfr ≥ <dig>  along with combined bulk allele ratio < <dig>  and the fdr adjusted p ≤  <dig>  for significant snps.
as shown in figure  <dig>  eleven lg potentially harbor genes with significant snps.
however, with only two families used in the study, if a gene is truly involved in resistance, long stretches of genomic segments are expected to be in ld because of genetic linkage.
therefore, the locations and distributions of genes containing significant snps with the following characteristics were further analyzed: 1) at least  <dig> genes were involved in the lg with significant snps; 2) or if the genes with significant snps were fewer than  <dig>  at least one snp has a bfr equal or greater than  <dig>  using this set of criteria, eight lg appeared to harbor qtls involved in esc disease resistance.
these lg were lg  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>  and  <dig>  of these eight lg, six contained genes with snps having bfr >  <dig> .
of the eight lgs, lg <dig>   <dig>  and  <dig> had the largest numbers of genes with significant snps.
detailed distributions of genes within lg were determined by locating the genes on the scaffolds along the linkage group.
as shown in figure  <dig>  in lg <dig>  snps within the von willebrand factor a domain-containing protein 7-like  gene had a very high bfr , and the nearby genes had a much lower bfr, suggesting that the resistant gene were near the g7c-like gene.
in lg  <dig>  the snp with the highest bfr was found in the middle of lg <dig> within the acidic chitinase-like gene, and bfr values in surrounding genes were gradually lower along both sides of the lg, suggesting tight linkage of the resistance gene near the acidic chitinase gene .
in lg <dig>  one snp was found to have extremely high bfr  within the dnaj subfamily a member  <dig> gene .
bfr values of snps within genes left of the dnaj subfamily a member  <dig> gene dropped sharply, but no genes with high bfr values were found on the right side of the dnaj subfamily a member  <dig> gene .
parental origin of highly expressed alleles
as blue catfish is generally more resistant to esc disease than channel catfish, an interesting question to ask is which allele was preferentially expressed after esc infection within the backcross progenies used in this study.
we therefore attempted to analyze the parental origin of the alleles for the genes with high combined allele ratios.
as shown in figure  <dig>  a total of  <dig> genes harboring snps with combined allele ratio of  <dig> or greater were identified.
of these, the parental origins could be determined for  <dig> genes with existing genome information, while the parental origin of the remaining  <dig> could not be determined.
of the  <dig> genes,  <dig> were of channel catfish origin and  <dig> were of blue catfish origin.
of the  <dig> genes preferentially expressed with channel catfish alleles, six were expressed high in resistant fish and five were expressed highly in susceptible fish.
similarly, of the  <dig> genes preferentially expressed with the blue catfish alleles, four were expressed highly in resistant fish while  <dig> were expressed highly in susceptible fish .
discussion
in this study, we conducted bsr-seq  <cit>  by combining the ngs-based rna-seq  <cit>  with bulk segregant analysis  <cit>  for the analysis of genes involved in disease response and disease resistance against enteric septicemia of catfish .
such a simple combination of rna-seq and bsa analysis allowed identification of differentially expressed genes, localization of disease resistance-related genes in linkage groups by mapping snps on the whole genome, and analysis of allele-specific expression.
bsr-seq carried the full capability of rna-seq that allowed identification of differentially expressed genes.
comparison of expression in resistant fish pool and susceptible fish pool with the control allowed the identification of differentially expressed genes after infection.
a total of  <dig>  and  <dig> genes were identified to be differentially expressed after esc infection in susceptible fish and resistant fish, respectively.
in the susceptible fish, many of the up-regulated genes represent the acute phase response protein genes, as previously reported  <cit> .
apparently, microarray studies were limited to the gene probes existing on the array, while rna-seq analysis has the ability to detect all induced genes, depending on their expression levels.
clearly, the greater numbers of genes identified from this study after infection indicated that rna-seq is more sensitive than the microarray analysis.
although rna-seq analysis was also previously conducted after esc infection  <cit> , the tissues were different in these studies.
in the work of li et al.
<cit> , intestine tissue was used, while liver was used in this study.
nevertheless, many of the differentially expressed genes identified here in the liver were among the differentially expressed genes in the intestine, as well as those identified in the microarray studies.
for instance, the acute phase response  genes such as cc chemokines, toll-like receptors, complement component proteins, catechol-o-methyltransferase domain containing  <dig>  apolipoprotein proteins, fibrinogens, angiotensinogen, mhc class i and ii, tumor necrosis factors were all found to be up-regulated, as found previously  <cit> .
the use of phenotypic extremes, resistant and susceptible fish, allowed comparison of expression patterns of genes involved in immune responses, although the time point was quite different.
for instance, a number of immune-related genes were found to be expressed higher in resistant fish than in susceptible fish including apolipoprotein a iv, apolipoprotein ab, apolipoprotein eb, apolipoprotein bb, apolipoprotein b <dig>  apolipoprotein m and complement component 1q , complement component 1 s , complement component  <dig>  complement component 3a, fibrinogen alpha, fibrinogen gamma, mhc class i, and mhc ii.
these genes expressed at higher levels in resistant fish could indicate their importance in the related disease resistance.
apolipoproteins have been shown to be important for disease resistance in mice  <cit>  and chickens  <cit> .
apparently, a much larger number of differentially expressed genes were identified in susceptible fish than in resistant fish as compared with the control.
we believe that much of this difference was caused by the different sampling time  between these two groups.
as heavy mortalities occurred 4–6 days after infection, the vast majority of differentially expressed genes in susceptible fish represented the massive responses of the host against the infection including the acute phase proteins and genes involved in inflammation and immune responses  <cit> .
in contrast, the resistant fish samples were collected two weeks after infection from the survivors.
as such, these resistant fish were either “resistant” or not infected, or may have cleared the bacteria from their system.
the massive host responses to infection may have been over by the time of two weeks after infection.
in addition to its ability to identify differentially expressed genes, bsr-seq carried the ability to identify significant snps between the pools of very strong phenotype contrast.
localization of the significant snps along the chromosomes should allow identification of positional candidate genes responsible for the trait.
in this study, use of f2-generation backcross progenies, when coupled to pooling of samples from phenotypic extremes, allowed analysis of significant snps between the resistant and susceptible fish.
we initially identified a set of significant snps between the two bulks using fisher’s exact test.
it’s the first step to remove the snps not relevant to the target trait in this study, by setting the fdr adjusted p-value ≤ <dig>  as the cutoff to identify the significant snps and reduce the computation intensive in the downstream analysis.
and the bfr of the snps with fdr adjusted p-value >  <dig>  were from  <dig> to  <dig> , which indicated that non-significant snps won’t affect the assignment of gene-level bfr.
although fisher’s exact test is technically suitable for the binary snp markers, it is compromised by the expression levels.
because the “allele frequencies” were called at the rna level, highly expressed genes had much lower p-values .
therefore, p-values were used only as the initial step for the identification of positional candidate snps linked with resistance qtls.
one gene can harbor many snps, but not all of them are relevant to disease resistant.
the low bfr of a snp means the allele frequency of that snp is similar between the susceptible group and resistant group, and such snps are irrelevant snps in relation to disease resistance.
as to the fact that different snps within a single gene can have different bfr, there may be several explanations including: 1) although the number of sequences from each allele at an snp site should be proportional to the “allele frequency”, that may not be the case practically, simply because the sequencing depth is limited; 2) we are dealing with pooled samples from two families that were derived from interspecific hybrids.
therefore, there are snp sites that are polymorphic only within channel catfish, only within blue catfish, polymorphic in both channel catfish and blue catfish, or not polymorphic within channel catfish or blue catfish, but are polymorphic between the two species.
thus allelic ratios at different snp sites are expected to be different.
given such complexities, it is reasonable to use snps with the highest bfr within each gene.
in contrast, either use of median or mean bfr or use of all the snps for each gene may bring irrelevant snps into the consideration for the analysis of snp localization, which can lead to the underestimate of the bfr in the candidate region or even miss the candidate region due to the incorrect decay pattern of ld.
in addition, it is worth to mention that we used pooled sample for the bsr-seq, which could induce the difficulty to the assessment of the within-group variance, however, this is the innate limitation of the bulked segregant analysis.
currently, there is no optimal resolution to avoid this limitation, but some studies claimed that this flaw will not cause a serious bias in the pooled sample rna-seq analysis  <cit> .
bulk frequency ratio  was previously used as an effective parameter for genetic analysis in bsa or bsr-seq  <cit> .
however, in those cases, genotypes were determined using dna.
here in this study, the “allele frequency” was calculated from the mapped reads of rna samples, and thus the calculated bfr could be compromised by allele-specific expression.
in order to identify positional candidates for resistance using transcriptome datasets generated from bsr-seq, we need to differentiate allele frequencies caused by genetic segregation and those caused by allele-specific induction/suppression: significant snps with large bfr and small combined allele ratio are likely to be caused by genetic segregation; significant snps with small bfr but large combined allele ratio are likely caused by allele-specific expression, while significant snps with large bfr and large combined allele ratios may have been caused by both genetic segregation and allele-specific expression.
for instance, if the allele ratio at an a/g snp site is  <dig> to  <dig> in resistance pool, and  <dig> to  <dig> in susceptible pool, the bfr should be  <dig>  when the two bulks were combined, now the allele ratio of a/g is 11:11 =  <dig>  this snp, with a high bfr and a low combined allele ratio, should be a snp with allele frequency difference between the bulks caused by genetic segregation.
in contrast, when one of the two alleles is differentially up-regulated, the combined allele ratio will stay large.
for instance, at an a/g snp site, if a is significantly up-regulated in the resistant fish, say 100a:5g, and in susceptible fish, a and g are roughly expressed equally, both at low levels, say 5a:5g.
in this case, the bfr = / =  <dig> ; when the two bulks are combined, the combined allele ratio would be  <dig> .
apparently in this case, the large allele ratio is caused by allele-specific expression.
in this study, each bulk was made of  <dig> fish with  <dig> fish from one family and  <dig> fish from a second family.
as the exact allele ratio at each snp site is unknown in the two families, we made several assumptions for the analysis of snps due to genetic segregation and those due to allele-specific expression.
at an a/g snp site, the parent in one family could be aa x gg, aa x ag, or ag x ag.
in these cases, the largest allele ratio can be 3: <dig>  at the dna level.
when two families were used, as in this study, the largest allele ratio at the dna level could be 7: <dig>  i.e., aa x ag in one family, and homozygous aa x aa in another family.
any snps with significantly larger combined allele ratio than 7: <dig> would suggest allele-specific expression.
we identified snps with combined allele ratio of greater than  <dig>  as being allele-specifically expressed.
apparently, many snps fell between the two possibilities, and to obtain reasonable assessment of those caused by genetic segregation and those caused by allele-specific expression, we neglected those ones between the two possibilities .
true linked snps are characteristic in that many significant snps can be observed in nearby genomic locations because of genetic linkage.
in catfish, the whole genome sequence assembly is still under way.
we therefore, mapped the significant snps to scaffolds and then examined the patterns of the snp distribution.
within a qtl region, statistical significance should be the highest with the gene underlining the performance trait, and ld should decay gradually on both sides of the chromosome surrounding the gene  <cit> .
in our study, quite many significant snps were located on unmapped scaffolds, but many were also mapped to genetic linkage groups including lg <dig>  lg <dig>  and lg <dig> that included at least  <dig> genes with high bfr  and low combined allele ratio .
as shown in figures  <dig>   <dig>  and  <dig>  the ld appeared to be decaying around the most significant snps, suggesting that these genomic regions indeed contain resistance-related genes.
for instance, in lg <dig>  the gene containing the most significant snp was protein g7c-like gene located at the 24 mb position, and the bfr values on both sides of this gene were significant, but lower than bfr within the g7c-like gene .
similarly, the gene containing the highest bfr was acidic chitinase-like gene that located in the middle of a 23 mb dna in lg <dig>  and the bfrs were lower on both sides along this lg .
in lg <dig>  the snps with highest bfr was located within the dnaj subfamily a member  <dig> gene close to the end of the 12 mb dna in lg <dig>  and the bfr on both sides were lower .
it was unknown if the detected genes with the highest bfr were themselves involved in disease resistance.
this was because some linked genes with even greater bfr were not yet mapped to the linkage group, staying as isolated scaffolds, and therefore cannot be viewed under the same “manhattan plot”, or because the expression level of the real disease resistance gene was so low that it was not detected under the bsr-seq analysis.
one logical thought is that if one gene is truly involved in disease resistance, it should be correlated with positional candidate genes as well as expression candidate genes.
in other words, it should be differentially expressed either between resistant and susceptible fish or differentially expressed after infection, and located at a genomic location the resistance phenotype is mapped.
a cross examination of differentially expressed genes and genes with high bfr but low combined allele ratio allowed identification of  <dig> genes .
only four of these  <dig> genes have been mapped to linkage groups, whereas the remaining  <dig> were just mapped to isolated scaffolds.
further mapping studies are required to determine how many qtls are involved.
these genes included some important immune function related genes such as the nlr, mhc-related genes, and mannan-binding lectin serine peptidase.
regardless of their direct involvement in the disease resistance, the importance of these genes in resistance and immune responses to esc should not be neglected.
additional research is warranted to determine if these genes are responsible for the resistance against esc.
a total of  <dig> genes were identified as genes with allele-specific expression .
one obvious question is what causes allele-specific expression .
two hypotheses were previously proposed to account for ase  <cit> .
in the first hypothesis, mutations in cis-acting dna elements can cause differences in binding of trans-acting factors, especially when such mutations are located in the promoter or enhancer regions.
inversely, mutations in the trans-acting factors would also cause their differences in binding to their target sequence.
in the second hypothesis, epigenetic factors such as differential methylation of the two alleles can cause differences in expression levels.
it has long been known that mutations in non-coding regions which affect gene expression can cause human genetic disease  <cit> .
a differentially expressed gene exhibits cis-acting variation when the differential expression is caused by factors linked to the differentially expressed alleles, such as differences in promoter sequences or chromatin state.
the list of examples in which cis-acting regulatory variation plays a key role in phenotypic variation are increasing  <cit> .
in vitro experiments prove that variants in gene promoter regions effect rates of transcription and these variants may also lead to a significant proportion of differential allelic expression  <cit> .
in addition, expressed genes contain trans-acting variation when the differential expression is caused by factors unlinked to the differentially expressed alleles, such as differences caused by genetic background and regulatory networks.
the importance of dna methylation as a driving factor in allele specific expression has been claimed and proved by a number of studies  <cit> .
in these studies, a direct correlation between allele specific expression and methylation was observed.
clearly, the different epigenetic state of each haploid genome is a major factor in the expression of the two alleles.
although x-chromosome inactivation and silencing are usually considered to be mainly related to epigenetic effect  <cit> , some studies also suggested that the change of gene regulation caused by epigenetic modification of sequence variation might be a common pathogenic mechanism in mammals  <cit> .
one hypothetical role for epigenetics is genetic imprinting leading to mono-allelic expression  <cit> .
both in vivo and in vitro experiments indicate that allele-specific differences in the rate of transcription are common existed, if not all genes are likely to show differential allelic expression in different individuals  <cit> .
however, the role of ase genes in complex traits is still not clear.
most ase studies have been conducted in humans, and no studies have been conducted in fish.
in addition to the above discussed possibilities, it is noteworthy that we used an interspecific hybrid system in this study.
in the channel catfish backcrossed progenies, overall 50% of chromosomes are “homozygous” from channel catfish while 50% chromosomes are heterozygous with one chromosome being from channel catfish and the other from blue catfish.
if a trans-acting factor is transcribed from channel catfish genes, it would bind the cis-acting elements from channel catfish with greater affinity, causing allelic expression.
this could be another explanation for the observed ase.
the significance of the observed allele-specific expression in relation to phenotype is unknown at present.
to date, the majority of expression analysis focus on the total amount of the transcripts.
however, emerging evidence underlies the importance of understanding the allele specific transcript in cases of disease  <cit> .
in our studies here, among the  <dig> allele-specifically expressed genes, parental origin of alleles can be determine for only  <dig> genes.
among these  <dig> genes, the channel catfish allele was expressed higher in  <dig> genes , and blue catfish allele was expressed higher in  <dig> genes .
there was no correlation of resistance with a specific parent.
however, it is possible that certain combinations of alleles would warrant resistance and certain combinations of alleles would lead to susceptibility.
this clearly warrants future studies.
CONCLUSIONS
in this study, we demonstrated the application of bsr-seq to study disease resistance by combining rna-seq with bulk segregant analysis.
it has the full capacity for the identification of differentially expressed genes, the capacity to identify significant snps between phenotypic bulks, the capacity to potentially identify the positional candidate genes, and the ability to identify allelic expressed genes.
among many differentially expressed genes,  <dig> genes were also among the genes that had high bfr and low combined allele ratio, suggesting that these genes could be potentially involved in disease resistance.

 performance is related to the subtle variation in gene expression and this relationship differs among individuals. in well-defined families, the first level of variation comes from genetic segregation and recombination of chromosomes. as a result of segregation and chromosomal recombination, each individual has different genetic makeup. upon a given genetic background, genetic potential carried on dna can only be realized when the genes are expressed. at the whole genome level, expression of each gene is affected by its genetic regulatory element as well as trans-acting factors including the impact of environment. a composite of genes, transcriptional regulation, post-transcriptional modification and regulation, translational regulation and post-translational modification and regulation, along with environmental impact and genotype-environment interactions eventually determines the phenotypes of individuals. when considered at the whole genome level, expression of tens of thousands of genes and combination of these genes make the variation of performance traits extremely complex with huge variability. the task of modern agricultural genomics is to gain understanding of such variations and their relationship in determination of production and performance traits. traditionally, genetic and molecular biological studies are conducted to dissect these variables at different levels. for instance, the effect of alleles can be dissected through genetic and qtl mapping analysis. gene expression can be analyzed using high throughput methodologies such as microarrays and rna-seq analysis. various epigenetic regulations have also been studied to understand the differences in gene expression with similar genetic background. such analyses have been very powerful in determination of genetic and epigenetic factors affecting performance and production traits. however, performance and production traits are often highly complex and the outcome of agricultural operations is affected by variations at all levels. for example, genetic background is very important because disease resistance genes allow the organism to survive the serious infections. in most cases, disease resistance genes have been studied through genetic linkage and qtl analyses that allow the identification of genomic regions containing disease resistance genes to be identified. even with the most powerful molecular approaches, analysis of complex traits such as disease resistance can be extremely challenging. in michelmore et al. developed a method called bulked segregant analysis to study disease resistance in plants. the basic idea of bsa was that phenotypic extremes should have drastic differences in genotypes. when samples are selected from phenotypic extremes, say the best and the worst performers, and their genotypes are analyzed in bulk, a correlation of genotypes with phenotypes can be attained. in other word, the variation among individuals may be quite subtle and difficult to detect; however, the pooled samples of the phenotypic extremes should pose a strong contrast in their genotypes at the genomic location linked to the trait. bsa has been used in numerous studies to associate phenotypes with related genomic locations. bsa has been evolving along with various types of molecular markers including restriction fragment length polymorphisms, random amplified polymorphic dnas, simple sequence repeats, amplified fragment length polymorphisms and single nucleotide polymorphisms. with the development of next-generation sequencing technologies, bsa was first enhanced by the application of sequence-based markers such as restriction-site associated dna markers and whole genome sequencing. in recent years, the application of rna-seq has allowed rapid and comprehensive understanding of transcriptome level of variations. coupling of bsa with rna-seq has led to the development of bulked segregant rna-seq, and it has been successfully applied in plants, but not yet demonstrated in animals. apparently, bsr-seq possesses the advantages of both bsa and rna-seq, with the high throughput for deep coverage of the transcriptome as well as the strong ability to detect genetic differences underlining the traits. such a technique is best suited to organisms with high fecundity such as many species of fish. catfish is the major aquaculture species in the united states, accounting for over 60% of all us aquaculture production. the two major cultivated catfish species are channel catfish and blue catfish. an inter-specific hybrid has been popular for aquaculture because of strong heterosis. not only is the interspecific hybrid is popular for aquaculture, it is also a superior system to study disease resistance because of their strong phenotypic difference. blue catfish is extremely resistant against esc disease while channel catfish is relatively susceptible. genetic linkage analysis of f generation of the interspecific hybrids would allow identification of disease resistance/susceptibility genes. in this study, we take advantage of bsr-seq for the analysis of disease resistant genes using the f generation backcross progenies of the interspecific hybrids. here we demonstrate that bsr-seq is capable of 1) revealing differentially expressed genes; 2) revealing positional candidates containing genes related to disease resistance after mapping snps on the whole genome; and 3) revealing allele-specific expression after bacterial infection.  sequence assembly and analysis rna-seq was conducted using illumina sequencing with three pooled samples of resistant fish, susceptible fish, and control fish. each pooled sample contained equal amount of rna collected from individuals. each pooled sample was barcoded such that reads from each of the three samples can be traced and analyzed separately. a total of over million reads were generated with million reads from the resistant fish, million from the susceptible fish, and million from the control fish. after quality trimming, a total of million reads were carried forward for analysis, with an average read length of  bp. raw read data are archived at the ncbi sequence read archive under accession srp the reads after quality trimming were assembled de novo into non-redundant contigs with a n contig length of  bp and an average contig length of 825 bp. of the assembled contigs, over had a length of over  bp. the assembled contigs were analyzed by blast searches to determine the gene identities of the contigs. of the contigs, had hits to unigenes, of which were known genes and were hypothetical genes. de novo assembly of the catfish liver transcriptome with infection of esc generated by illumina sequencing and assembled with trinity contigs were used as queries for blast searches against zebrafish and nr protein databases. the cut-off value for gene identity was set at e-value ≤ 1e- differentially expressed genes after infection differentially expressed genes after infection in resistant fish and susceptible fish were determined by comparing their expression levels in rpkm with that of the control group. as summarized in table  a total of genes were differentially expressed in the resistant group as compared with the control group. a significant fraction of the differentially expressed genes in resistant fish were highly regulated, with genes being up- or down-regulated 10-fold or more, and additional genes up- or down-regulated 5– fold. edwardsiella ictaluri although a relatively small number of genes were differentially expressed in resistant fish, a large number of genes were differentially expressed in susceptible fish, with a total of genes being differentially expressed. not only the number of differentially expressed genes was drastically more in susceptible fish, the number of highly regulated genes was also much greater in susceptible fish, with genes being up- or down-regulated folds or more, and additional genes were up- or down-regulated 5– folds. comparison of gene expression in resistant fish and susceptible fish after infection although a small number of genes were differentially expressed in resistant fish when being compared with the control fish, a large number of genes exhibited differential expression in resistant fish when being compared with the susceptible fish. a total of genes were differentially expressed in resistant and susceptible fish, with genes expressed significantly higher in resistant fish than in susceptible fish and genes expressed significantly lower in resistant fish than in susceptible fish. of the genes expressed significantly higher in resistant fish, genes were expressed over times more than in susceptible fish; were expressed 50– times more than in susceptible fish; were expressed 10– times more than in susceptible fish. of the genes expressed significantly lower in resistant fish, genes were expressed over times less in resistant fish than in susceptible fish; were expressed 50– times less in resistant fish than in susceptible fish; were expressed 10– times less in resistant fish than in susceptible fish. edwardsiella ictaluri these genes are listed in additional file 1: table s identification of snps and significant snps snps were identified by alignment of short reads to the reference assembly of the rna-seq. in order to be qualified for snps, at least reads were required for each group and a total of minor allele reads count must be greater than among three groups. using software popoolation a total of snps were identified. in order to be sure the identified snps were properly identified, a second software package, varscan was also used. a total of snps were identified. the difference was caused by differences in cutoff of quality scores of sequence reads by the two programs. although the total numbers of snps identified by the two software differed by approximately 6%, the vast majority of snps identified by the two software were identical. these snps were located within contigs. in order to determine snps with significant difference in allele frequencies between the resistant and susceptible fish, fisher’s exact test was initially performed. as shown in table  snps were identified as significant snps. these significant snps were located in contigs. of the contigs, had significant hits to known genes, and the remaining probably represented contigs assembled from 5′- and 3′-untranslated regions or from long non-coding rnas. the contigs with hits to known genes represented unique genes harboring snps. bulk frequency ratios although significant snps identified through rna-seq analysis reflect the final ratios of different alleles at the rna level in the two bulked samples, statistical analysis for significant snps using fisher’s exact test was only the first step for screening snps that may be significantly associated with the trait. in order to compare the snp allele frequencies more directly, bulk frequency ratios were generated from the rna-seq data between the two bulks, the resistant fish and the susceptible fish. the bfr of genes were determined by the maximum bfr of the significant snps generated from the fisher’s exact test located in this gene. as shown in figure  large proportion of genes containing the significant snps had a bfr of at least a total of genes had a bfr equal or greater than among these genes, genes had a bfr of 4–16; genes had a bfr over 16; and genes had a bfr over. the four genes had the highest bfr are multidrug resistance-associated protein tumor suppressor candidate homolog, uncharacterized protein loc and dnaj subfamily a member the additional genes with the largest bfr of over are myoglobin, cytosolic phospholipase a metalloreductase steap suppression of tumorigenicity protein, protein g7c-like, beta-tubulin, purine nucleoside phosphorylase 4b, ef-hand domain-containing protein d serum amyloid p-component precursor, plasminogen activator inhibitor precursor, si:dkey-269d, transmembrane protein c9orf mhc class ii beta, bone morphogenetic protein 1a, interferon-induced very large gtpase 1-like, aldehyde dehydrogenase family member a1-a, cc chemokine scya rna-binding protein 8a and 40s ribosomal protein s3a. a list of the genes with various bfr values were listed in additional file 3: table s genes with large bfr caused by genetic segregation as rna-seq data is analyzed in terms of rpkm at the rna level, the allele ratios obtained by rna-seq are compounded by two factors: the genotype allele frequencies at the dna level and the relative expression levels of the two alleles at the rna level. for instance, the two alleles may have very different genotype allele frequencies in the two bulked samples, and in these cases, even if the expression is not regulated at the transcriptional level, the final ratio of the two alleles between the bulked samples are expected to be different. however, if one of the two alleles is differentially regulated, the final allele ratio at the rna level would be different from the allele ratio at the dna level. in order to differentiate snps with large bfr caused by genotype allele frequency difference from those with large bfr caused by allele-specific expression, the ratio of the two alleles was analyzed with combined bulk of both resistance and susceptible bulks. theoretically, now the combined bulk should include alleles at expected segregation ratios without any connection with the phenotypes. if the large bfr is caused by different allele frequencies at the dna level in the two bulks but not by allele-specific expression, the ratio of the two alleles in the combined bulk should be relatively small and predictable with the family structure. for instance, at an aa x ag snp site in a single family, the progenies should have a 3a:1g allele frequency at the genomic level. however, the situation is more complex if more than one family is used. nevertheless, the largest allele ratio at the dna level can still be predicted. for instance, if two families are used as in this study, at an snp site, the largest possible allele ratio at the dna level is 7: i.e., aa x ag in one family, and aa x aa in the second family, where the largest possible allele ratio at the dna level is 7a:1g. any other combinations would result in a smaller combined allele ratio at the dna level. therefore, we differentiated the large bfr caused by genetic segregation and those caused by allele-specific expression based on the combined allele ratio of the resistant and susceptible bulks. while those large bfrs with very large combined allele ratios are likely caused, at least in part, by allele-specific expression, and those large bfrs with small combined allele ratios are likely caused by genetic segregation, the bfrs in the transitional zone could be caused by both genetic segregation and allele-specific expression. of the genes with large bfr, had a combined bulk allele ratio of or less. the vast majority of these had a combined bulk allele ratio of 1– suggesting that the large bfr of these genes were not caused by allele-specific expression, and likely caused by genetic segregation. genes with large bfr caused by allele-specific expression as shown in figure  a large number of genes harboring significant snps had a significantly higher combined allele ratio, with genes had an allele ratio of greater than 9:. considering that most of the genes with bfr > had a combined allele ratio of 1–, many of the genes with combined allele ratio of greater than could be caused by allele-specific expression. on the cautious side, even in the extreme cases of 7: allele ratios at the dna levels, twice the largest possible allele ratio at the dna level should be 14: as shown in figure  after the fisher’s exact test on the significant different level between the two alleles on the genes with allele ratio ≥  genes had a combined allele ratio of greater than with fdr p-value smaller than, indicating that these large allele ratios are caused, at least in part, by allele-specific expression. a list of these allele-specific expressed genes was provided in additional file 4: table s of the genes with high combined allele ratios, genes were with bfr higher than and allele ratio higher than they are plasminogen activator inhibitor interferon-induced very large gtpase 1-like, uncharacterized protein loc and cc chemokine scya apparently, these large combined allele ratios were caused both by genetic segregation and by allele-specific expression. location of genes with high bulk frequency ratio to determine the genomic location of snps with high levels of bfr, genes containing snps with high bfr were initially used as query for blast searches against the draft catfish genome sequence scaffolds in relation to the linkage map. a total of genes were identified to contain significant snps with high bfr. blastn searches were conducted to determine the locations of the genes on the scaffolds of the catfish genome draft sequence, and they were found to be within genomic sequence scaffolds. of the scaffolds, can be located in linkage groups, and the remaining cannot be located on linkage groups because no markers from these scaffolds were mapped. in order to further analyze the linkage disequilibrium and to reduce the possibility of false positives, we placed additional stringent parameters including snps with bfr ≥ along with combined bulk allele ratio < and the fdr adjusted p ≤  for significant snps. as shown in figure  eleven lg potentially harbor genes with significant snps. however, with only two families used in the study, if a gene is truly involved in resistance, long stretches of genomic segments are expected to be in ld because of genetic linkage. therefore, the locations and distributions of genes containing significant snps with the following characteristics were further analyzed: 1) at least genes were involved in the lg with significant snps; 2) or if the genes with significant snps were fewer than at least one snp has a bfr equal or greater than using this set of criteria, eight lg appeared to harbor qtls involved in esc disease resistance. these lg were lg and of these eight lg, six contained genes with snps having bfr > . of the eight lgs, lg and had the largest numbers of genes with significant snps. detailed distributions of genes within lg were determined by locating the genes on the scaffolds along the linkage group. as shown in figure  in lg snps within the von willebrand factor a domain-containing protein 7-like gene had a very high bfr, and the nearby genes had a much lower bfr, suggesting that the resistant gene were near the g7c-like gene. in lg the snp with the highest bfr was found in the middle of lg within the acidic chitinase-like gene, and bfr values in surrounding genes were gradually lower along both sides of the lg, suggesting tight linkage of the resistance gene near the acidic chitinase gene. in lg one snp was found to have extremely high bfr within the dnaj subfamily a member gene. bfr values of snps within genes left of the dnaj subfamily a member gene dropped sharply, but no genes with high bfr values were found on the right side of the dnaj subfamily a member gene. parental origin of highly expressed alleles as blue catfish is generally more resistant to esc disease than channel catfish, an interesting question to ask is which allele was preferentially expressed after esc infection within the backcross progenies used in this study. we therefore attempted to analyze the parental origin of the alleles for the genes with high combined allele ratios. as shown in figure  a total of genes harboring snps with combined allele ratio of or greater were identified. of these, the parental origins could be determined for genes with existing genome information, while the parental origin of the remaining could not be determined. of the genes, were of channel catfish origin and were of blue catfish origin. of the genes preferentially expressed with channel catfish alleles, six were expressed high in resistant fish and five were expressed highly in susceptible fish. similarly, of the genes preferentially expressed with the blue catfish alleles, four were expressed highly in resistant fish while were expressed highly in susceptible fish. discussion in this study, we conducted bsr-seq by combining the ngs-based rna-seq with bulk segregant analysis for the analysis of genes involved in disease response and disease resistance against enteric septicemia of catfish. such a simple combination of rna-seq and bsa analysis allowed identification of differentially expressed genes, localization of disease resistance-related genes in linkage groups by mapping snps on the whole genome, and analysis of allele-specific expression. bsr-seq carried the full capability of rna-seq that allowed identification of differentially expressed genes. comparison of expression in resistant fish pool and susceptible fish pool with the control allowed the identification of differentially expressed genes after infection. a total of and genes were identified to be differentially expressed after esc infection in susceptible fish and resistant fish, respectively. in the susceptible fish, many of the up-regulated genes represent the acute phase response protein genes, as previously reported. apparently, microarray studies were limited to the gene probes existing on the array, while rna-seq analysis has the ability to detect all induced genes, depending on their expression levels. clearly, the greater numbers of genes identified from this study after infection indicated that rna-seq is more sensitive than the microarray analysis. although rna-seq analysis was also previously conducted after esc infection, the tissues were different in these studies. in the work of li et al., intestine tissue was used, while liver was used in this study. nevertheless, many of the differentially expressed genes identified here in the liver were among the differentially expressed genes in the intestine, as well as those identified in the microarray studies. for instance, the acute phase response genes such as cc chemokines, toll-like receptors, complement component proteins, catechol-o-methyltransferase domain containing apolipoprotein proteins, fibrinogens, angiotensinogen, mhc class i and ii, tumor necrosis factors were all found to be up-regulated, as found previously. the use of phenotypic extremes, resistant and susceptible fish, allowed comparison of expression patterns of genes involved in immune responses, although the time point was quite different. for instance, a number of immune-related genes were found to be expressed higher in resistant fish than in susceptible fish including apolipoprotein a iv, apolipoprotein ab, apolipoprotein eb, apolipoprotein bb, apolipoprotein b apolipoprotein m and complement component 1q, complement component 1 s, complement component complement component 3a, fibrinogen alpha, fibrinogen gamma, mhc class i, and mhc ii. these genes expressed at higher levels in resistant fish could indicate their importance in the related disease resistance. apolipoproteins have been shown to be important for disease resistance in mice and chickens. apparently, a much larger number of differentially expressed genes were identified in susceptible fish than in resistant fish as compared with the control. we believe that much of this difference was caused by the different sampling time between these two groups. as heavy mortalities occurred 4–6 days after infection, the vast majority of differentially expressed genes in susceptible fish represented the massive responses of the host against the infection including the acute phase proteins and genes involved in inflammation and immune responses. in contrast, the resistant fish samples were collected two weeks after infection from the survivors. as such, these resistant fish were either “resistant” or not infected, or may have cleared the bacteria from their system. the massive host responses to infection may have been over by the time of two weeks after infection. in addition to its ability to identify differentially expressed genes, bsr-seq carried the ability to identify significant snps between the pools of very strong phenotype contrast. localization of the significant snps along the chromosomes should allow identification of positional candidate genes responsible for the trait. in this study, use of f2-generation backcross progenies, when coupled to pooling of samples from phenotypic extremes, allowed analysis of significant snps between the resistant and susceptible fish. we initially identified a set of significant snps between the two bulks using fisher’s exact test. it’s the first step to remove the snps not relevant to the target trait in this study, by setting the fdr adjusted p-value ≤ as the cutoff to identify the significant snps and reduce the computation intensive in the downstream analysis. and the bfr of the snps with fdr adjusted p-value >  were from to, which indicated that non-significant snps won’t affect the assignment of gene-level bfr. although fisher’s exact test is technically suitable for the binary snp markers, it is compromised by the expression levels. because the “allele frequencies” were called at the rna level, highly expressed genes had much lower p-values. therefore, p-values were used only as the initial step for the identification of positional candidate snps linked with resistance qtls. one gene can harbor many snps, but not all of them are relevant to disease resistant. the low bfr of a snp means the allele frequency of that snp is similar between the susceptible group and resistant group, and such snps are irrelevant snps in relation to disease resistance. as to the fact that different snps within a single gene can have different bfr, there may be several explanations including: 1) although the number of sequences from each allele at an snp site should be proportional to the “allele frequency”, that may not be the case practically, simply because the sequencing depth is limited; 2) we are dealing with pooled samples from two families that were derived from interspecific hybrids. therefore, there are snp sites that are polymorphic only within channel catfish, only within blue catfish, polymorphic in both channel catfish and blue catfish, or not polymorphic within channel catfish or blue catfish, but are polymorphic between the two species. thus allelic ratios at different snp sites are expected to be different. given such complexities, it is reasonable to use snps with the highest bfr within each gene. in contrast, either use of median or mean bfr or use of all the snps for each gene may bring irrelevant snps into the consideration for the analysis of snp localization, which can lead to the underestimate of the bfr in the candidate region or even miss the candidate region due to the incorrect decay pattern of ld. in addition, it is worth to mention that we used pooled sample for the bsr-seq, which could induce the difficulty to the assessment of the within-group variance, however, this is the innate limitation of the bulked segregant analysis. currently, there is no optimal resolution to avoid this limitation, but some studies claimed that this flaw will not cause a serious bias in the pooled sample rna-seq analysis. bulk frequency ratio was previously used as an effective parameter for genetic analysis in bsa or bsr-seq. however, in those cases, genotypes were determined using dna. here in this study, the “allele frequency” was calculated from the mapped reads of rna samples, and thus the calculated bfr could be compromised by allele-specific expression. in order to identify positional candidates for resistance using transcriptome datasets generated from bsr-seq, we need to differentiate allele frequencies caused by genetic segregation and those caused by allele-specific induction/suppression: significant snps with large bfr and small combined allele ratio are likely to be caused by genetic segregation; significant snps with small bfr but large combined allele ratio are likely caused by allele-specific expression, while significant snps with large bfr and large combined allele ratios may have been caused by both genetic segregation and allele-specific expression. for instance, if the allele ratio at an a/g snp site is to in resistance pool, and to in susceptible pool, the bfr should be when the two bulks were combined, now the allele ratio of a/g is 11:11 =  this snp, with a high bfr and a low combined allele ratio, should be a snp with allele frequency difference between the bulks caused by genetic segregation. in contrast, when one of the two alleles is differentially up-regulated, the combined allele ratio will stay large. for instance, at an a/g snp site, if a is significantly up-regulated in the resistant fish, say 100a:5g, and in susceptible fish, a and g are roughly expressed equally, both at low levels, say 5a:5g. in this case, the bfr = / = ; when the two bulks are combined, the combined allele ratio would be. apparently in this case, the large allele ratio is caused by allele-specific expression. in this study, each bulk was made of fish with fish from one family and fish from a second family. as the exact allele ratio at each snp site is unknown in the two families, we made several assumptions for the analysis of snps due to genetic segregation and those due to allele-specific expression. at an a/g snp site, the parent in one family could be aa x gg, aa x ag, or ag x ag. in these cases, the largest allele ratio can be 3: at the dna level. when two families were used, as in this study, the largest allele ratio at the dna level could be 7: i.e., aa x ag in one family, and homozygous aa x aa in another family. any snps with significantly larger combined allele ratio than 7: would suggest allele-specific expression. we identified snps with combined allele ratio of greater than as being allele-specifically expressed. apparently, many snps fell between the two possibilities, and to obtain reasonable assessment of those caused by genetic segregation and those caused by allele-specific expression, we neglected those ones between the two possibilities. true linked snps are characteristic in that many significant snps can be observed in nearby genomic locations because of genetic linkage. in catfish, the whole genome sequence assembly is still under way. we therefore, mapped the significant snps to scaffolds and then examined the patterns of the snp distribution. within a qtl region, statistical significance should be the highest with the gene underlining the performance trait, and ld should decay gradually on both sides of the chromosome surrounding the gene. in our study, quite many significant snps were located on unmapped scaffolds, but many were also mapped to genetic linkage groups including lg lg and lg that included at least genes with high bfr and low combined allele ratio. as shown in figures  and the ld appeared to be decaying around the most significant snps, suggesting that these genomic regions indeed contain resistance-related genes. for instance, in lg the gene containing the most significant snp was protein g7c-like gene located at the 24 mb position, and the bfr values on both sides of this gene were significant, but lower than bfr within the g7c-like gene. similarly, the gene containing the highest bfr was acidic chitinase-like gene that located in the middle of a 23 mb dna in lg and the bfrs were lower on both sides along this lg. in lg the snps with highest bfr was located within the dnaj subfamily a member gene close to the end of the 12 mb dna in lg and the bfr on both sides were lower. it was unknown if the detected genes with the highest bfr were themselves involved in disease resistance. this was because some linked genes with even greater bfr were not yet mapped to the linkage group, staying as isolated scaffolds, and therefore cannot be viewed under the same “manhattan plot”, or because the expression level of the real disease resistance gene was so low that it was not detected under the bsr-seq analysis. one logical thought is that if one gene is truly involved in disease resistance, it should be correlated with positional candidate genes as well as expression candidate genes. in other words, it should be differentially expressed either between resistant and susceptible fish or differentially expressed after infection, and located at a genomic location the resistance phenotype is mapped. a cross examination of differentially expressed genes and genes with high bfr but low combined allele ratio allowed identification of genes. only four of these genes have been mapped to linkage groups, whereas the remaining were just mapped to isolated scaffolds. further mapping studies are required to determine how many qtls are involved. these genes included some important immune function related genes such as the nlr, mhc-related genes, and mannan-binding lectin serine peptidase. regardless of their direct involvement in the disease resistance, the importance of these genes in resistance and immune responses to esc should not be neglected. additional research is warranted to determine if these genes are responsible for the resistance against esc. a total of genes were identified as genes with allele-specific expression. one obvious question is what causes allele-specific expression. two hypotheses were previously proposed to account for ase. in the first hypothesis, mutations in cis-acting dna elements can cause differences in binding of trans-acting factors, especially when such mutations are located in the promoter or enhancer regions. inversely, mutations in the trans-acting factors would also cause their differences in binding to their target sequence. in the second hypothesis, epigenetic factors such as differential methylation of the two alleles can cause differences in expression levels. it has long been known that mutations in non-coding regions which affect gene expression can cause human genetic disease. a differentially expressed gene exhibits cis-acting variation when the differential expression is caused by factors linked to the differentially expressed alleles, such as differences in promoter sequences or chromatin state. the list of examples in which cis-acting regulatory variation plays a key role in phenotypic variation are increasing. in vitro experiments prove that variants in gene promoter regions effect rates of transcription and these variants may also lead to a significant proportion of differential allelic expression. in addition, expressed genes contain trans-acting variation when the differential expression is caused by factors unlinked to the differentially expressed alleles, such as differences caused by genetic background and regulatory networks. the importance of dna methylation as a driving factor in allele specific expression has been claimed and proved by a number of studies. in these studies, a direct correlation between allele specific expression and methylation was observed. clearly, the different epigenetic state of each haploid genome is a major factor in the expression of the two alleles. although x-chromosome inactivation and silencing are usually considered to be mainly related to epigenetic effect, some studies also suggested that the change of gene regulation caused by epigenetic modification of sequence variation might be a common pathogenic mechanism in mammals. one hypothetical role for epigenetics is genetic imprinting leading to mono-allelic expression. both in vivo and in vitro experiments indicate that allele-specific differences in the rate of transcription are common existed, if not all genes are likely to show differential allelic expression in different individuals. however, the role of ase genes in complex traits is still not clear. most ase studies have been conducted in humans, and no studies have been conducted in fish. in addition to the above discussed possibilities, it is noteworthy that we used an interspecific hybrid system in this study. in the channel catfish backcrossed progenies, overall 50% of chromosomes are “homozygous” from channel catfish while 50% chromosomes are heterozygous with one chromosome being from channel catfish and the other from blue catfish. if a trans-acting factor is transcribed from channel catfish genes, it would bind the cis-acting elements from channel catfish with greater affinity, causing allelic expression. this could be another explanation for the observed ase. the significance of the observed allele-specific expression in relation to phenotype is unknown at present. to date, the majority of expression analysis focus on the total amount of the transcripts. however, emerging evidence underlies the importance of understanding the allele specific transcript in cases of disease. in our studies here, among the allele-specifically expressed genes, parental origin of alleles can be determine for only genes. among these genes, the channel catfish allele was expressed higher in genes, and blue catfish allele was expressed higher in genes. there was no correlation of resistance with a specific parent. however, it is possible that certain combinations of alleles would warrant resistance and certain combinations of alleles would lead to susceptibility. this clearly warrants future studies.  in this study, we demonstrated the application of bsr-seq to study disease resistance by combining rna-seq with bulk segregant analysis. it has the full capacity for the identification of differentially expressed genes, the capacity to identify significant snps between phenotypic bulks, the capacity to potentially identify the positional candidate genes, and the ability to identify allelic expressed genes. among many differentially expressed genes, genes were also among the genes that had high bfr and low combined allele ratio, suggesting that these genes could be potentially involved in disease resistance. 
detailed analysis of these significant snps allowed differentiation of significant snps caused by genetic segregation and those caused by allele-specific expression. in addition, snps residing on unique genes were identified as significant snps between susceptible and resistant fish. coupling of rna-seq with bulked segregant analysis should allow correlation of expression patterns and associated snps with the phenotypes. in this study, we demonstrated the use of bulked segregant rna-seq for the analysis of differentially expressed genes and associated snps with disease resistance against enteric septicemia of catfish. bsr-seq is applicable to analysis of genes underlining various performance and production traits without significant investment in the development of large genotyping platforms such as snp arrays. mapping of the significant snps, along with analysis of differentially expressed genes, allowed identification of candidate genes underlining disease resistance against esc disease. a total of differentially expressed genes were found between resistant and susceptible fish. this study demonstrated the use of bsr-seq for the identification of genes involved in disease resistance against esc through expression profiling and mapping of significantly associated snps. the application of rna-seq has accelerated gene expression profiling and identification of gene-associated snps in many species. however, the integrated studies of gene expression along with snp mapping have been lacking.
 the application of rna-seq has accelerated gene expression profiling and identification of gene-associated snps in many species. however, the integrated studies of gene expression along with snp mapping have been lacking. coupling of rna-seq with bulked segregant analysis should allow correlation of expression patterns and associated snps with the phenotypes.  in this study, we demonstrated the use of bulked segregant rna-seq for the analysis of differentially expressed genes and associated snps with disease resistance against enteric septicemia of catfish. a total of differentially expressed genes were found between resistant and susceptible fish. in addition, snps residing on unique genes were identified as significant snps between susceptible and resistant fish. detailed analysis of these significant snps allowed differentiation of significant snps caused by genetic segregation and those caused by allele-specific expression. mapping of the significant snps, along with analysis of differentially expressed genes, allowed identification of candidate genes underlining disease resistance against esc disease.  this study demonstrated the use of bsr-seq for the identification of genes involved in disease resistance against esc through expression profiling and mapping of significantly associated snps. bsr-seq is applicable to analysis of genes underlining various performance and production traits without significant investment in the development of large genotyping platforms such as snp arrays. bulk segregant analysisrna-seqdisease resistancecatfishallele-specific expression 