Author Response: Genetic Basis of Srna Quantitative Variation Analyzed Using an Experimental Population Derived from an Elite Rice Hybrid
Jia Wang,Wen Yao,Dan Zhu,Weibo Xie,Qifa Zhang
DOI: https://doi.org/10.7554/elife.03913.036
2015-01-01
Abstract:Article Figures and data Abstract eLife digest Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract We performed a genetic analysis of sRNA abundance in flag leaf from an immortalized F2 (IMF2) population in rice. We identified 53,613,739 unique sRNAs and 165,797 sRNA expression traits (s-traits). A total of 66,649 s-traits mapped 40,049 local-sQTLs and 30,809 distant-sQTLs. By defining 80,362 sRNA clusters, 22,263 sRNA cluster QTLs (scQTLs) were recovered for 20,249 of all the 50,139 sRNA cluster expression traits (sc-traits). The expression levels for most of s-traits from the same genes or the same sRNA clusters were slightly positively correlated. While genetic co-regulation between sRNAs from the same mother genes and between sRNAs and their mother genes was observed for a portion of the sRNAs, most of the sRNAs and their mother genes showed little co-regulation. Some sRNA biogenesis genes were located in distant-sQTL hotspots and showed correspondence with specific length classes of sRNAs suggesting their important roles in the regulation and biogenesis of the sRNAs. https://doi.org/10.7554/eLife.03913.001 eLife digest Genes within the DNA of a plant or animal contain instructions to make molecules called RNAs. Some RNA molecules can be decoded to make proteins, whereas others have different roles. A single gene often contains the instructions to make both protein-coding RNAs and non-coding RNAs. Molecules called small RNAs (or sRNAs) do not code for proteins. Instead, sRNAs can control protein-coding RNA molecules or chemically alter the DNA itself; this allows them to perform many different roles in living organisms. In plants, for example, these molecules affect how the plant grows, the shapes and structures it forms, and how likely it is to survive challenges such as drought and diseases. Often different plants of the same species have different amounts of sRNAs, but the reasons for this remain unclear. Now, Wang, Yao et al. have made use of a technique called ‘expression quantitative locus’ analysis to look at how sRNAs in rice plants are controlled by additional information encoded within DNA. The analysis identified over 53 million sRNA molecules from a population of rice plants. Many of these sRNAs varied in their abundance between different plants within the population. Wang, Yao et al. also found many thousands of individual instructions within the DNA of the rice that can either increase or reduce the abundance of their associated sRNA. Some of the abundant sRNAs were influenced by instructions within their own genes; some were influenced by instructions from other genes; and some were influenced by both. Wang, Yao et al. also found that the control of protein-coding RNAs was not necessarily related to the control of sRNAs encoded by the same gene. Further work is now needed to identify which specific DNA sequences regulate the abundance of sRNA molecules in plants and other organisms. https://doi.org/10.7554/eLife.03913.002 Introduction Small RNAs (sRNAs) are non-coding RNAs mainly 18–30 nt in length that regulate a wide range of biological processes in eukaryotic organisms (Carthew and Sontheimer, 2009; Axtell, 2013). According to their origin, sRNAs can be grouped into two major types: hpRNAs that are derived from single-stranded precursors with a hairpin structure (such as microRNAs [miRNAs]) and short interfering RNA (siRNAs) that are derived from double-stranded RNA precursors such as heterochromatic small interfering RNAs (hc-siRNA) and trans-acting siRNAs (ta-siRNA). There has been an explosion of interest in recent years in studies of miRNAs and siRNAs on their identification, biogenesis, and functioning in diverse biological processes. In plants, sRNAs function in regulating growth, development (Juarez et al., 2004; Zhu and Helliwell, 2011), architecture (Jiao et al., 2010; Miura et al., 2010), yield (Zhang et al., 2013), and response to biotic and abiotic stresses (Lu et al., 2008; Shukla et al., 2008). Such regulations are usually achieved by mediating endogenous mRNA cleavage and decay, DNA methylation of source and target loci, and chromatin modification and transcriptional silencing (Arikit et al., 2013). Although the sRNAs differ in length, sequences, and functions, the pathways for their biogenesis and functioning from precursor transcription, processing, maturation, and action are relatively conserved, which involve the activities of a number of enzymes including RNA polymerase II (Pol II), RNA-dependent RNA polymerases (RDRs), Dicer-like proteins (DCLs), and Argonautes (AGOs) (Chen, 2009; Ghildiyal and Zamore, 2009). It is also known that the abundance of sRNA species can be highly variable among individuals within a species, and the possible regulatory role of such quantitative difference has been assumed (He et al., 2010; Groszmann et al., 2011b). It is not known whether the quantitative variation of sRNA species between genotypes is related to the biological machinery, as quantitative variation of sRNAs has not been assayed at the population level and their genetic control has yet to be elucidated. The recently developed expression quantitative trait locus (eQTL) analysis has provided an approach for determining the genetic control of the expression level of a gene, including cis- and trans-eQTLs, as well as epistatic effects (Becker et al., 2012). This approach can also be applied to the genetic analysis of quantitative variation of sRNAs by regarding the abundance of the sRNAs in the population as quantitative traits. Once the QTLs are identified, subsequent studies can be pursued very much the same way as the analysis of genes and regulatory networks underpinning phenotypic QTLs (Xing and Zhang, 2010). There have been several studies focusing on the genetic regulation of known and validated miRNAs and small nucleolar RNAs in specific tissues/cells from samples of the human population (Borel et al., 2011; Gamazon et al., 2012; Parts et al., 2012; Civelek et al., 2013; Jin and Lee, 2013; Siddle et al., 2014). These studies detected a number of cis- and trans-miQTLs that could also influence the expression of the mRNA targets, which may be associated with phenotype difference. Here, we performed a whole genome QTL analysis of the entire sRNA kingdom consisted of 18-nt to 26-nt sRNAs from flag leaf of rice using an experimental genetic population. The analysis revealed features of the genetic controls of sRNA abundance showing both commonality and distinction with their precursor transcripts. It was also shown that the abundance of sRNAs is probably related to proteins constituting the machinery for sRNA biogenesis and functioning. Results Patterns and distributions of sRNAs The genetic materials consisted of 98 hybrids obtained by paired crosses of 196 recombinant inbred lines (RILs) derived by single seed descent from a cross between Zhenshan 97 and Minghui 63, the parents of Shanyou 63 that was the most widely cultivated rice hybrid in the 1980s and 1990s and still used with reduced area in recent years in China. Because the genetic composition of these crosses resembled individuals in an F2 population, they were referred to as an immortalized F2 (IMF2) population (Hua et al., 2002, 2003), and each hybrid in the population was hereafter referred to as an IMF2 for ease of description. An sRNA library was constructed using RNA extracted from flag leaf at the day of full expansion for each IMF2, and two biological replicates were obtained for each of the two parental lines and their F1 hybrid, producing a total of 104 libraries. The read size of the raw sequencing data obtained using Illumina Hiseq2000 varied from 3 nt to 44 nt (Figure 1—figure supplement 1). Sequences of 18–26 nt in length that appeared to be the more abundant than others were kept for the analyses after filtering out low-quality reads and eliminating ones matching tRNAs, rRNAs, snRNAs, and snoRNAs (Figure 1—figure supplement 2). The numbers of resulting reads varied from 14.52 million to 27.73 million per library (2.08 billion in total) with non-redundant reads ranging from 3.00 million to 6.86 million per library (0.52 billion in total) (Supplementary file 1 in Dryad [Wang et al., 2015], Table 1, Figure 1—figure supplement 2). The 24 nt sRNAs were the most numerous in both redundant and distinct reads (Figure 1—figure supplement 2). Table 1 Average number of reads obtained for the 104 libraries from flag leaves of the IMF2 population and the parental lines https://doi.org/10.7554/eLife.03913.003 Quality filtering*Size filtering†NcRNA filtering‡Genome mapping§Total#Redundant reads22,864,06722,318,92819,968,5669,816,9331,020,961,021Distinct reads5,106,6545,035,0734,986,1842,499,23353,613,739 * Reads after filtering out low-quality reads. † Reads of 18–26 nt in length. ‡ Reads after eliminating ones matching tRNAs, rRNAs, snRNAs, and snoRNAs. § Reads mapped to the SNP-replaced reference genomes of the parents with unique locations allowing no mismatch. # The total reads of 104 libraries used to identify s-traits. We mapped the reads to the SNP-replaced reference genomes of the parents (‘Materials and methods’), with unique location allowing no mismatch. The reads could be divided into four categories: (1) approximately 82.28% of the sRNAs had identical sequences between parents, (2) 0.71% of the sRNAs had SNPs between the two parents, (3) 8.35% were only mapped to the Zhenshan 97 genome, and (4) 8.66% were specifically mapped to the Minghui 63 genome (Supplementary file 2 in Dryad [Wang et al., 2015]). A total of 53,613,739 unique sRNA sequences including ones with SNPs between the parents were identified by combining all the reads from 104 libraries (Table 1). Approximately 84.7% of sRNAs were found in no more than 5 IMF2s, while only 0.13% were present in all 98 IMF2s (Figure 1A). Figure 1 with 5 supplements see all Download asset Open asset The distribution of sRNAs and s-traits in different genomic regions across the IMF2 population. (A) The distribution of sRNAs across the IMF2 population. (B) The distribution of sRNAs aligned to the genic region, 2-kb upstream, and 500-bp downstream of annotated genes, as well as intergenic regions. TE: transposons; NTE: non-transposon genes. (C) The percentages of sRNAs of different sizes. (D) The distribution of sRNAs of different sizes in different genomic regions. (E) The distribution of sRNAs of different sizes in different portions of genic regions. (F) The percentages of s-traits of different sizes. (G) The distribution of s-traits aligned to the genic region, 2-kb upstream, and 500-bp downstream of annotated genes, as well as intergenic regions. The legend is the same as in (B). (H) The distribution of s-traits of different size in different genomic regions. The legend is the same as in (D). (I) The distribution of s-traits of different sizes in different portions of genic regions. The legend is the same as in (E). https://doi.org/10.7554/eLife.03913.004 The distributions of the sRNAs in different portions of the genome were not random. Approximately 51.1% of sRNAs originated from the 2-kb upstream and the genic regions of non-transposon genes (Figure 1B). In addition, sRNAs occurred in high frequencies near the transcription start sites compared to other regions of the promoters (Figure 1—figure supplement 3). About 34.4% of sRNAs originated from intergenic regions, which account for no more than 20% of the whole genome length (Figure 1B). About 43.9% of sRNAs were 24 nt, while 12.1% were 21 nt (Figure 1C). Different species of sRNAs differed in their origins of genomic regions. The genic regions of non-transposon genes held the highest number of 21 nt sRNAs, while 24 nt sRNAs were most common in the intergenic regions (Figure 1D). sRNAs of 21 nt in genic regions were mostly derived from the exon of non-transposon genes, while 24 nt sRNAs in genic regions were mainly from the intron of non-transposon genes (Figure 1E). The 24-nt sRNAs mainly consist of endogenous heterochromatic siRNAs (Axtell, 2013) and tend to be produced from repeats and transposable elements as well as intergenic regions (Chen, 2009). The 24-nt sRNAs are components of the epigenome that target the homologous genomic regions for de novo DNA methylation through RNA-directed DNA methylation to maintain genome stability by transcriptional gene silencing (Groszmann et al., 2011a, 2013; Pooggin, 2013). Next, the abundance of each sRNA in a library was normalized to number of reads per millions (RPM) (‘Materials and methods’) to quantify the sRNA levels in each library, which were subject to analyses and comparisons. An sRNA with RPM value ≥0.6 was regarded as expressed, and an sRNA identified as expressed in more than 25 of the 98 IMF2s was regarded as an sRNA expression trait (s-trait). In this way, a total of 165,797 s-traits were recovered (Supplementary file 3 in Dryad [Wang et al., 2015]). A total of 87.9% of s-traits were 24 nt (Figure 1F). The s-traits mainly originated from the intergenic regions, 2-kb upstream, and genic regions of non-transposon genes (Figure 1G). Although the 2-kb upstream and the genic regions of non-transposon genes contained almost equal number of sRNAs (Figure 1B), the number of s-traits derived from the 2-kb upstream of non-transposon genes was more than twice of the number of s-traits in genic regions of non-transposon genes (Figure 1G). The distribution of different species s-traits in different genomic regions was similar to that of sRNAs (Figure 1D,H). Although sRNAs of various lengths from genic regions were mostly from exons of non-transposon genes (Figure 1E), the s-traits of different sizes were mainly derived from introns of non-transposon genes (Figure 1I). The coefficients of variation of the 165,797 s-traits (Figure 1—figure supplement 4) showed a very wide range of distribution of the sRNA levels in the population. We selected 9 s-traits randomly from all the 165,797 s-traits, and the distributions of the expression values for most of these s-traits across the IMF2 population were more or less normal (Figure 1—figure supplement 5). Expression correlations between sRNAs that originated from the same genes To test whether s-traits originated from the same gene were transcribed together, the correlation coefficients between the expression values of s-traits originated from the same gene were calculated. A total of 3,739,873 correlation coefficients involving 15,078 genes and 96,914 s-traits were obtained, 2,593,503 (69.3%) and 507,528 (13.6%) of which were contributed by 2278 s-traits from genes LOC_Os03g01360 and 1008 s-traits from LOC_Os07g01240, respectively (Figure 2). The correlation coefficients naturally fell into three classes: strong negative correlations, no correlation, and strong positive correlations with obvious dividing points at −0.3 and 0.3 (Figure 2A), which were far above 0.25, the threshold for statistical significance at p < 0.05 determined by simulation data. While 85.8% of these strong correlations (positive or negative) were contributed by s-traits originated from LOC_Os03g01360 (Figure 2B), most of s-traits originated from gene LOC_Os07g01240 showed low (or no) correlations with each other (Figure 2C). Although the expression values of s-traits originated from the same regions of the remaining 15,076 genes in the whole genome were mostly slightly positively correlated, the correlation coefficients were significantly distinct from that of the simulation of random data (Mann–Whitney U test, p-value <2.2e-16, Figure 2D), indicating various degrees of correlations. We found that the correlations were affected by the sizes of s-traits (Figure 2—figure supplement 1). Expression correlations between s-traits of the same size were stronger than that between s-traits of different sizes. In particular, the expression correlations between 21-nt sRNAs and 25-nt sRNAs from the same genes were very weak. The correlations between overlapped s-traits were slightly stronger than that between s-traits not overlapping with each other (Figure 2—figure supplement 2A,B). The correlations between s-traits from the same transposons were slightly different from that of s-traits from the same non-transposons genes (Figure 2—figure supplement 2C,D). Most of the s-traits originating from LOC_Os03g01360 were significantly correlated (positive or negative) with each other (Figure 2—figure supplement 3), and the correlations between different s-traits were mainly dependent on their genomic positions rather than their sizes, which were mostly 21, 22, and 24 nt, accounting for 40.2%, 34.4%, and 19.4% of the 2278 s-traits. This gene could be separated as three parts, the 2 kb-upstream, 5′ UTR, the first intron and the first coding DNA sequence (CDS) as the first part, the most of the second intron comprising the second part, and the second CDS and part of the second intron as the third part (Figure 2—figure supplement 3). s-traits within each part were positively correlated with each other, while s-traits from the first part were mostly negatively correlated with those from the other two parts. The correlations between the expression levels of s-traits from the second parts and s-traits from the third parts were mostly positive. Thus, although very strong positive and negative correlations were found between sRNAs from the same genic regions of specific genes, this was not a common phenomenon. Taken together, this correlation analysis indicated that most of the sRNAs derived from the same genes were correlated with each other to various extents, implying that nearby sRNAs were possibly jointly transcribed. Consecutively overlapped sRNAs were found from introns of specific genes and intergenic regions, indicating that sRNAs were not necessarily jointly transcribed from the genes. Figure 2 with 5 supplements see all Download asset Open asset Expression correlations between s-traits originating from the same genes. (A) The distribution of expression correlations between s-traits originating from the same genes. (B) The distribution of expression correlations between s-traits originating from LOC_Os03g01360. (C) The distribution of expression correlations between s-traits originating from LOC_Os07g01240. (D) Comparison of the distribution of expression correlations for the simulation data and the true data excluding s-traits from LOC_Os03g01360 and LOC_Os07g01240. https://doi.org/10.7554/eLife.03913.010 Expression correlations between sRNAs that originated from the same sRNA cluster As many sRNAs may be derived from a single primary transcript, we defined sRNA clusters to quantify sRNA expression (Castel and Martienssen, 2013), based on our data of all 104 libraries. An sRNA island was defined as a genomic region composed of consecutive genomic positions matched by at least 30× sRNA read coverage. Nearby sRNA islands resided within 1000 nt from each other were merged and regarded as the same sRNA cluster. As a result, 80,362 sRNA clusters were recovered with an average size of 1927 bp (ranging from 51 bp to 47,021 bp) (Figure 2—figure supplement 4, Supplementary file 4 in Dryad [Wang et al., 2015]). A total of 46,293 (57.6%) sRNA clusters overlapped with genic regions, 23,143 (28.8%) of which were entirely derived from genic regions. In addition, the distances between another 19,164 (23.8%) sRNA clusters and the genic regions were smaller than 1 kb. The abundance of each sRNA cluster was normalized using DEseq (‘Materials and methods’) (Anders and Huber, 2010). An sRNA cluster with normalized expression value ≥6 was regarded as expressed, and an sRNA cluster expressed in more than 25 of all the 98 IMF2s was regarded as an sRNA cluster expression trait (sc-trait). A total of 50,139 (62.4%) sRNA clusters were determined as sc-traits, 34,762 of which were expressed in all of the 98 IMF2s. We next calculated the correlation coefficients between the expression values of s-traits originating from the same sRNA clusters to infer whether they were transcribed together. A total of 4,282,587 correlation coefficients involving 21,571 sRNA clusters and 157,994 s-traits were obtained, 72.5% of which were contributed by s-traits from clusters chr03-261504-273325 and chr07-136331-141317 (Figure 2—figure supplement 5). These two sRNA clusters corresponded to the two genes, LOC_Os03g01360 and LOC_Os07g01240, mentioned above. Of all the correlations, 2,776,296 were relatively strong (above 0.3 or below −0.3 based on a threshold of 0.25 determined by simulation data with p <0.05). However, 81.5% of these strong correlations were contributed by the 2280 s-traits from cluster chr03-261504-273325 (Figure 2—figure supplement 5B). A total of 1008 s-traits originated from cluster chr07-136331-141317, most of which were only slightly correlated with each other (Figure 2—figure supplement 5C). The distribution of correlation coefficients of s-traits from the same sRNA clusters was quite similar to that of s-traits from the same genes (Figure 2—figure supplement 5D). Although the s-traits originating from the same region of the remaining 215,69 sRNA clusters in the whole genome were only slightly positively correlated, the correlation coefficients were significantly different from that of simulated random data (Mann–Whitney U test, p-value <2.2e-16, Figure 2—figure supplement 5E). Thus, again very strong positive and negative expression correlations were only observed for sRNAs from certain sRNA clusters, while most sRNAs originating from the same sRNA clusters were only slightly positively correlated with each other. Expression correlations between sRNAs and their mother genes mRNA sequencing was conducted using the same samples as in the sRNA sequencing with one biological replicate for the IMF2 population and two biological replicates for each of the two parental lines and their F1 hybrid, producing a total of 104 libraries. After removal of low-quality sequencing data, reads were mapped to the Nipponbare reference genome using TopHat (Figure 3—figure supplement 1, Supplementary file 5 in Dryad [Wang et al., 2015]) (Trapnell et al., 2009). Most of the reads were mapped to the CDS and UTR regions of non-transposon genes (Figure 3—figure supplement 2). The abundance of each mRNA in a library was normalized to the number of fragments per kilobase of transcript per million mapped fragments (FPKM) using Cufflinks (Trapnell et al., 2010) to quantify the expression levels of mRNAs in each library. An mRNA with FPKM value ≥1 was regarded as expressed; 35,819 (54.2%) mRNAs were expressed in none of the 98 IMF2s, while 16,849 (25.5%) were expressed in all the 98 IMF2s. An mRNA identified as expressed in more than 25 of the 98 IMF2s was regarded as an e-trait. In this way, a total of 24,987 e-traits were obtained from a total of 66,123 mRNAs in the whole genome. For ease of description, we hereafter refer to genic sequences including both exons and introns that encoded the s-traits as ‘mother genes’. To assess whether sRNAs and their mother genes were transcribed together, the correlation coefficients of the expression values of s-traits and their mother genes were calculated. A total of 77,640 correlation coefficients were obtained, involving 9968 mRNAs and 56,185 s-traits. Among all these correlations, 13,749 (17.7%) were strong correlations (above 0.3 or below −0.3, based on a threshold of 0.23 determined by simulation data with p < 0.05). However, 45.6% of these strong correlations were again contributed by gene LOC_Os03g01360 (Figure 3A). Correlations involving LOC_Os07g01240 were mostly weak (Figure 3B). The correlation coefficients involving rest of the genes in the whole genome were centered around zero (Figure 3C). Figure 3 with 2 supplements see all Download asset Open asset Correlations between the s-traits and the transcripts of their mother genes. (A) Correlations between s-traits and e-traits derived from LOC_Os03g01360. (B) Correlations between s-traits and e-traits derived from LOC_Os07g01240. (C) Correlations between s-traits and e-traits originating from the same genes excluding the above two loci. (D) Correlations between s-traits and different e-traits derived from LOC_Os03g01360. https://doi.org/10.7554/eLife.03913.016 In-depth analysis of s-traits from gene LOC_Os03g01360 revealed that the expression values of s-traits from the first part of LOC_Os03g01360 were positively correlated with two transcripts, LOC_Os03g01360.1 and LOC_Os03g01360.4, and negatively correlated with the other LOC_Os03g01360.2 (Figure 3D). On the other hand, the expression values of most s-traits from the third part of LOC_Os03g01360 were negatively correlated with LOC_Os03g01360.1 and LOC_Os03g01360.4 but positively correlated with LOC_Os03g01360.2 (Figure 3D). Whereas both negative and positive expression correlations were found between s-traits from the second part of LOC_Os03g01360 and these three transcripts (Figure 3D). Thus, although strong expression correlations, either positive or negative, were observed between specific sRNAs and their mother gene, the expression levels of most sRNAs were not correlated with their mother genes, suggesting independent transcription of sRNAs and mRNAs. QTL analysis of s-traits, sc-traits, and e-traits using the IMF2 population QTL analysis was performed for the 165,797 s-traits based on the ultrahigh-density SNP map, which contained 1556 recombination events and was composed of 1568 bins with an average size of 238 kb (ranging from 6 kb to 7947 kb) (Figure 4—figure supplement 1, Figure 4—figure supplement 2, Supplementary file 6 in Dryad [Wang et al., 2015], ‘Materials and methods’) (Xie et al., 2010; Yu et al., 2011), using composite interval mapping (CIM) in R/qtl with 1000 permutations (Haley and Knott, 1992; Broman and Speed, 2002; Manichaikul et al., 2009). With a false discovery rate (FDR) set at 5%, a total of 70,858 sQTLs were recovered for 66,649 s-traits (Supplementary file 7 in Dryad [Wang et al., 2015]). These sQTLs were classified as local-sQTLs and distant-sQTLs, which were also referred to as cis- and trans-QTLs, respectively, in previous studies (Wang et al., 2010, 2014), according to their locations relative to the s-traits. A local-sQTL indicated the existence of local functional polymorphism(s) that could influence the abundance of the s-trait, and a distant-sQTL meant the s-trait expression variation in the population was controlled by regulatory element(s) distant from the s-trait precursor sequence. For defining a local-sQTL, we adopted a 1.5 LOD-drop support interval of the corresponding sQTL or no more than 250 kb from the closest marker in a recombination sparse region. Otherwise, it was regarded as a distant-sQTL (Yvert et al., 2003; Morley et al., 2004; Keurentjes et al., 2007). Such classification resulted in 40,049 local-sQTLs (shown in the diagonal of Figure 4) and 30,809 distant-sQTLs (off-diagonal of Figure 4) (Supplementary file 7 in Dryad [Wang et al., 2015]). Expression variations of s-traits explained by local-sQTLs were substantially higher than that by distant-sQTLs (Figure 4—figure supplement 3). The LOD values for 80.4% of the local-sQTLs and 34.5% of distant-sQTLs were larger than 10. Figure 4 with 8 supplements see all Download asset Open asset sQTLs for the 66,649 s-traits. The color key shows the LOD value. X-axis, the physical position of sQTLs along the genome of 12 chromosomes. Y-axis, the physical position of s-traits. QTLs with LOD value <5 are not included in the presentation. https://doi.org/10.7554/eLife.03913.019 Next, we performed QTL analysis for 50,139 sc-traits based on the same genetic map as in sQTL analysis using CIM in R/qtl with 1000 permutations. A total of 22,263 scQTLs including 11,476 local-scQTLs and 10,787 distant-sQTLs were recovered for 20,049 sc-traits utilizing the same definitions of local- and distant-QTLs (Figure 4—figure supplement 4, Supplementary file 7 in Dryad [Wang et al., 2015]). Again, the LOD values of local-scQTLs and expression variations explained by local-scQTLs were much higher than distant-scQTLs (Figure 4—figure supplement 5). We also performed QTL analysis for the 24,987 e-traits based on the same genetic map, using the same algorithm and parameters as in the QTL analysis of s-traits and sc-traits. With an FDR set at 5%, a total of 6423 eQTLs were recovered for 6123 e-traits (Supplementary file 7 in Dryad [Wang et al., 2015]). These eQTLs were classified as local-eQTLs and distant-eQTLs according to the same rule applied in the sQTL analysis. As a result, 2964 local-eQTLs (shown in the diagonal of Figure 4—figure supplement 6) and 3459 distant-eQTLs (off-diagonal of Figure 4—figure supplement 6) were detected (Supplementary file 7 in Dryad [Wang et al., 2015]). Expression variations of mRNAs explained by local-eQTLs were significantly higher than that by distant-eQTLs (Figure 4—figure supplement 7). Local-QTLs regulating the expression of large numbers of traits were observed for s-traits and sc-traits but not for e-traits (Figure 4, Figure 4—figure supplement 4 and Figure 4—figure supplement 6), probably due to the consecutive transcription of sRNAs and sRNA clusters, which was distinct from the transcription of mRNAs. In the extreme case, a total of 1660 s-traits were regulated by a single local-sQTL (Bin358)