A recent and powerful technique is to obtain transcriptomes from rare cell populations, such as single neurons in Caenorhabditis elegans, by enriching dissociated cells using fluorescent sorting. However, these cell samples often have low yields of RNA that present challenges in library preparation. This can lead to PCR duplicates, noisy gene expression for lowly expressed genes, and other issues that limit endpoint analysis. Furthermore, some common resources, such as sequence-specific kits for removing ribosomal RNA, are not optimized for nonmammalian samples. To advance library construction for such challenging samples, we compared two approaches for building RNAseq libraries from less than 10 nanograms of C. elegans RNA: SMARTSeq V4 (Takara), a widely used kit for selecting poly-adenylated transcripts; and SoLo Ovation (Tecan Genomics), a newly developed ribodepletion-based approach. For ribodepletion, we used a custom kit of 200 probes designed to match C. elegans rRNA gene sequences. We found that SoLo Ovation, in combination with our custom C. elegans probe set for rRNA depletion, detects an expanded set of noncoding RNAs, shows reduced noise in lowly expressed genes, and more accurately counts expression of long genes. The approach described here should be broadly useful for similar efforts to analyze transcriptomics when RNA is limiting.

During the life cycle of shellfish, larval development, especially metamorphosis, has a vital influence on the dynamics, distribution, and recruitment of natural populations, as well as seed breeding. Rapana venosa, a carnivorous gastropod, is an important commercial shellfish in China, and is an ecological invader in the United States, Argentina, and France. However, information about the mechanism of its early development is still limited, because research in this area has long suffered from a lack of genomic resources. In this study, 15 digital gene expression (DGE) libraries from five developmental stages of R. venosa were constructed and sequenced on the IIIumina Hi-Sequation 2500 platform. Bioinformaticsanalysis identified numerous differentially and specifically expressed genes, which revealed that genes associated with growth, nervous system, digestive system, immune system, and apoptosis participate in important developmental processes. The functional analysis of differentially expressed genes was further implemented by gene ontology, and Kyoto encyclopedia of genes and genomes enrichment. DGE profiling provided a general picture of the transcriptomic activities during the early development of R. venosa, which may provide interesting hints for further study. Our data represent the first comparative transcriptomic information available for the early development of R. venosa, which is a prerequisite for a better understanding of the physiological traits controlling development.

The Spt-Ada-Gcn5 Acetyltransferase (SAGA) complex is a transcriptional coactivator with histone acetylase and deubiquitinase activities that plays an important role in visual development and function. In Drosophila melanogaster, four SAGA subunits are required for the deubiquitination of monoubiquitinated histone H2B (ubH2B): Nonstop, Sgf11, E(y)2, and Ataxin 7. Mutations that disrupt SAGA deubiquitinase activity cause defects in neuronal connectivity in the developing Drosophila visual system. In addition, mutations in SAGA result in the human progressive visual disorder spinocerebellar ataxia type 7 (SCA7). Glial cells play a crucial role in both the neuronal connectivity defect in nonstop and sgf11 flies, and in the retinal degeneration observed in SCA7 patients. Thus, we sought to identify the gene targets of SAGA deubiquitinase activity in glia in the Drosophila larval central nervous system. To do this, we enriched glia from wild-type, nonstop, and sgf11 larval optic lobes using affinity-purification of KASH-GFP tagged nuclei, and then examined each transcriptome using RNA-seq. Our analysis showed that SAGA deubiquitinase activity is required for proper expression of 16% of actively transcribed genes in glia, especially genes involved in proteasome function, protein folding and axon guidance. We further show that the SAGA deubiquitinase-activated gene Multiplexin (Mp) is required in glia for proper photoreceptor axon targeting. Mutations in the human ortholog of Mp, COL18A1, have been identified in a family with a SCA7-like progressive visual disorder, suggesting that defects in the expression of this gene in SCA7 patients could play a role in the retinal degeneration that is unique to this ataxia.

Kinship analyses are important pillars of ecological and conservation genetic studies with potentially far-reaching implications. There is a need for power analyses that address a range of possible relationships. Nevertheless, such analyses are rarely applied, and studies that use genetic-data-based-kinship inference often ignore the influence of intrinsic population characteristics. We investigated 11 questions regarding the correct classification rate of dyads to relatedness categories (relatedness category assignments; RCA) using an individual-based model with realistic life history parameters. We investigated the effects of the number of genetic markers; marker type (microsatellite, single nucleotide polymorphism SNP, or both); minor allele frequency; typing error; mating system; and the number of overlapping generations under different demographic conditions. We found that (i) an increasing number of genetic markers increased the correct classification rate of the RCA so that up to >80% first cousins can be correctly assigned; (ii) the minimum number of genetic markers required for assignments with 80 and 95% correct classifications differed between relatedness categories, mating systems, and the number of overlapping generations; (iii) the correct classification rate was improved by adding additional relatedness categories and age and mitochondrial DNA data; and (iv) a combination of microsatellite and single-nucleotide polymorphism data increased the correct classification rate if <800 SNP loci were available. This study shows how intrinsic population characteristics, such as mating system and the number of overlapping generations, life history traits, and genetic marker characteristics, can influence the correct classification rate of an RCA study. Therefore, species-specific power analyses are essential for empirical studies.

The usefulness of genomic prediction in crop and livestock breeding programs has prompted efforts to develop new and improved genomic prediction algorithms, such as artificial neural networks and gradient tree boosting. However, the performance of these algorithms has not been compared in a systematic manner using a wide range of datasets and models. Using data of 18 traits across six plant species with different marker densities and training population sizes, we compared the performance of six linear and six non-linear algorithms. First, we found that hyperparameter selection was necessary for all non-linear algorithms and that feature selection prior to model training was critical for artificial neural networks when the markers greatly outnumbered the number of training lines. Across all species and trait combinations, no one algorithm performed best, however predictions based on a combination of results from multiple algorithms (i.e., ensemble predictions) performed consistently well. While linear and non-linear algorithms performed best for a similar number of traits, the performance of non-linear algorithms vary more between traits. Although artificial neural networks did not perform best for any trait, we identified strategies (i.e., feature selection, seeded starting weights) that boosted their performance to near the level of other algorithms. Our results highlight the importance of algorithm selection for the prediction of trait values.

Exposure to the mycotoxin aflatoxin B1 (AFB1) strongly correlates with hepatocellular carcinoma (HCC). P450 enzymes convert AFB1 into a highly reactive epoxide that forms unstable 8,9-dihydro-8-(N7-guanyl)-9-hydroxyaflatoxin B1 (AFB1-N7-Gua) DNA adducts, which convert to stable mutagenic AFB1 formamidopyrimidine (FAPY) DNA adducts. In CYP1A2-expressing budding yeast, AFB1 is a weak mutagen but a potent recombinagen. However, few genes have been identified that confer AFB1 resistance. Here, we profiled the yeast genome for AFB1 resistance. We introduced the human CYP1A2 into ∼90% of the diploid deletion library, and pooled samples from CYP1A2-expressing libraries and the original library were exposed to 50 μM AFB1 for 20 hs. By using next generation sequencing (NGS) to count molecular barcodes, we initially identified 86 genes from the CYP1A2-expressing libraries, of which 79 were confirmed to confer AFB1 resistance. While functionally diverse genes, including those that function in proteolysis, actin reorganization, and tRNA modification, were identified, those that function in postreplication DNA repair and encode proteins that bind to DNA damage were over-represented, compared to the yeast genome, at large. DNA metabolism genes also included those functioning in checkpoint recovery and replication fork maintenance, emphasizing the potency of the mycotoxin to trigger replication stress. Among genes involved in postreplication repair, we observed that CSM2, a member of the CSM2(SHU) complex, functioned in AFB1-associated sister chromatid recombination while suppressing AFB1-associated mutations. These studies thus broaden the number of AFB1 resistance genes and have elucidated a mechanism of error-free bypass of AFB1-associated DNA adducts.

Structural rearrangements like copy number variations in the male-specific Y chromosome have been associated with male fertility phenotypes in human and mouse but have been sparsely studied in other mammalian species. Here, we designed digital droplet PCR assays for 7 horse male-specific Y chromosome multicopy genes and SRY and evaluated their absolute copy numbers in 209 normal male horses of 22 breeds, 73 XY horses with disorders of sex development and/or infertility, 5 Przewalski's horses and 2 kulans. This established baseline copy number for these genes in horses. The TSPY gene showed the highest copy number and was the most copy number variable between individuals and breeds. SRY was a single-copy gene in most horses but had 2-3 copies in some indigenous breeds. Since SRY is flanked by 2 copies of RBMY, their copy number variations were interrelated and may lead to SRY-negative XY disorders of sex development. The Przewalski's horse and kulan had 1 copy of SRY and RBMY. TSPY and ETSTY2 showed significant copy number variations between cryptorchid and normal males (P < 0.05). No significant copy number variations were observed in subfertile/infertile males. Notably, copy number of TSPY and ETSTY5 differed between successive male generations and between cloned horses, indicating germline and somatic mechanisms for copy number variations. We observed no correlation between male-specific Y chromosome gene copy number variations and male-specific Y chromosome haplotypes. We conclude that the ampliconic male-specific Y chromosome reference assembly has deficiencies and further studies with an improved male-specific Y chromosome assembly are needed to determine selective constraints over horse male-specific Y chromosome gene copy number and their relation to stallion reproduction and male biology.

The double-stranded RNA-binding protein Staufen has been implicated in various posttranscriptional gene regulatory processes. Here, we demonstrate that the Caenorhabditis elegans homolog of Staufen, STAU-1, functionally interacts with microRNAs. Loss-of-function mutations of stau-1 significantly suppress phenotypes of let-7 family microRNA mutants, a hypomorphic allele of dicer, and a lsy-6 microRNA partial loss-of-function mutant. Furthermore, STAU-1 modulates the activity of lin-14, a target of lin-4 and let-7 family microRNAs, and this modulation is abolished when the 3' untranslated region of lin-14 is removed. Deep sequencing of small RNA cDNA libraries reveals no dramatic change in the levels of microRNAs or other small RNA populations between wild-type and stau-1 mutants, with the exception of certain endogenous siRNAs in the WAGO pathway. The modulation of microRNA activity by STAU-1 does not seem to be associated with the previously reported enhanced exogenous RNAi (Eri) phenotype of stau-1 mutants, since eri-1 exhibits the opposite effect on microRNA activity. Altogether, our results suggest that STAU-1 negatively modulates microRNA activity downstream of microRNA biogenesis, possibly by competing with microRNAs for binding on the 3' untranslated region of target mRNAs.

Estimation of genetic population structure based on molecular markers is a common task in population genetics and ecology. We apply a generalized linear model with LASSO regularization to infer relationships between individuals and populations from molecular marker data. Specifically, we apply a neighborhood selection algorithm to infer population genetic structure and gene flow between populations. The resulting relationships are used to construct an individual-level population graph. Different network substructures known as communities are then dissociated from each other using a community detection algorithm. Inference of population structure using networks combines the good properties of: (i) network theory (broad collection of tools, including aesthetically pleasing visualization), (ii) principal component analysis (dimension reduction together with simple visual inspection), and (iii) model-based methods (e.g., ancestry coefficient estimates). We have named our process CONE (for community oriented network estimation). CONE has fewer restrictions than conventional assignment methods in that properties such as the number of subpopulations need not be fixed before the analysis and the sample may include close relatives or involve uneven sampling. Applying CONE on simulated data sets resulted in more accurate estimates of the true number of subpopulations than model-based methods, and provided comparable ancestry coefficient estimates. Inference of empirical data sets of teosinte single nucleotide polymorphism, bacterial disease outbreak, and the human genome diversity panel illustrate that population structures estimated with CONE are consistent with the earlier findings.

Malaria continues to be a major global health problem, where disease transmission is deeply linked to the repeated blood feeding nature of the anautogenous mosquito. Given the tight link between blood feeding and disease transmission, understanding basic biology behind mosquito physiology is a requirement for developing effective vector-borne disease control strategies. In the mosquito, numerous loss of function studies with notable phenotypes demonstrate microRNAs (miRNAs) play significant roles in mosquito physiology. While the field appreciates the importance of a handful of miRNAs, we still need global mosquito tissue miRNA transcriptome studies. To address this need, our goal was to determine the miRNA transcriptome for multiple tissues of the pre-vitellogenic mosquito. To this end, by using small RNA-Seq analysis, we determined miRNA transcriptomes in tissues critical for mosquito reproduction and immunity including (i) fat body-abdominal wall enriched tissues, (ii) midguts, (iii) ovaries, and (iv) remaining tissues comprised of the head and thorax. We found numerous examples of miRNAs exhibiting pan-tissue high- or low- expression, tissue exclusion, and tissue enrichment. We also updated and consolidated the miRNA catalog and provided a detailed genome architecture map for the malaria vector, Anopheles gambiae This study aims to build a foundation for future research on how miRNAs and potentially other small RNAs regulate mosquito physiology as it relates to vector-borne disease transmission.

In flowering plants, gene body methylation (gbM) is associated with a subset of constitutively expressed genes. It has been proposed that gbM modulates gene expression. Here, we show that there are no consistent and direct differences to expression following the loss of gbM. By comparing expression of gbM genes in Arabidopsis thaliana accessions to orthologous genes in two Eutrema salsugineum genotypes, we identified both positive and negative expression differences associated with gbM loss. However, expression is largely unaffected by gbM loss in E. salsugineum Expression differences between species were within the variation of expression observed within A. thaliana accessions that displayed variation in gbM. Furthermore, experimentally induced loss of gbM did not consistently lead to differences in expression compared to wild type. To date, there is no convincing data to support a direct causal link between the presence/absence of gbM and the modulation of expression in flowering plants.

Conventional karyotypes and various genetic linkage maps have been established in sunflower (Helianthus annuus L., 2n = 34). However, the relationship between linkage groups and individual chromosomes of sunflower remains unknown and has considerable relevance for the sunflower research community. Recently, a set of linkage group-specific bacterial /binary bacterial artificial chromosome (BAC/BIBAC) clones was identified from two complementary BAC and BIBAC libraries constructed for cultivated sunflower cv. HA89. In the present study, we used these linkage group-specific clones (~100 kb in size) as probes to in situ hybridize to HA89 mitotic chromosomes at metaphase using the BAC-fluorescence in situ hybridization (FISH) technique. Because a characteristic of the sunflower genome is the abundance of repetitive DNA sequences, a high ratio of blocking DNA to probe DNA was applied to hybridization reactions to minimize the background noise. As a result, all sunflower chromosomes were anchored by one or two BAC/BIBAC clones with specific FISH signals. FISH analysis based on tandem repetitive sequences, such as rRNA genes, has been previously reported; however, the BAC-FISH technique developed here using restriction fragment length polymorphism (RFLP)-derived BAC/BIBAC clones as probes to apply genome-wide analysis is new for sunflower. As chromosome-specific cytogenetic markers, the selected BAC/BIBAC clones that encompass the 17 linkage groups provide a valuable tool for identifying sunflower cytogenetic stocks (such as trisomics) and tracking alien chromosomes in interspecific crosses. This work also demonstrates the potential of using a large-insert DNA library for the development of molecular cytogenetic resources.

The cost-effectiveness of sequencing pools of individuals (Pool-Seq) provides the basis for the popularity and widespread use of this method for many research questions, ranging from unraveling the genetic basis of complex traits, to the clonal evolution of cancer cells. Because the accuracy of Pool-Seq could be affected by many potential sources of error, several studies have determined, for example, the influence of sequencing technology, the library preparation protocol, and mapping parameters. Nevertheless, the impact of the mapping tools has not yet been evaluated. Using simulated and real Pool-Seq data, we demonstrate a substantial impact of the mapping tools, leading to characteristic false positives in genome-wide scans. The problem of false positives was particularly pronounced when data with different read lengths and insert sizes were compared. Out of 14 evaluated algorithms novoalign, bwa mem and clc4 are most suitable for mapping Pool-Seq data. Nevertheless, no single algorithm is sufficient for avoiding all false positives. We show that the intersection of the results of two mapping algorithms provides a simple, yet effective, strategy to eliminate false positives. We propose that the implementation of a consistent Pool-Seq bioinformatics pipeline, building on the recommendations of this study, can substantially increase the reliability of Pool-Seq results, in particular when libraries generated with different protocols are being compared.

The rapid evolution of insecticide resistance remains one of the biggest challenges in the control of medically and economically important pests. Insects have evolved a diverse range of mechanisms to reduce the efficacy of the commonly used classes of insecticides, and finding the genetic basis of resistance is a major aid to management. In a previously unstudied population, we performed an F2 resistance mapping cross for the common bed bug, Cimex lectularius, for which insecticide resistance is increasingly widespread. Using 334 SNP markers obtained through RAD-sequencing, we constructed the first linkage map for the species, consisting of 14 putative linkage groups (LG), with a length of 407 cM and an average marker spacing of 1.3 cM. The linkage map was used to reassemble the recently published reference genome, facilitating refinement and validation of the current genome assembly. We detected a major QTL on LG12 associated with insecticide resistance, occurring in close proximity (1.2 Mb) to a carboxylesterase encoding candidate gene for pyrethroid resistance. This provides another example of this candidate gene playing a major role in determining survival in a bed bug population following pesticide resistance evolution. The recent availability of the bed bug genome, complete with a full list of potential candidate genes related to insecticide resistance, in addition to the linkage map generated here, provides an excellent resource for future research on the development and spread of insecticide resistance in this resurging pest species.

Metazoan embryos undergo a maternal-to-zygotic transition (MZT) during which maternal gene products are eliminated and the zygotic genome becomes transcriptionally active. During this process, RNA-binding proteins (RBPs) and the microRNA-induced silencing complex (miRISC) target maternal mRNAs for degradation. In Drosophila, the Smaug (SMG), Brain tumor (BRAT), and Pumilio (PUM) RBPs bind to and direct the degradation of largely distinct subsets of maternal mRNAs. SMG has also been shown to be required for zygotic synthesis of mRNAs and several members of the miR-309 family of microRNAs (miRNAs) during the MZT. Here, we have carried out global analysis of small RNAs both in wild-type and in smg mutants. Our results show that 85% of all miRNA species encoded by the genome are present during the MZT. Whereas loss of SMG has no detectable effect on Piwi-interacting RNAs (piRNAs) or small interfering RNAs (siRNAs), zygotic production of more than 70 species of miRNAs fails or is delayed in smg mutants. SMG is also required for the synthesis and stability of a key miRISC component, Argonaute 1 (AGO1), but plays no role in accumulation of the Argonaute family proteins associated with piRNAs or siRNAs. In smg mutants, maternal mRNAs that are predicted targets of the SMG-dependent zygotic miRNAs fail to be cleared. BRAT and PUM share target mRNAs with these miRNAs but not with SMG itself. We hypothesize that SMG controls the MZT, not only through direct targeting of a subset of maternal mRNAs for degradation but, indirectly, through production and function of miRNAs and miRISC, which act together with BRAT and/or PUM to control clearance of a distinct subset of maternal mRNAs.

Cucumber (Cucumis sativus L.) has a large, paternally transmitted mitochondrial genome. Cucumber plants regenerated from cell cultures occasionally show paternally transmitted mosaic (MSC) phenotypes, characterized by slower growth, chlorotic patterns on the leaves and fruit, lower fertility, and rearrangements in their mitochondrial DNAs (mtDNAs). MSC lines 3, 12, and 16 originated from different cell cultures all established using the highly inbred, wild-type line B. These MSC lines possess different rearrangements and under-represented regions in their mtDNAs. We completed RNA-seq on normalized and non-normalized cDNA libraries from MSC3, MSC12, and MSC16 to study their nuclear gene-expression profiles relative to inbred B. Results from both libraries indicated that gene expression in MSC12 and MSC16 were more similar to each other than MSC3. Forty-one differentially expressed genes (DEGs) were upregulated and one downregulated in the MSC lines relative to B. Gene functional classifications revealed that more than half of these DEGs are associated with stress-response pathways. Consistent with this observation, we detected elevated levels of hydrogen peroxide throughout leaf tissue in all MSC lines compared to wild-type line B. These results demonstrate that independently produced MSC lines with different mitochondrial polymorphisms show unique and shared nuclear responses. This study revealed genes associated with stress response that could become selection targets to develop cucumber cultivars with increased stress tolerance, and further support of cucumber as a model plant to study nuclear-mitochondrial interactions.

The basidiomycete Moniliophthora roreri causes frosty pod rot of cacao (Theobroma cacao) in the western hemisphere. Moniliophthora roreri is considered asexual and haploid throughout its hemibiotrophic life cycle. To understand the processes driving genome modification, using long-read sequencing technology, we sequenced and assembled 5 high-quality M. roreri genomes out of a collection of 99 isolates collected throughout the pathogen's range. We obtained chromosome-scale assemblies composed of 11 scaffolds. We used short-read technology to sequence the genomes of 22 similarly chosen isolates. Alignments among the 5 reference assemblies revealed inversions, translocations, and duplications between and within scaffolds. Isolates at the front of the pathogens' expanding range tend to share lineage-specific structural variants, as confirmed by short-read sequencing. We identified, for the first time, 3 new mating type A locus alleles (5 in total) and 1 new potential mating type B locus allele (3 in total). Currently, only 2 mating type combinations, A1B1 and A2B2, are known to exist outside of Colombia. A systematic survey of the M. roreri transcriptome across 2 isolates identified an expanded candidate effector pool and provided evidence that effector candidate genes unique to the Moniliophthoras are preferentially expressed during the biotrophic phase of disease. Notably, M. roreri isolates in Costa Rica carry a chromosome segment duplication that has doubled the associated gene complement and includes secreted proteins and candidate effectors. Clonal reproduction of the haploid M. roreri genome has allowed lineages with unique genome structures and compositions to dominate as it expands its range, displaying a significant founder effect.

The Western European house mouse (Mus musculus domesticus) is a widespread human commensal that has recently been introduced to North America. Its introduction to the Americas is thought to have resulted from the transatlantic movements of Europeans that began in the early 16th century. To study the details of this colonization history, we examine population structure, explore relevant demographic models, and infer the timing of divergence among house mouse populations in the eastern United States using published exome sequences from five North American populations and two European populations. For North American populations of house mice, levels of nucleotide variation were lower, and low-frequency alleles were less common than for European populations. These patterns provide evidence of a mild bottleneck associated with the movement of house mice into North America. Several analyses revealed that one North American population is genetically admixed, which indicates at least two source populations from Europe were independently introduced to eastern North America. Estimated divergence times between North American and German populations ranged between ∼1,000 and 7,000 years ago and overlapped with the estimated divergence time between populations from Germany and France. Demographic models comparing different North American populations revealed that these populations diverged from each other mostly within the last 500 years, consistent with the timing of the arrival of Western European settlers to North America. Together, these results support a recent introduction of Western European house mice to eastern North America, highlighting the effects of human migration and colonization on the spread of an invasive human commensal.

Traits involved in reproduction evolve rapidly and show great diversity among closely related species. However, the genetic mechanisms that underlie the diversification of courtship traits are mostly unknown. Japanese medaka fishes (Oryzias latipes) use anal fins to attract females and to grasp females during courtship; the males have longer anal fins with male-specific ossified papillary processes on the fin rays. However, anal fin morphology varies between populations: the southern populations tend to have longer anal fins and more processes than the northern populations. In the present study, we conducted quantitative trait locus (QTL) mapping to investigate the genetic architecture underlying the variation in the number of papillary processes of Japanese medaka fish and compared the QTL with previously identified QTL controlling anal fin length. First, we found that only a few QTL were shared between anal fin length and papillary process number. Second, we found that the numbers of papillary processes on different fin rays often were controlled by different QTL. Finally, we produced another independent cross and found that some QTL were repeatable between the two crosses, whereas others were specific to only one cross. These results suggest that variation in the number of papillary processes is polygenic and controlled by QTL that are distinct from those controlling anal fin length. Thus, different courtship traits in Japanese medaka share a small number of QTL and have the potential for independent evolution.

While quantitative PCR (qPCR) is widely recognized as being among the most accurate methods for quantifying gene expression, it is highly dependent on the use of reliable, stably expressed reference genes. With the increased availability of high-throughput methods for measuring gene expression, whole-transcriptome approaches may be increasingly utilized for reference gene selection and validation. In this study, RNA-seq was used to identify a set of novel qPCR reference genes and evaluate a panel of traditional "housekeeping" reference genes in two species of the evolutionary model plant genus Mimulus More broadly, the methods proposed in this study can be used to harness the power of transcriptomes to identify appropriate reference genes for qPCR in any study organism, including emerging and nonmodel systems. We find that RNA-seq accurately estimates gene expression means in comparison to qPCR, and that expression means are robust to moderate environmental and genetic variation. However, measures of expression variability were only in agreement with qPCR for samples obtained from a shared environment. This result, along with transcriptome-wide comparisons, suggests that environmental changes have greater impacts on expression variability than on expression means. We discuss how this issue can be addressed through experimental design, and suggest that the ever-expanding pool of published transcriptomes represents a rich and low-cost resource for developing better reference genes for qPCR.