De Novo mPing and Ping Insertion Sites.

The genomic locations of mPing elements in each RIL were determined with RelocaTE2 (25). A total of 87,450 mPing insertions were identified across the 272 RILs and simplified to 16,914 unique mPing loci in the population (Table 1 and SI Appendix, Supplementary Data 2). Of these 16,914 loci, 3% were parental while the remaining 16,448, or 97% of mPing loci in the population were nonparental and defined herein as de novo insertions. Among the de novo insertions, 88% (14,534 of 16,448) are unique to a single individual, indicating that most new insertions occurred during or after the F2 generations that produced the recombinant inbreds (Table 1). The rest of the de novo insertions (1,914 of 16,448) are shared among 2 to 145 RILs, representing either early insertions that occurred in F1 or late insertions that share target sites in the genome (SI Appendix, Fig. S2). Of those found in a single individual, 72% (10,527 of 14,534) are homozygous and likely transposed in early generations of RIL self-pollination while the remaining 28% (4,007 of 14,534) are heterozygous (Table 1) and represent either more recent mPing insertions or mPing alleles that are lethal or less fit when homozygous. To distinguish these possibilities, SNPs flanking heterozygous mPing loci were analyzed to determine genotypes of genomic regions underlying mPing insertions as recent insertions are more likely to be in homozygous (inbred) regions. Of the 4,007 heterozygous insertions, we were able to unambiguously genotype 2,964, of which 96% (2,864 of 2,964) were homozygous (SI Appendix, Table S2). These data suggest that the vast majority of heterozygous mPing insertions are late events. Finally, like previously reported de novo mPing insertions (7, 8), 45% of all insertions (both homozygous and heterozygous) identified in this study are within 5 kb upstream of protein-coding genes (SI Appendix, Fig. S3).

Table 1.

Classification of parental and de novo mPing, Ping, and Pong loci in RILs

Classification	No. of mPing loci	No. of Ping loci	No. of Pong loci
Parental	466	8	6
Shared	7	0	5
Unique HEG4	415	7	0
Unique NB	44	1	1
De novo	16,448	17	0
Shared	1,914	0	0
Unique homozygous	10,527	17	0
Unique heterozygous	4,007	0	0

Open in a separate window

Analysis of Ping elements identified only 17 nonparental (de novo) Pings, and over half of these (8 of 17) are within 5 kb of protein-coding genes (Table 1 and SI Appendix, Table S3). Among the population, RIL270 had two de novo Ping insertions (SI Appendix, Supplementary Data 2) while each of the 15 RILs had a single new Ping insertion. De novo Ping insertions were only found in RILs with two or more parental Ping loci (SI Appendix, Supplementary Data 2).

Correlation between Ping Copy Number and Number of New mPing Insertions.

Among the 272 RILs, significant variation (0 to 138) was found in the number of unique mPing insertions (from here on restricted to de novo unique homozygous mPing insertions) (SI Appendix, Supplementary Data 2), suggesting that multiple genetic loci control mPing transposition activity. To identify these loci, QTL mapping was performed using the number of unique mPing insertions as the trait (SI Appendix, Fig. S4). Three major loci were identified with the logarithm of the odds (LOD) scores greater than 3.31 and accounting for 49% of the phenotypic variation (SI Appendix, Fig. S5 and Table S4). Two of these loci contain multiple Ping elements (PingA and PingB on chromosome 1 and PingE, PingF, and PingG on chromosome 9), confirming prior data that all bursting accessions had multiple Pings (7–10 Pings) (9). A third QTL located on chromosome 4 (accounting for 11.82% of phenotypic variance) does not contain a Ping element (SI Appendix, Fig. S5 and Table S4), suggesting that additional factors may contribute to the rate of transposition. These potential host factors are beyond the scope of this study and will not be further discussed.

To date, quantification of the impact of Ping loci on mPing transposition has been restricted to accessions with either the single Ping locus in NB (PingH) or the collective impact of the 7 to 10 Ping loci in the bursting accessions (HEG4, EG4, A123, and A119) (9). In NB, mPing rarely transposes, and Ping transcript levels are very low (9). In contrast, there is significant amplification of mPing and transcription of Ping in the accessions with 7 to 10 Pings, including the RIL parent, HEG4 (9). This population allows for higher resolution analyses of Ping copy number and transposition and transcription activity because the eight parental Pings are segregating and, theoretically, should contain most combinations of Ping loci (SI Appendix, Fig. S6). In fact, several RILs were found to contain from 0 to 8 Ping loci, and a positive correlation was found to exist between Ping copy number and the number of new mPing insertions (Fig. 2A) (two-tailed Pearson’s correlation test, r = 0.71, P = 1e−42). For example, there are approximately 5 unique mPing insertions in RILs with one Ping whereas RIL with seven Ping elements have, on average, 65 unique mPing insertions. In addition, transcript levels of both ORF1 (Fig. 2B, blue) and TPase (Fig. 2B, red) increased linearly with Ping copy number (Fig. 2B; two-tailed Pearson’s correlation test, ORF1: r = 0.9, P = 6e−14; TPase: r = 0.85, P = 4e−11), suggesting that a simple dosage relationship exists between Ping expression and transposition activity.

Open in a separate window

Fig. 2.

Accumulation of unique mPing insertions is dependent on Ping dosage. (A) Correlation between unique mPing insertions and Ping copy number in the RILs. The 272 RILs were grouped by Ping copy number ranging from 0 to 7, and the scatterplot shows the number of unique mPing insertions in each group. The number of RILs in each group is in parentheses (n=). Green lines are the group median. The significance of correlation was tested by a two-tailed Pearson’s correlation test and is indicated by P value. (B) qRT-PCR analysis (three replicates) of Ping ORF1 and TPase transcription levels. RILs with Ping copy numbers of 1 to 7 were randomly selected. Relative transcription levels of ORF1 (blue) and TPases (red) were normalized with rice actin. Colored lines are the best-fit line of the linear regression model. The significance of correlation was tested by a two-tailed Pearson’s correlation test and is indicated by P values for both ORF1 (blue) and TPase (red). Sampled RILs with 1 to 7 Ping copies used in B were as follows: 1 (RIL12, -100, -166, -179); 2 (RIL16, -19, -36, -123); 3 (RIL5, -18, -22, -37, -111); 4 (RIL8, -13, -87, -119, -177, -219); 5 (RIL10, -15, -23, -58, -73, -92, -94); 6 (RIL30, -44, -54, -118, -134); and 7 (RIL11, -21, -34, -69, -198).

Ping Activity at Eight Different Loci.

The RIL population provided a unique opportunity to explore the activity of various combinations of the eight segregating Ping loci present in the two parents. This was of particular interest because, although all Ping elements are identical, one Ping locus has been implicated in initiating the bursts (9, 15). Specifically, PingA_Stow on chromosome 1 (called PingA in this study) was shown previously to be the only Ping shared by the four bursting accessions (9) and found preferentially in the genomes of rice accessions with higher than background mPing copies (15). Taken together, these correlative data suggest that the PingA locus differs from the other Ping loci, perhaps by catalyzing more transposition.

The activity of individual Ping loci was quantified by first identifying RIL with one Ping. Among the RIL population are 23 lines with a single Ping that, collectively, harbor all parental Pings except PingF. PingA is present as a single Ping locus in RIL60. Because mPing rarely transposes in accessions with one Ping (like NB), two high-resolution independent methodologies, transposon display and deep sequencing, were employed to assess the number of new mPing insertions.

To identify new mPing insertions in individual progeny, transposon display (TD) was performed using DNA isolated from eight plants derived by self-pollination from each analyzed single-Ping RIL (Fig. 3A and SI Appendix, Fig. S7 and Table S5). The mPing transposition activity was estimated by counting new amplicons present in one individual but absent from its siblings. New mPing insertions were detected in all single Ping RILs analyzed except RIL179 (PingB) (Fig. 3A); however, the progeny of RIL60 (with PingA), had significantly more (4 to 8 vs. 0 to 3 in the other single-Ping RIL) (Fig. 3A) (one-way ANOVA with Tukey’s honestly significant difference [HSD] test, P < 9.48e−8).

Open in a separate window

Fig. 3.

New mPing insertions in single-Ping RILs. (A) New mPing insertions estimated by transposon display. New insertions were estimated by counting the number of new bands displayed by single plants. Eight sibling plants from a single-seed descent were used for each selected single-Ping RIL. Error bars show SE of 2 to 3 independent biological replicates. (B) New mPing insertions estimated by a pooled-sequencing approach with multiple progeny from a single-Ping RIL. RIL179 (PingB) was excluded due to contamination. Error bars show SE of 3 to 4 independent biological replicates represented by black diamonds. Differences between PingA and other Ping loci were tested by a one-way ANOVA followed by Tukey’s HSD post hoc test.

In the second approach, new insertions were directly counted following deep sequencing of DNA samples isolated from eight-pooled progeny of each RIL. Detection of new mPing insertions for these RILs confirmed that all Pings tested are active, thus providing a simple explanation for the Ping dosage series reported above. Similarly, the finding that PingA promotes more transposition than the other Pings (Fig. 3B and SI Appendix, Table S6) (one-way ANOVA with Tukey’s HSD test, P < 3.09e−3) is consistent with the hypothesis that it initiated the bursts, but the underlying mechanism is beyond the scope of this study (Discussion).

Loss of Parental mPing Insertion Loci in the RIL Population.

Ping and mPing elements, like other DNA transposons, can both insert into new loci and excise from existing sites. By tracking the fate of the 466 parental loci in the RIL population (44 from NB, 415 from HEG4, and 7 shared) (Materials and Methods), we sought to determine whether some mPing loci were lost at a higher frequency than others. Loss of TEs over time has been attributed to many factors, the two most prominent being negative selection and structural instability (26, 27). To discriminate between these possibilities, we developed methodologies to identify all parental mPing excision events among the RIL population, determine excision frequencies at each locus, and analyze the genome context of loci that excised frequently. To this end, we characterized empty (excised) mPing sites in RILs within haplotype blocks classified by parent of origin and expected to contain parental mPing loci (Fig. 4A). On average, 127 RILs (ranging from 71 to 159) were analyzed for each parental mPing.

Open in a separate window

Fig. 4.

Identification of mPing loci with high excision frequencies. (A) Workflow for identification of mPing excisions (see details in Materials and Methods). Red and blue colors represent NB and HEG4 genotypes in parental accessions. mPing insertions (orange boxes) are exclusively in HEG4 genotypes. Sequencing reads aligned to pseudoreference sequences with or without the mPing insertion are indicated by light blue bars. (B) Example illustrating the method used to detect independent excision events at mPing locus Chr1:6816415 to 6816417. Red letters in NB indicate the 3-bp TSDs upon mPing insertion. Footprints (deletions and insertions) are indicated by dots or arrows. CIGAR-like strings are used to record the footprints. Excision footprints are characterized as deletions (D) or insertions (A) with the length (e.g., D-X bp) and, in parenthesis, the number of RILs with this footprint. (C) Distribution of parental mPing loci and number of independent excision events. The P value is based on a two-tailed binomial test. The dashed line is the cutoff for significant high frequency excision. The number of different mPing loci with that many excisions is at the top of each bar. (D) mPing loci with more excision events are in close proximity to another mPing. The box shows a zoom of the 0.04- to 100-kb range.

In total, 742 excisions were detected from 177 parental mPing loci (17 NB, 160 HEG4), and no excisions were detected from the remaining 289 loci (34 NB, 255 HEG4). To validate the accuracy of the computer-assisted excision data, 30 loci were randomly selected and assayed by PCR with DNA isolated from the RIL where the excision was detected. Using this approach, 26 of 30 events were confirmed (SI Appendix, Fig. S8). Next, to assess the excision frequency at each of the 177 parental loci, it was necessary to establish the independence of the excision events. For example, an excision event occurring early in the generation of the RIL population might be propagated in several RILs. Two approaches were used to eliminate dependent events: 1) a standard method where all excisions were analyzed after removal of early events, and 2) a conservative method where excision footprints were exploited to identify independent events. The latter is likely to underestimate independent events because perfect excision, which is common for mPing, will only be counted once per locus per RIL. Because both approaches produced very similar results, only the conservative method is described below. However, the results of both are found in Materials and Methods.

For the conservative method, 367 of the 742 excisions from parental mPing loci have excision footprints (SI Appendix, Supplementary Data 3). This reduced to 322 independent events following removal of excisions from the same locus with the same footprint (Fig. 4B and SI Appendix, Supplementary Data 3). Of the 177 loci, 114 had only a single excision, 51 had fewer than five (Fig. 4C), and 12 loci experienced five or more independent excision events (Fig. 4C and SI Appendix, Table S7), which is significantly more than expected by chance (two-tailed binomial test, P = 7.1e−4). Significantly, a majority of these loci (7 of 12) have another mPing insertion within 10 kb. Of these seven loci, five have a nearby mPing in HEG4, including two pairs of high frequency excision loci located close to each other (2.8 kb between Chr1:36267659 and Chr1:36270511; 2.6 kb between Chr5:15210313 and Chr5:15213006) and one (Chr1:6806761) located 9.6 kb from another parental mPing (Chr1:6816415) with four independent excisions (Fig. 4D and SI Appendix, Table S7). Two of the seven loci have a nearby mPing (1.9 kb and 2.3 kb away) that is not a parental locus but rather present in 45% of the RILs (SI Appendix, Table S7). For the remaining 5 of 12 high excision mPing loci, the nearest mPing is from 22.3 to 335.8 kb (SI Appendix, Table S7). Taken together, these findings indicate that the loss of mPing among the RILs can be attributable to structural instability, not negative selection. In support of this claim is the finding that expression analysis of the protein-coding genes associated with the 12 high frequency excision loci identified only a single differentially expressed gene between the two parental accessions (LOC_Os01g07300) (SI Appendix, Table S8).

Assessing Ping Activity.

Genome Instability as the Burst Continues.

If the burst is not likely to end with Ping silencing, how will it end? Data from this study strongly suggest that increases in mPing copy number destabilize the genome by elevating the chances of potentially deleterious SV. TEs have long been known to generate SVs. The first TE discovered by McClintock, Dissociation (Ds), was initially detected as a site of chromosome breakage (29). These so-called “breaking” Ds elements are a structurally distinct minority of characterized elements, often composed of tightly linked Ds copies (30, 31). Double-Ds elements form frequently because of the preference of the Ac transposase to catalyze local transposition (32, 33). Similar local hops and chromosome breakage are associated with TEs from all three kingdoms, including the Drosophila P and the bacterial Tn10 elements (34, 35).

In contrast, Ping does not promote local transposition of mPing as new insertions from a single integrated donor site in transgenic Arabidopsis thaliana, yeast, and soybean were not linked to that site (16, 36, 37). In operational terms, this means that mPing bursts scatter copies randomly throughout the genome, with only a few near existing insertion sites. However, as the burst continues, the probability that new insertions will be linked to existing mPing elements increases. Whether linked mPing elements are more likely to generate SVs was investigated in the RIL population in two ways: 1) by quantifying the relative frequencies of parental mPing excisions, and 2) by identifying SVs in the RIL population and determining whether any had mPing at one or both break points. For the analysis of excisions, the RILs were scored for the frequency of excision among the 466 parental mPing loci. The fact that 7 of the 12 loci with five or more excisions had another mPing within 10 kb suggests that tightly linked elements are unusually active with regard to transposition (SI Appendix, Table S7).

Whether higher transposition activity also contributes to the formation of a subset of SVs was assessed by identifying deletions and small duplications among the RILs. As mentioned in Results, a limitation of this analysis is that the use of short sequence reads precludes the identification of inversions and many larger SVs. Overall, 16 SVs were detected, and 9 were shown to have mPing at one (4 of 9) or both (5 of 9) break points (Table 2). Furthermore, the positions of the break points—at mPing termini—implicate Ping transposase activity in SV formation. This situation is reminiscent of the occurrence of Ds breakage only if Ac is in the genome (29), and that aberrant transposition is the likely mechanism (38, 39).

RNA Extraction and qRT-PCR.

RILs were grown as described above. Seedlings of 3-wk-old plants were used to extract total RNAs using the RNeasy Plant Mini Kit (Qiagen). Total RNAs were treated with amplification-grade RNase-free DNase I (Qiagen) to remove contaminating DNAs. The treated DNA-free total RNAs were reverse transcribed into complementary DNAs (cDNAs) using SuperScript III first-strand supermix (Invitrogen). The resulting cDNAs were used as templates to perform qRT-PCR with iQ SYBR Green Supermix (Bio-Rad) on a CFX96 system (Bio-Rad). Transcript levels of ORF1 and TPase were normalized to rice actin. Primers for qRT-PCR are in SI Appendix, Supplementary Data 4.

Construction of Recombination Bin and Linkage Map.

A total of 105,900 SNPs between parental accessions HEG4 and NB were identified and filtered using the Genome Analysis Toolkit (GATK) (46) as described previously (9). High-quality SNPs were used to genotype the RILs and a hidden Markov model (HMM) algorithm was applied to construct recombination bins (22). Briefly, paired-end reads of each RIL were mapped to the NB reference genome (MSUv7) by Burrows–Wheeler Aligner (BWA) aligner v0.6.2 (21, 47). Read alignments were sorted, and PCR duplicates were removed using MarkDuplicate as implemented in PicardTools (broadinstitute.github.io/picard/). Potential SNP sites were identified by the GATK toolkit following the best practice workflow of GATK (46). Only SNP sites matching parental SNPs were considered as informative sites. A genotyping error rate of 1.4% was estimated by counting the number of incorrectly genotyped SNPs from independently sequenced parental samples, NB and HEG4. Of these 105,900 SNPs, 1,311 were misgenotyped as HEG4 or heterozygous SNPs in the NB sample (864 HEG4 and 447 heterozygous, respectively) whereas 1,656 were misgenotyped as NB or heterozygous SNPs in the HEG4 sample (333 NB and 1,323 heterozygous, respectively). The resulting RIL SNP genotypes were imported into the R environment to construct a recombination bin map using the HMM approach as described in the Maximum Parsimony of Recombination (MPR) package (22). The expected proportion of genotypes for each locus was 49.75:0.50:49.75, with an error rate of 1.4%. Recombinant bins were used as markers, and the missing data were imputed using an argmax method in the R/qtl package (48). A linkage map was built using the haldane function in the R/qtl package (48).

Identification of mPing and Ping Copies in the RILs.

mPing copies, including elements present in the NB reference genome and elements present in the RILs or HEG4, were identified by RelocaTE2 (25). Briefly, RelocaTE2 identified paired-end read pairs that match the ends of mPing (junction reads) or internal sequences of mPing (supporting reads) using BLAT (49). After trimming mPing sequences, the remaining sequences from the original read pairs were searched against MSUv7 by BWA v0.6.2 to identify potential mPing insertions either present in the reference genome or in the RILs. Raw results from RelocaTE2 were filtered by removing mPing insertions with junction reads from only one end or mPing insertions without 3 bp target site duplications (TSDs), which is a hallmark of mPing transposition. Heterozygous mPing insertions were characterized by analyzing BAM files of raw read pairs aligned to MSUv7 (50). These mPing insertions with reads aligned to both sides of the flanking sequences were characterized as heterozygous mPing insertions. Ping copy number in RILs was characterized as described for mPing. Read pairs from each mPing/Ping insertion site were analyzed to distinguish mPing and Ping sequence differences. All candidate Ping insertions were confirmed by PCR or Southern blot.

Determination of mPing Insertion Sites in Single-Ping RILs.

Two experimental approaches (transposon display and pooled DNA sequencing) were performed on multiple individuals to determine de novo mPing insertions. Because mPing transpositions are rare in single-Ping RILs, multiple individuals were pooled to maximize detection of new insertions catalyzed by each Ping.

Analysis of mPing Excision Events.

Parental mPing insertions segregating in the RILs were identified. First, 51 mPing insertions exclusively in NB or shared between HEG4 and NB were identified as parental mPing loci. These loci have allele frequencies of ~50% (35 to 55%) for 44 NB-specific mPing loci or ~90% (89 to 95%) for 7 shared mPing loci. Second, 415 nonreference mPing insertions were identified from four independent sequenced HEG4 libraries (Sequence Read Archive [SRA] accession numbers: SRR1619147, SRR833485, SRR833529, and SRR1619373). Determination of the allele frequency of the 415 mPing loci in the RILs revealed that they were present in more than 10% of the 272 RILs (ranging from 16 to 74%), suggesting these mPing loci were segregating in the population. Together, these 466 mPing loci were used as parental mPing loci in the analysis.

The recombination bin map and SNPs that distinguish the parental genotypes were used to phase the haplotypes of mPing loci in each RIL. Empty mPing loci within an NB or HEG4 haplotype were considered as potential excision sites and further confirmed by analyzing the alignments of reads covering the mPing insertion sites. Excisions were characterized only if there were two or more reads covering the 20 bp flanking an mPing insertion. Footprints of excision events were characterized by small insertion/deletions (indels) in the Compact Idiosyncratic Gapped Alignment Report (CIGAR) alignments. Two approaches were applied to identify independent excision events: 1) a standard method that counts each excision as an independent event assuming that all RILs are independent; 2) a conservative method that counts independent excisions from each locus according to the configuration of footprints in each RIL. Excisions that have the same footprint were counted as dependent events. A binomial test was used to test for the significance of high excision frequency mPing insertions based on a frequency of 0.0051 excision events per locus per plant. Two methods identified similar sets of mPing loci with high frequency excision events (SI Appendix, Table S13). Excisions of high excision frequency mPing loci identified by both the standard method and conservative method were manually inspected and confirmed by a PCR approach.

Structural Variation in the RILs and Proximity to mPing Elements.

A bioinformatics approach was employed to detect potential sequence rearrangements in the RILs. Only deletions could be identified with confidence because of the limitations associated with the resolution of short reads (75 bp and 100 bp in this study) and small insertion size libraries (~200 bp in this study). CNVnator was used to analyze the read depth of bam alignment files in 272 RILs (28). A series of filters were used to obtain high-confident deletions from the raw results of CNVnator: 1) considering only those regions with normalized read depth less than 0.05 as candidate deletions; 2) excluding centromeres, telomere chromosome ends, and low coverage regions in NB and HEG4 genomes; and 3) any candidate regions should have at least one read supporting the junction of at least one end of deletions. Candidate SVs were further analyzed by a PCR approach or sequencing with additional primers (SI Appendix, Supplementary Data 4).

Differential Gene Expression Analysis.

RNAseq reads from NB and HEG4 were obtained from SRA BioProject PRJNA264731. Each sample has three biological replicates with 51 nucleotide (nt) single-end reads. Reads were aligned to MSUv7 with tophat v2.1.1 (52) using default parameters. Read counts of protein-coding genes were estimated by the htseq-count tool from HTSeq v0.7.2 (53). Differential gene expression analysis was performed with edgeR v2.99.8 (54) using the generalized linear model (GLM). Genes with a false discovery rate (FDR) ≤ 0.05 were reported as differentially expressed genes.

We thank Dr. Anna M. McClung (US Department of Agriculture) for maintaining the RILs; and Dr. Venkateswari J. Chetty for preparation of plant materials and sequencing libraries. This work was supported by NSF Grant IOS-1027542 (to S.R.W. and J.E.S.).