The account of the degree of damage to plants (in percentages) was carried out on the 14th day after inoculation, in accordance with the express method for assessing rice varietal resistance to blast. The assessment was carried out by taking two indicators into account: the type of reaction (in points) and the degree of damage (in percentages), using the ten-point scale of the International Rice Research Institute [12]:

The intensity of disease development (IDD, %) was calculated by the formula (Equation (1)):

$$\text{ADD} = \sum \text{(a} \times \text{b)} / \text{n} \times 9 \tag{1}$$

where IDD is the intensity of disease development (%), (a × b) is the sum of the products of the number of infected plants multiplied by the corresponding damage point, and n is the number of recorded plants (pcs).

Depending on the damage points, all varieties wee conventionally divided into 4 groups: resistant, intermediate, susceptible, and strongly susceptible.

#### **3. Results and Discussion**

The development of blast-resistant varieties and their rapid introduction into production is the most promising solution in the fight against this disease. However, the development of resistant varieties is one of the most difficult areas of breeding. The causative agents of the disease have a great potential for variability, which, combined with its colossal reproduction capabilities, provides the pathogen with the highest adaptive capabilities. Combining several effective genes of resistance on a genetic basis of the best varieties widely used in production is an effective breeding strategy for long-term resistance to variable fungal pathogens.

Based on the use of DNA marker breeding (marker-assisted selection (MAS)—breeding with use of DNA markers towards genes of interest), we introduced 5 blast resistance genes into domestic rice varieties adapted to the agro-climatic conditions of rice cultivation in southern Russia.

A series of crosses made it possible to obtain rice lines based on the varieties Boyarin, Flagman and Virage with the introgressed and pyramided blast resistance genes *Pi-1, Pi-2, Pi-33, Pi-ta*, and *Pi-b* in a homozygous state. During all cycles of backcrossing, the transfer of the dominant alleles of each such gene in the offspring was controlled by closely linked molecular markers. Plants with no resistance alleles in the genotype were discarded.

At the first stage of work in 2005 at Agrarian research center "Donskoy", 6 hybrids were obtained from crossing the varieties Boyarin and Virage with three donors of blast resistance carrying the *Pi-l, Pi-2*, and *Pi-33* genes. After analysis at the Federal Scientific Rice Centre, homozygous forms were identified for the dominant alleles.

At the second stage of work (2008), after crossing the *Pi-1* + *33* × Boyarin and *Pi-2* × Boyarin hybrids between themselves, it was possible to obtain forms with three pyramided genes simultaneously: *Pi-1, Pi-2*, and *Pi-33* in a homozygous state.

At the third stage of work (2010), they were hybridized with varieties—donors of the *Pi-ta* and *Pi-b* genes—for combining 5 genes. There were two types of crosses: ((*Pi-1* + *2* + *33*) × *Pi-ta*) × *Pi-b* and *Pi-1* + *2* + *33* × (*Pi-ta* × *Pi-b*).

Leaves were selected from the best F2 hybrid plants for DNA analysis at All-Russian Research Institute of Agricultural Biotechnology and the Federal Scientific Rice Centre using one marker for each gene. Based on the analysis results, it was possible to isolate two rice samples that were homozygous for all five dominant alleles. Reanalysis of the leaves of these samples confirmed last year's results, i.e., homozygosity for the dominant alleles of all five loci.

Figures 1 and 2 show the panicles of two lines, 1225/13 and 1396/13, which show the presence of dominant alleles at five loci in the homozygous state: *Pi-1*, *Pi-2*, *Pi-33*, *Pi-b*, and *Pi-ta*.

**Figure 1.** Panicle of the early-ripening line 1225/13.

**Figure 2.** Panicle of the mid-ripening line 1396/13.

Line 1225/13 is early maturing, matures in 110 days, and dwarfish (80 cm), with a small panicle (15 cm) (Figure 1).

The second line 1396/13 is mid-ripening, the period to maturity is 120 days, and it is taller (100 cm), with a large long panicle (22 cm) (Figure 2).

Against the infectious background in the Federal Scientific Rice Centre, the index of disease development (IDD) in this line was only 1.4%, while the variety Novator was damaged by 67.7%. The results of the analysis made it possible to send these lines to the breeding nursery in 2014–2015 for testing for yield and blast resistance. The variety **Pentagen** (1396/13), carrying 5 genes for blast resistance, is widely used in hybridization with high-yielding Russian varieties.

In the process of work at the Federal Scientific Rice Centre in 2007–2008, crosses were carried out and F1 hybrids were obtained from the combination (Flagman × C101-Lac) × (Flagman × C-101-A-51), which have the blast resistance genes *Pi-33* and *Pi-2* in their genotypes, respectively. The resulting F1 generation was used in backcrosses with the recipient parental forms. It should be noted that the F1 plants had a high degree of sterility (up to 95%). After the first series of backcrosses in 2008, the BC1 and BC2 generations were obtained in artificial climate chambers. In BC1 populations, fertility increased and averaged about 50%. Starting from the first backcrossing, marker control was carried out for the presence of transferred donor alleles in the hybrid offspring. In 2009, plants of the BC3 and BC4 generation were obtained. Among these plants, we selected the forms with the shortest growing

season and the highest panicle fertility. From the BC4F1 stage (the first self-pollination of rice plants, which makes it possible to transfer the donor allele to a homozygous state), individual selection was carried out. Segregation for Pi-2, Pi-33, and Sub1A genes fit into the Mendelian framework: in the second generation as a result of DNA analysis of the obtained plants, the ratio was 1:2:1 by genotype and 3:1 by phenotype.

Plants were selected that were closest in morphotype to the recipient parental form and had donor genes for resistance to the pathogen *Pyricularia oryzae* Cav. in their genotype in a homozygous state [22].

Figure 3 shows the results of PCR analysis for identification of the *Pi-33* blast resistance gene in the BC4F3 hybrid material.

**Figure 3.** Electrophoregram of genomic DNA amplification products at the loci RM310 and RM72: 1–4, 7–12, analyzed hybrid BC4F3 plants; 5, donor line of the Pi-33 gene C101-Lac; Flagman, maternal form.

The figure shows that plants Number 2, 4, and 7–12 are homozygotes for the dominant allele; plants Number 1 and 3 are heterozygous. The size of the PCR product in varieties with the dominant allele of the Pi-33 gene, which determines the resistance, is 198 bp; in varieties with a recessive allele, it is 152 bp.

In 2015–2016, the resulting rice lines with introgressed blast resistance genes *Pi-2* and *Pi-33* were crossed with the variety Khan Dan (Vietnamese breeding): the donor of the *Sub1A* gene. This work was performed for obtaining breeding material with combined genes for disease resistance and tolerance to prolonged immersion of plants under water. In 2017–2020, F4 and BC2F3 generations were obtained using climate chambers at the Federal Scientific Rice Centre (All-Russian Rice Research Institute, Krasnodar, Russia).

To increase economic efficiency and reduce labor costs, multi-primer systems have been developed to identify two genes (*Pi* and *Sub1A*) in a hybrid material simultaneously.

At the first stage, when we selected DNA markers for reliable interpretation of PCR products and identification of non-specifically amplified fragments, the following parameters were taken into account: the annealing temperature of specific pairs introduced into the reaction mixture, the difference in the size of PCR products synthesized during amplification with specific primer pairs (at least 100 base pairs), and the self-complementarity of the primer sequences.

The results of testing the combination of primer pairs flanking the marker regions of the *Pi-2* + *Sub1A* genes are shown in Figure 4.

**Figure 4.** Electrophoregram of multiplex PCR of genomic DNA amplification products at the loci RM527 and SSR140 for the Pi-2 gene and at the Sub1A203 locus for the *Sub1A* gene: 1–5, 9–13, analyzed hybrid plants of the BC2F3 generation; 6, Khan Dan, donor of the *Sub1A* gene; 7, Flagman, maternal form; 8, C101Lac-A-51, donor line of the *Pi-2* gene.

The electrophoregram shows that when PCR with such a combination of molecular markers is carried out, the target products specific for DNA markers of the desired genes are reliably amplified. Samples Number 3 and 12 have dominant alleles of the genes *Pi-2* and *Sub1A* in a homozygous state in their genotype; Samples 1, 4, 5, and 9 are homozygous for the *Sub1A* gene and have the *Pi-2* gene in the genotype in a heterozygous state; Sample 10 is a recessive homozygote for two target genes and was rejected. The size of the PCR product in varieties with the dominant allele of the *Pi-2* gene, which determines the resistance, is 233 bp. The size of the PCR product in varieties with the dominant allele of the *Sub1A* gene, which determines the resistance, is 118 bp. Clear identification on the electrophoregram makes it possible to reliably identify the presence of dominant alleles of the target genes.

The introduction of such varieties into production will allow us to avoid epiphytic development of the disease, preserving the biological rice productivity, and obtaining environmentally friendly agricultural products.

**Magnat** is the first cultivar in Russia created at the Agrarian Research Center Donskoy together with the Federal Scientific Rice Centre by the method of marker selection from a hybrid population (C101A-51 × Boyarin) × (C101 LAC × Boyarin) with transfer of blast resistance genes. Sample C101 LAC is a donor of the genes *Pi-1* and *Pi-33*, and C101A-51 is a *Pi-2*. The growing season is 125 ± 1 days and the plant height is 96 ± 2 cm. The panicle is erect and compact, 17.5 ± 0.5 cm long, and bears 185 ± 5 spikelets. The grain is oval, 8.3 ± 0.2 mm long, 3.1 ± 0.1 mm wide, and 2.2 ± 0.1 mm thick and weighs 24.0 ± 2.0 mg. The yield of the Magnat variety was 8.25 t/ha, which is 1.1 t/ha higher than that of the Boyarin standard.

The rice variety **Pirouette** was bred at the Agrarian Research Center Donskoy, together with the Federal Scientific Rice Centre, by the method of stepwise hybridization and marker breeding from a hybrid population (C101-A-51 (*Pi-2*) × Boyarin) × (C101-Lac (*Pi-1* + *33*) × Virazh). It contains three blast resistance genes: *Pi-1, Pi-2,* and *Pi-33.* The variety is mid-ripening, the growing season from flooding to full ripeness is 124 ± 1 days. The average yield of the variety Pirouette was 9.57 t/ha, which is 1.13 t/ha higher than that of the standard variety Yuzhanin. Plant height is 88 ± 2 cm; the panicle is erect, compact, and 17.5 ± 0.5 cm long and carries 165 ± 5 spikelets. The spikelets are oval, 8.9 ± 0.2 mm long, and 3.7 ± 0.1 mm wide. The weight of 1000 grains is 31.6 ± 2.0 g. The variety is resistant to lodging and shedding, is cold-tolerant, and germinates well from under a layer of water. It has been ncluded in the Register of Breeding Achievements of the Russian Federation for the North Caucasus region since 2020.

The rice variety **Kapitan** was bred at the Agrarian Research Center Donskoy in cooperation with the Federal Scientific Rice Centre by the method of triple backcrossing and marker breeding from the Flagman × IR-36 hybrid population. The variety is mid-ripening and the growing season from the flooding to full ripeness is 120 ± 1 days. On average, over the years of competitive testing, the yield of the variety Kapitan was 8.13 t/ha, which is 0.64 t/ha higher than that of the variety Yuzhanin. A higher yield of this variety is formed due to more grain in the panicle and an increased weight of the caryopsis. The average height of plants is 90 ± 2 cm; the panicle is semi-inclined, compact, and 18.5 ± 0.5 cm long; and the average number of spikelets is 140 ± 10 pieces (Figure 5). The grains are oval, 9.5 ± 0.2 mm long, and 3.6 ± 0.1 mm wide. The average weight of 1000 grains is 35.0 ± 2.1 g. The variety carries the *Pi-ta* gene and is resistant to blast, lodging, and shedding. The variety has been under state testing since 2019.

**Figure 5.** Rice panicles of the variety Kapitan.

The rice variety **Argamak** was bred at the Agrarian Research Center Donskoy by individual multiple selection of plants with the largest panicles from a hybrid population Il. 14 (*Pi-1, Pi-2, Pi-33*) × Kuboyar. The variety belongs to the mid-ripening group, and the growing season from flooding to full ripeness is 119 days. On average, over the years of competitive testing (2017–2019), the yield of the variety was 8.79 t/ha, which is 1.59 t/ha higher than that of the variety Yuzhanin. The maximum yield was formed in 2019: 10.1 t/ha, 2.55 more than the standard. The high yield of this variety was formed due to the greater grain content of the panicle than that of the standard and the increased density of the stem. Plant height is 93 ± 2 cm on average; the panicle is erect, compact, and 16 ± 0.5 cm long; the average number of spikelets is 142 ± 6 pieces. The grains are oval,

8.4 ± 0.2 mm long, 3.3 ± 0.1 mm wide. Weight of 1000 grains—31.1 ± 1.9 g. The variety is resistant to blast, lodging, and shedding. It has been tested at state varietal testing since 2020.

The rice variety **Lenaris** (Federal Scientific Rice Centre) had shown high adaptability, non-lodging, and the possibility for straight combine harvesting. Its yield was 10.6 t/ha. Plants had high spikelet fertility and short stems (77 ± 5 cm) and were resistant to the Krasnodar population of *P. oryzae* as well. Their panicle is slightly drooping and compact; its length is 18 ± 1.0 cm. The mass of 1000 grains is about 30.4 ± 1.8 g.

In 2013–2014, the Agrarian Research Center Donskoy conducted crosses and obtained F1–F2 hybrids of the variety Novator with Asian donor rice varieties carrying the *SalTol* and *Sub1A* genes. The hybrids of the second generation varied significantly in terms of quantitative traits: growing season (from early ripening to non-flowering), plant height (75–122 cm), panicle length (14–25 cm), number of filled grains (80–206 pcs), number of spikelets (99–300 pcs), panicle density (4.4–16.6 pcs/cm), 1000-grain weight (26.3–34.9 g), grain weight per panicle (2.1–5.5 g), etc.

Hybridization of the salt-sensitive domestic variety Novator with the lines IR52713-2B-8-2B-1-2, IR74099-3R-3-3, and IR61920-3B-22-2-1 (NSIC Rc 106)—*SalTol* locus donors—was carried out. The first generation of hybrids was used to generate an F2 hybrid population. From the populations of plants of the second generation, 90 early-ripening samples with well-ripened grains (30 in each combination of crossing) were selected, which were analyzed by PCR for the presence of introduced *SalTol* alleles. As an example, Figure 2 shows the data of electrophoretic analysis of PCR products with the Rm493 marker. The donor allele of the parental line NSIC Rc 106, designated as 2.2, was found in a homozygous state in Sample 282. The rest of the plants, whose amplification spectra are presented in this electrophoregram, carried the alleles of the donor and the variety Novator; that is, they were heterozygous for the *SalTol* locus (Figure 2). Similar results were obtained during DNA analysis of the studied rice samples with the RM7075 marker (Figures 6 and 7).

**Figure 6.** Electrophoregram of genomic DNA amplification products with RM 493: 1.1, Novator; 1.2, NSIC Rc 106; 17–296, hybrid plants NSIC Rc 106 × Novator; DNA marker (100–1500 bp).

**Figure 7.** Electrophoregram of genomic DNA amplification products with RM 7075: 1.1, Novator; 2.1, NSIC Rc 106; 17–286, hybrid plants NSIC Rc 106 × Novator; DNA marker (100–1500 bp).

In general, according to the results of DNA analysis of F2 hybrids, 17 plants homozygous for the dominant allele of the *SalTol* locus were identified, 29 samples carried *SalTol* in a heterozygous state, and 44 plants showed only recessive alleles inherited from the variety Novator.

Segregation for *SalTol* genes did not fit into the Mendelian framework, since the sample was unrepresentative due to selection. Plants with recessive alleles of the gene prevailed, and the number of salt-tolerant dominant homozygotes was less than the expected number. This is due to the linkage of *SalTol* genes with genes unfavorable for plants in our conditions: photosensitivity, late maturity, spikelet shedding, and spinosity.

Testing plants under salinity in the early stages of development is a quick, common method based on simple criteria. It was shown that at the initial vegetation stage, the length of the root and shoots and seed germination are potential indicators of resistance to the effects of increased salt concentrations [18,19]. Evaluation of the potential salt tolerance of the studied rice hybrids and their parental forms revealed significant variations in salinity tolerance depending on the genotype. The greatest decrease in seed germination—52%—was found in the salt-sensitive variety Novator. The line NSIC Rc 106 and second-generation plants, which were homo- and heterozygous for the *SalTol* locus, showed the highest resistance by seed germination (germination decrease of less than 5%). The donor lines IR52713-2B-8-2B-1-2, and IR74099-3R-3-3 and hybrid combinations obtained on their basis with the *SalTol* gene in a homozygous state also showed high resistance for this trait.

The least suppression of growth indices, as well as in the case of seed germination, was recorded in the lines NSIC Rc 106, IR52713-2B-8-2B-1-2, IR74099-3R-3-3, and *SalTol* homozygous plants from the F2 generation; the greatest decrease in the length of shoots and roots under salt stress was shown in the variety Novator and in hybrid plants that did not inherit the *SalTol* locus according to molecular analysis data. Thus, DNA analysis made it possible to simplify the breeding scheme and obtain salt-tolerant F2 hybrids carrying the *SalTol* locus in a homozygous state. These results indicate that the developed codominant markers of the *SalTol* locus RM 493 and RM 7075 are an effective tool for marker-assisted selection of salt-tolerant forms based on domestic rice genotypes.

Rice samples with *SalTol* genes in 2018–2020 were studied in a control nursery and in competitive variety testing; productive forms were identified.

At the same time, in 2013, hybrids were obtained by crossing the variety Novator with donors of the *Sub1A* gene. The Asian varieties turned out to be late-ripening and photosensitive and did not flower under our conditions. Hybridization was carried out only with the help of artificial climate chambers. The first generation in 2013 was characterized by a high degree of sterility (90–95%) and brown color of the flowering scales during maturation, which indicates significant genetic differences between the parental forms. In the second generation in 2014, a very large spectrum of splitting was observed in terms of the growing season, plant height, panicle length and shape, number of spikelets, and spinosity (Table 2).

**Table 2.** Variations in the quantitative traits in F2 hybrids from crossing submergence-resistant samples with the variety Novator, 2014.

Note: The average value is indicated in brackets.

This wide range of variability is not observed in other crops. This is due to the genetic and ecological-geographical remoteness of the crossed forms. In each combination, about 400 plants were selected for morphometric and genetic analysis. Among the F2 hybrids, we managed to select the best plants according to many traits, combining early maturity, optimal plant height, good grain size in panicles, non-shattering, and fertility of spikelets (Table 3).

**Table 3.** Selected F2 hybrid plants from crossing submergence-resistant samples with Novator, 2014.

Note \*: (BxN), BR 11 × Novator; (CxN), CR-1009 × Novator; (IxN, Inbara 3 × Novator; (TxN), TDK-1 × Novator.

PCR analysis of leaves was carried out in 20 plants of each of the four hybrids, as a result of which, forms with the *Sub1A* flood resistance gene were isolated. The electrophoretic analysis of PCR products with the RM 7481 marker is shown in Figure 8. The donor allele of the parental line CR-1009 was homozygous in Samples 2, 3, 5, 9, 13 and 17. Plants 2, 4, 6–8, 10, 11, 16, 18, and 19 were heterozygous at the *Sub1A* locus; that is, they carried both the alleles of the donor and the recessive alleles inherited from the variety Novator. Thus, according to the results of PCR analysis with the RM 7481 marker, 14 homozygotes of F2 plants at the *Sub1A* locus were identified, 40 samples carried *Sub1A* in a heterozygous state, and 22 plants inherited only the recessive allele from the variety Novator.

Of the analyzed BR-11 × Novator hybrid plants, the Sub1A gene (in homo- and heterozygous state) was present in nine, i.e., in a ratio of 9:11, although with monohybrid segregation, it should have been 15:5. In the hybrid combination CR-1009 × Novator, F2 segregated in a ratio of 18:2, i.e., almost all of the selected plants had the *Sub1A* gene. In the hybrids Inbara-3 × Novator and TDK-1 × Novator, segregation took place in a ratio of 14:6 or approximately 3:1, i.e., close to Mendelian.

The deviations in segregation of the two combinations can be explained by the influence of selection and gene linkage. A total of 55 plants with the target gene in the homo- and heterozygous state were isolated from 80 plants of four hybrids. The selected samples with the *Sub1A* gene in 2015 were reproduced in the Federal State Unitary Enterprise "Proletarskoye"of the Rostov Region, where the best F3 plants were selected from them for DNA analysis.

In F3 plants, significant morphological and biological segregation continued. Significant variation was noted for the growing season, plant height, size of panicles and caryopses, fertility, grain shedding, etc. The best forms were selected from them and leaves were taken for DNA analysis. At the next stage of work, in 2016–2020, constant lines carrying the *Sub1A* gene in a homozygous state were selected and tested for yield and resistance to prolonged water flooding. As a result, rice varieties for herbicide-free technologies will be developed, vigorously overcoming a deep layer of water in the germination phase with minimal seed loss.

#### **4. Conclusions**

1. As a result of the studies carried out using molecular marking based on PCR in combination with traditional breeding, early-maturing rice lines with genes for resistance to salinization (*SalTol*) and to flooding (*Sub1A*), which are suitable for cultivation in the south of Russia, were isolated.

2. Rice lines have been developed, the genotype of which contains five effective blast resistance genes (*Pi-1, Pi-2, Pi-33, Pi-ta,* and *Pi-b*). The introduction of such varieties into production will allow us to avoid epiphytotic development of the disease, preserving the biological productivity of rice and obtaining environmentally friendly agricultural products.

3. Samples of the F4 and BC2F3 generations were obtained with combined blast resistance (*Pi*) and prolonged flooding tolerance (*Sub1A*) genes as a factor in the control of weeds in the homo- and heterozygous state, which was confirmed by the data of their DNA analysis. The testing of the obtained rice breeding resources for submergence tolerance under laboratory conditions made it possible to select tolerant rice forms that will be studied in the breeding process for a complex of agronomically valuable traits. Their use will reduce the use of chemical plant protection products against diseases and weeds, thereby increasing the ecological status of the rice-growing industry.

The research was carried out with the financial support of the Kuban Science Foundation in the framework of the scientific project № 20.1/1.

**Author Contributions:** Conceptualization, P.K.; methodology, E.D. and K.A.; validation, P.K.; formal analysis, E.D.; investigation, E.D., P.K., M.R., S.L. and E.K.; resources, E.D.; writing—original draft preparation, E.D., P.K. and K.A.; project administration, E.D.; funding acquisition, E.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Federal Scientific Rice Centre.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### **Nanopore RNA Sequencing Revealed Long Non-Coding and LTR Retrotransposon-Related RNAs Expressed at Early Stages of Triticale SEED Development**

**Ilya Kirov 1,2,\*, Maxim Dudnikov 1,2, Pavel Merkulov 1, Andrey Shingaliev 1, Murad Omarov 1,3, Elizaveta Kolganova 1, Alexandra Sigaeva 1, Gennady Karlov <sup>1</sup> and Alexander Soloviev <sup>1</sup>**

Received: 24 November 2020; Accepted: 15 December 2020; Published: 17 December 2020

**Abstract:** The intergenic space of plant genomes encodes many functionally important yet unexplored RNAs. The genomic loci encoding these RNAs are often considered "junk", DNA as they are frequently associated with repeat-rich regions of the genome. The latter makes the annotations of these loci and the assembly of the corresponding transcripts using short RNAseq reads particularly challenging. Here, using long-read Nanopore direct RNA sequencing, we aimed to identify these "junk" RNA molecules, including long non-coding RNAs (lncRNAs) and transposon-derived transcripts expressed during early stages (10 days post anthesis) of seed development of triticale (AABBRR, 2*n* = 6*x* = 42), an interspecific hybrid between wheat and rye. Altogether, we found 796 lncRNAs and 20 LTR retrotransposon-related transcripts (RTE-RNAs) expressed at this stage, with most of them being previously unannotated and located in the intergenic as well as intronic regions. Sequence analysis of the lncRNAs provide evidence for the frequent exonization of Class I (retrotransposons) and class II (DNA transposons) transposon sequences and suggest direct influence of "junk" DNA on the structure and origin of lncRNAs. We show that the expression patterns of lncRNAs and RTE-related transcripts have high stage specificity. In turn, almost half of the lncRNAs located in Genomes A and B have the highest expression levels at 10–30 days post anthesis in wheat. Detailed analysis of the protein-coding potential of the RTE-RNAs showed that 75% of them carry open reading frames (ORFs) for a diverse set of GAG proteins, the main component of virus-like particles of LTR retrotransposons. We further experimentally demonstrated that some RTE-RNAs originate from autonomous LTR retrotransposons with ongoing transposition activity during early stages of triticale seed development. Overall, our results provide a framework for further exploration of the newly discovered lncRNAs and RTE-RNAs in functional and genome-wide association studies in triticale and wheat. Our study also demonstrates that Nanopore direct RNA sequencing is an indispensable tool for the elucidation of lncRNA and retrotransposon transcripts.

**Keywords:** long non-coding RNAs; seed development; Nanopore sequencing; retrotransposons; triticale

#### **1.  Traditionally, lncRNAs are considered to be protein non-coding, although many of them carry small open reading frames and encode functional peptides in different plants [1]. A broad range of functionality has been described for lncRNAs in plants, including the miRNA sponge, protein scaffolding, and post-transcriptional regulation of target genes via antisense pairing followed by mRNA decoy [2]. Based on their genomic localization regarding other genes, lncRNAs are broadly grouped into different classes: (1) lincRNAs or long non-coding intergenic RNAs; (2) NAT lncRNAs or natural anti-sense lncRNAs; (3) intronic lncRNAs, located in the introns; and (4) sense lncRNAs [3]. LncRNAs have been identified in many plant species and their expression in response to various abiotic and biotic stresses has been investigated, although our knowledge about lncRNA functions is still very limited [4–9]. However, the catalogue of lncRNAs for some plants with complex polyploidy and repeat-rich genomes, including wheat (*Triticum aestivum* L.), remains mainly underexplored [10].

Intergenic space and introns are frequently sources of lncRNA origin, although these genomic regions are enriched by insertions or remnants of transposable elements (TEs) [11]. Corroborating this, lncRNA sequences are more often associated with TEs than protein-coding genes and the bulk of them possess TE-related sequences [11–16]. For example, up to 75% of human lncRNAs have at least one exon with sequences originating from TEs [12,17]. A similar trend was demonstrated in some plant species, including maize, where 65% of lncRNAs had similarities to TEs [18]. More intriguingly, the TE-derived sequences may trigger the origin of new lncRNAs, providing a positive feedback loop with the evolution of lncRNAs [11,19]. TE-derived lncRNAs can have important and conserved biological functions [20]. For example, the rice lncRNA MIKKI is derived from LTR retrotransposons and has been shown to sequestrate miR171 and prevent degradation of its targets, mRNAs of SCARECROW-Like (SCL) transcription factors, in roots [20]. It is important to note that TEs can become transcriptionally and transpositionally active under certain circumstances, including stressful conditions, and in some developmental stages (e.g., microsporogenesis) and tissues (e.g., developing endosperm) [21–23]. However, individual TE transcripts and their coding potential have only been studied episodically in plants. The association of lncRNAs with repeat sequences like TEs makes the annotation of many lncRNAs challenging because of the ambiguity in the mapping of repeat-derived short RNAseq reads. RNAseq reads with multiple mapping positions in the genome sequence are frequently discarded from further analysis. Although some tools have been developed so far to overcome these obstacles, most of the lncRNA identification pipelines still ignore ambiguously mapped reads [24,25]. Thus, long-read sequencing technologies provide a great opportunity for transcriptome exploration, including the identification of transcribed repeats (e.g., transposable elements) or repeat-related transcripts [26]. Moreover, Panda and Slotkin [26] showed that by using Nanopore long cDNA reads, it is even possible to trace the expression of individual TEs from multicopy families in *Arabidopsis* and maize. The application of long-read technologies to the exploration of the lncRNA repertoire in plants has been demonstrated for several plant species, including *Oryza sativa* L. ssp. *japonica* [27], *Populus simonii* [28], *P. qiongdaoensis* [29], poplar "Nanlin 895" [30], *Trifolium pratense* [31], *Cardamine violifolia* [32], and *Vigna angularis* [33]. These studies demonstrated that long-read sequencing can be used to obtain a comprehensive catalogue of lncRNAs.

Grain development is one of the most important and practically relevant biological processes. It involves massive biochemical, physiological, and transcriptomic changes [34–36]. Wheat grain development is divided into five stages: (i) undifferentiated embryo and cellularization of the endosperm (0–7 days post anthesis (dpa)); (ii) embryo differentiation the embryo with formation of

the main cell types (transfer cells, aleurone, starchy endosperm and the surrounding cells) (7–14 dpa); (iii) root and leaf primordia differentiation, full kernel development and the milk-ripe stage (14–21 dpa); (iv) further growth and differentiation of primary and lateral roots, and the dough stage of endosperm (21–31 dpa); (v) fully differentiated embryo (31–50 dpa) [34]. Studies of the transcriptome during wheat seed development have been extensively elucidated using RNAseq sequencing [35,37–43]. Global transcriptome analysis of developing seeds has shown that the expression of protein-coding genes is highly dynamic. Recently, Madhavan et al. (2020) used publicly available Illumina RNAseq reads and identified 443 lncRNAs expressed during the grain filling stage (14 and 30 dpa) [44]. It is currently unknown which types of lncRNA are expressed during other stages of wheat seed development.

Here, we used Nanopore long-read sequencing to discover lncRNAs, TE transcripts, and TE-related lncRNAs that are specifically expressed during the cell proliferation stage of seed development (10 dpa) in triticale (× *Triticosecale* Wittmack, AABBRR genome, 2*n* = 6*x* = 42) an interspecific hybrid between wheat and rye (*Secale cereale* L.). We identified 796 lncRNAs and 20 LTR retrotransposon-derived transcripts, with most of them being previously unannotated. The majority of the detected retrotransposon RNAs had a single intron, carried open reading frames (ORFs) encoding for a diverged set of GAG proteins, and were encoded by potentially autonomous and non-autonomous retrotransposons. The lncRNAs were also expressed during wheat seed development and had high stage specificity. Moreover, we found that lncRNA loci were biased toward frequent TE sequence exonization and were mainly located in the intergenic regions of A and B genomes of triticale. Our work explored the lncRNA landscape during the early stage of wheat and triticale seed development and provides a unique dataset for further functional studies of lncRNA and TEs, and their implications for seed development. Finally, the identified lncRNAs can be further incorporated into genome-wide association studies for marker-assisted improvement of the bread quality of modern triticale genotypes.

#### **2. Material and Methods**

#### *2.1. Plant Material and DNA Isolation*

For this study, the spring triticale line "L8665" obtained from the Department of Genetics, Russian State Agrarian University, was used. For DNA isolation, seeds of this line were germinated in the dark (room temperature) during 5–7 days and genomic DNA was isolated by the cetyltrimethylammonium bromide (CTAB) protocol [45].

#### *2.2. Sample Collection and RNA Isolation*

Plants of the spring triticale line "L8665" were grown in a greenhouse under natural light conditions. Developing seeds at 10 days post anthesis and flag leaves were dissected and placed into liquid nitrogen. RNA was isolated by the ExtractRNA kit (Evrogen, Moscow, Russia) following the manufacturer's instructions. The RNA concentration and integrity were estimated by Nanodrop (Nanodrop Technologies, Wilmington, CA, USA) and gel electrophoresis using an 1.2% agarose gel with ethidium bromide staining.

#### *2.3. RT-PCR*

For RT-PCR, RNA was treated by DNAse I (Qiagen, Hilden, Germany, Q-79254) following the manufacturer's instruction and used for cDNA synthesis. cDNA was synthesized using the MMLV RT kit (Evrogen, Moscow, Russia). Primers used for RT-PCR and the expected fragment lengths are listed in Table 1. The CDC (Cell division control protein, AAA-superfamily of ATPases; Ta54227) gene was used as a reference because of its high expression stability as shown by a previous study [46].

#### **Table 1.** Primers used for RT-PCR.

The PCR conditions were 94 ◦C for 1 min; 35 cycles of 94 ◦C for 1 min, 58 ◦C for 1 min, and 72 ◦C for 1 min; and a final elongation of 72 ◦C for 3 min.

#### *2.4. Nanopore Direct RNA Sequencing and Transcript Assembly*

Poly-A+ mRNA was purified from 100 μg of total RNA by the Dynabeads mRNA DIRECT Kit (ThermoFisher Scientific, Waltham, MA, USA) following the manufacturer's instructions. Final poly-A+ RNA concentration was measured using a Quantus Fluorometer (Promega Corporation, Madison, WI, USA) and checked by gel electrophoresis. For Nanopore sequencing, a library was prepared from 1 μg Poly-A+ using the nanopore Direct RNA Sequencing Kit SQK-RNA002 (Oxford Nanopore Technologies, Oxford, UK). Direct RNA Sequencing (DRS) was carried out using MinION and flow cell FLO-MIN106. Basecalling was performed by Guppy (Version 4.0.11).

For transcript assembly, sequences of the A and B genomes and unanchored scaffolds (Un) of *Triticum aestivum* ([48]) downloaded from the EnsemblPlants server (https://plants.ensembl.org/ Triticum\_aestivum/Info/Index) were combined with *Secale cereale* genome sequences [49] into a single fasta file. Reads were mapped to the obtained fasta file by minimap2 software [50] with the '-ax splice' argument. The obtained sam file was converted to a bam file by samtools [51] (samtools view -Sb) followed by sorting of the bam file by bamtools [52]. Transcript assembly was performed by StringTie2 [53] with the following arguments: –L –j 2 –f 0.05. The obtained gtf file was converted to gff format by the gffread tool [54]. The sorted bam and gff files were used for read mapping visualization by the locally installed JBrowse [55].

To obtain a high-confidence set of transcripts, we extracted transcript sequences using gtf from StringTie2 assembly by the gffread tool [54]. The reads were then mapped back to the transcripts using minimap2 software (settings: -ax map-ont), and the bam file with only primary alignments and mapping quality > 30 was obtained using samtools with –F 256 –q 30 –b parameters. This bam file was used to count the number of reads per transcripts and transcripts with > 5 DRS reads were selected.

#### *2.5. Long Non-Coding RNA Prediction*

To identify lncRNAs, all high-confidence transcripts with a length of >200 bp were selected using biopython [56]. Transcripts with open reading frame (ORF) lengths of >300 bp predicted by getorf were filtered out. Protein coding potential was calculated by three tools: LncFinder [57], PLEK [58], and CNCI [59]. LncFinder [57] was run in Rstudio Version 1.2.1335 (http://www.rstudio.com/) with R version 3.6.0. The following parameters were utilized: parallel.cores = 20, SS.features = TRUE, format = "DNA", frequencies.file = "wheat", svm.model = "wheat". PLEK [58] and CNCI [59] were run with

the default settings. Only transcripts identified by all three tools as non-coding and without similarity to Pfam domains were classified as lncRNAs.

To identify lncRNAs with exons overlapping the annotated TEs, we intersected CLARITE TE annotation (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC\_RefSeq\_Annotations/v1.0/) with the genomic coordinates of the assembled exons.

#### *2.6. Retrotransposon-Related Transcript Annotation*

LTR retrotransposons (RTEs) were predicted in the genome using LTRharvest 1.5.10 with the default parameters [60] and LTRdigest 1.5.10 [61] with the following parameters: -aaout yes -pptlen 10 30 -pbsoffset 0 3 -pdomevalcutoff 0.001. Hidden markov model (HMM) profiles of RTE domains were downloaded from the GyDB database [62]. The gff3 file from the LTRdigest analysis was treated by a custom Python script (https://github.com/Kirovez/LTR-RTE-analysis/blob/master/LtrDiParser\_v2.2.py) to extract sequences of LTR retrotransposons possessing similarity to any RTE domains including GAG, reverse transcriptase, RNAse H, aspartic protease, and integrase. To identify retrotransposon-related transcripts, the TEsorter tool (parameters: -eval 0.00001 -db gydb) was run with the set of all high-confidence transcripts. The transcripts with a similarity to RTE domains were manually checked in the locally installed JBrowse [55]. TEsorter data were also used for RTE classification. We also ran TEsorter with confident transcripts and the RExDB [63] database (-db rexdb) to identify the transcripts of Class II transposons but no transcripts were detected.

#### *2.7. GAG Protein Analysis*

To find GAG proteins, ORFs were predicted for RTE transcripts and Blastp with corresponding proteins was run followed. GAG proteins were aligned by MAFFT [64] with the standard parameters and a phylogenetic tree was constructed using iTOL [65]. The multiple alignment visualization was carried out in Jalview version <IP_ADDRESS> [66]. RNA-binding motifs (CX2CX4HX4C, where X is any amino acid) were identified by a custom Python script (https://github.com/Kirovez/LTR-RTE-analysis/blob/ master/RBM\_GAG\_screen.py).

#### *2.8. Gene Ontology Enrichment*

Gene ontology (GO) enrichment analyses was performed using ShinyGO v0.61 [67] (http:// bioinformatics.sdstate.edu/go/) with an false discovery rate (FDR) <0.01.

#### *2.9. Expression Analysis*

For RNA-Seq analysis of lncRNA and RTE transcript expression in different organs and development stages, publicly available data were used (Table 2). Reads were mapped on the de novo assembled transcriptome by Hisat2 [68] with the default options. The obtained files with alignments were used to calculate RPM (read per million reads) values for every transcript. For this purpose, we used Salmon v0.8.1 [69] and the quant command with the default parameters.

**Table 2.** Publicly available RNAseq data used in this study.

#### *2.10. Extrachromosomal Circular DNA Isolation*

Extrachromosomal circular DNA (eccDNA) was isolated and amplified according to the protocol of Lanciano et al. [72] with several modifications. Briefly, 5 μg of genomic DNA was treated by Plasmid-Safe ATP-Dependent DNAse (Epicenter, Madison, WI, USA) for 48 h according to the manufacturer's instructions. DNA precipitation was carried out by 0.1 volume 3 M sodium acetate and 2.5 volume absolute ethanol, followed by overnight incubation at −20 ◦C. After centrifugation, the eccDNA pellet was obtained and exposed to the rolling circle amplification (RCA) reaction by the Illustra TempliPhi 100 Amplification Kit (GE Healthcare, Chicago, IL, USA) for 65 h at 28 ◦C. Detection of the eccDNA of LTR retrotransposon TaeST2.45518.1 (named '*MIG*', location in wheat genome: 7B: 312,336,869 ... 312,341,902 (5 kb)) was performed by inverse PCR with specific primers: Forward: CACACCACTAGCAACCTCCA ; Reverse: TGCTTGTGACAAGATGGGCA. The PCR conditions were 94 ◦C for 1 min; 35 cycles of 94 ◦C for 1 min, 58 ◦C for 1 min, 72 ◦C for 1 min; and final elongation at 72 ◦C for 3 min.

#### *2.11. Statistics and Visualization*

Statistical analysis was done in Rstudio Version 1.2.1335 (http://www.rstudio.com/) with R version 3.6.0. Bar plots, density plots, and box plots were drawn by ggplot2 [73]. Heatmaps were constructed by the ComplexHeatmap [74] R package.

#### **3. Results**

#### *3.1. Direct Oxford Nanopore RNA Sequencing*

Total RNA was isolated from whole spikes of hexaploid triticale (AABBRR) cv. L8665 at 10 days post anthesis (dpa). This RNA was used for direct RNA sequencing by MinION (Oxford Nanopore). In total, 1,100,000 direct RNA sequencing (DRS) reads with N50~1.1 kb were obtained (Figure 1). To assemble the transcripts, reads were mapped to the genome sequence created artificially by combining the A and B genome sequences of the wheat chromosome-level assembly and rye draft genome contigs. Overall, 82,785 transcripts from 74,904 loci were assembled with 47,378 and 26,169 loci located in genomes A/B (AB lncRNAs) and R. A total of 36,490 transcripts had >5 mapped DRS reads, representing a set of transcripts with high confidence that was used for further analysis.

To estimate the triticale seed development phase used for Nanopore sequencing, we determined the expression values (reads per million mapped reads (RPM)) of the key genes involved in starch biosynthesis (expression started at Phase 1), the genes of storage proteins (high molecular weight glutenins and gliadins), and those of the *wbm* protein (Table 3), which are expressed during the grain filling stage. We identified the expression of starch biosynthesis genes, while no genes of storage proteins or the *wbm* protein were expressed. This suggests that the seeds used for Nanopore sequencing were in the early stages of seed development (before 14 dpa).

Thus, using direct Nanopore RNA sequencing, we assembled a high-confidence transcriptome of triticale seed at the early development stage and detected the expression of key genes known to be involved in the biological process (starch biosynthesis) occurring at this stage.

#### *3.2. Long Non-Coding RNA Prediction*

The assembled high-confidence set of transcripts was used for long non-coding RNA (lncRNA) prediction. The following criteria were applied to distinguish lncRNAs from protein-coding transcripts: (1) transcript length of >200 bp; (2) transcripts with an ORF length of <300 bp; (3) transcripts classified as non-coding by three tools including LncFinder [57], PLEK [58], and CNCI [59]; and (4) transcripts with no similarity to any Pfam domains. Using these criteria, we identified 796 triticale lncRNAs (Supplementary Files S1 and S2) encoded by 780 loci in Genomes A (167 lncRNAs), B (212 lncRNAs), and R (410 lncRNAs) and in unanchored wheat scaffolds (seven lncRNAs) (Figure 1A). Most of the lncRNAs had lengths of <1000 bp (Figure 1B). LncRNAs (386 transcripts) encoded by the loci of

Genomes A and B or unanchored wheat scaffolds (AB lncRNAs) were used for further analysis because of the significantly better annotation of these genomes compared with the R genome. Intersection of the genome position of the AB lncRNA loci with lncRNA and mRNA loci previously annotated in the A or B wheat genomes showed that 281 (73%) of the triticale AB lncRNAs are located in the previously unannotated genomic regions (Figure 1C). This number is significantly higher (Fisher's Exact Test, *<sup>p</sup>*-value <sup>&</sup>lt; 2.2 <sup>×</sup> <sup>10</sup><sup>−</sup>16) than that for non-lncRNA AB transcripts (10%, 2523), pointing to the underexplored nature of lncRNA loci. Moreover, only 13% (106) of the AB lncRNAs in our dataset were previously known.

**Figure 1.** Identification and classification of triticale long non-coding RNAs (lncRNAs). (**A**) Bar plot showing the number of triticale lncRNA loci located on different chromosomes of Genomes A, B, and R and unanchored wheat contigs (Un). Sc indicates loci mapped to *Secale cereale* contigs. (**B**) Density plot of the lncRNA length distribution. (**C**) Pie graph showing the portion of the AB lncRNAs located in the unannotated genomic regions. (**D**) The portion of genes encoding AB lncRNAs and other RNAs (not classified as lncRNAs) that have exons with transposable elements. Three stars (\*\*\*) indicate significant differences based on Fisher's Exact Test for Count Data, p-value < 0.001. (**E**) RT-PCR with specific primers on three genic lncRNAs of distinct types (**F**), and RNA isolated from developing grain (10 days post anthesis (dpa)) and flag leaves. CDC: reference gene (the cell division control protein). (**F**) Types of genic lncRNAs and the corresponding wheat genes.

LncRNAs are frequently associated with transposons. We found that 111 (29%) AB lncRNAs had exons that had transposon sequences with a length of >50 bp (Figure 1D). This was significantly more than that expected by chance (10%, Fisher's exact test, *<sup>p</sup>*-value <sup>&</sup>lt; 2.2 <sup>×</sup> <sup>10</sup><sup>−</sup>16). Of the TE-related lncRNAs, 61, 47, and 3 lncRNAs had exons with similarity to Class I, Class II, and unclassified TEs, respectively.

The AB lncRNAs were classified regarding the position of the annotated wheat protein-coding genes, resulting in 17 genic lncRNAs. Gene ontology analysis revealed that the genes overlapping with genic lncRNAs were significantly enriched (FDR < 0.01) in several Gene Ontology categories including "vesicle-mediated transport", "lipid modification", "ATPase activity", and "hydrolase activity". The genic lncRNAs were of different types, with Intronic-antisense (two lncRNAs, Figure 1F (top)), exonic-anti-sense (five lncRNAs, Figure 1F (middle)), and exonic-sense (nine lncRNAs, Figure 1F (bottom)) transcripts being the most common types. In addition, exonic (sense and anti-sense) lncRNAs were found. The expression levels of three types of genic lncRNA (depicted in Figure 1F) belonged to the distinct types estimated in developing grain (10 dpa) and flag leaves (Figure 1E). The RT-PCR results suggested that two lncRNAs (lnc001 and lnc003) were expressed in both samples, while the expression of one type of lncRNA (lnc002) only occurred in developing grain (10 dpa).

Altogether, we identified hundreds of previously unknown genic and intergenic lncRNAs of triticale and showed that they frequently possess the remnants of Class I and Class II transposable elements.

#### *3.3. AB lncRNAs Are Prone to Tissue-Specific Expression*

To find lncRNAs with possible specific roles during seed development, we used wheat RNAseq data to estimate lncRNA expression in several seed developmental stages, leaf tissues, and pistils of wheat.

A non-zero expression value in at least one condition was obtained for 351 AB lncRNAs. To estimate any biases in AB lncRNA expression compared with all high-confidence transcripts assembled from DRS reads, we calculated the tissue-specificity index (TSI). The results showed that the TSI was significantly (according to the Wilcoxon rank sum test with continuity correction, *<sup>p</sup>*-value <sup>&</sup>lt; 2.2 <sup>×</sup> <sup>10</sup><sup>−</sup>16) higher for lncRNAs, suggesting high tissue specificity of lncRNA expression (Figure 2A). We further found that 46% of the lncRNAs had their highest level of expression at 10, 20 or 30 dpa, with 107 lncRNAs having their maximum expression level at 10 dpa; this is in accordance with the type of our triticale sample (10 dpa) used for RNA isolation (Figure 2B). Moreover, we identified 95 AB lncRNAs for which >90% of the sum of RPKM (reads per kilobases per million reads) values across all the samples accounted for the 10, 20, and 30 DPA stages (Figure 2C).

Thus, the expression pattern of the identified lncRNAs was found to be tissue-specific, with almost half of the lncRNAs demonstrating the maximum expression level during seed development.

#### *3.4. Retrotransposon-Related Transcripts Encoding GAG Proteins Are Expressed during Early Seed Development in Triticale*

Early embryonic and endosperm development may be accompanied by the activation of transposable element (TE) activity [22]. Therefore, we analyzed the assembled transcripts for the presence of open reading frames encoding TE-related proteins. We focused on the transcripts of Genomes A/B because of the high quality of wheat genome assembly and annotation compared with the rye genome. No transcripts corresponding to DNA transposons were detected. However, we found 20 transcripts (RTE-RNAs) carrying a single ORF with similarities to distinct proteins of LTR retrotransposons. Surprisingly, we found no transcripts encoding for the full set of RTE proteins (GAG and POL). To check whether any RTE-RNAs were encoded by LTR retrotransposons with detectable LTR sequences (RTEs), we predicted RTEs in Genomes A and B. The results showed that 5 and 10 RTE-RNAs (15) were transcribed from full-length (potentially autonomous) or non-complete (one or more RTE domains were not detected while both LTRs were present) RTEs, respectively. Thus, almost 25% of the RTE-RNAs were found to be transcribed from potentially autonomous RTEs. For five RTE-RNAs, no associated RTEs were predicted, but 75% (15) of the RTE-RNAs were found to carry an ORF encoding a single GAG protein (GAG-RNAs). Of those, eight (53%) and five (33%) GAG-RNAs were encoded by full-length (Figure 3A) and non-complete copies of LTR retrotransposons, respectively. For two (14%) GAG-RNAs, no corresponding RTEs were identified (Table 4). It should be noted that four GAG-RNA genes were found to be located in the introns of three annotated protein-coding genes in the sense or anti-sense orientation, including TraesCS2B02G261900 (sense and anti-sense), TraesCS5A02G298800 (anti-sense), and TraesCS1B02G222500 (sense). Two GAG-RNA genes (TaeST2.11597.1 and TaeST2.11598.1) were found to be located in the introns of the same gene (TraesCS2B02G261900).

**Figure 2.** Expression patterns of triticale AB lncRNAs based on wheat RNAseq data. (**A**) Boxplot of the tissue-specificity index (TSI) for lncRNAs and all AB high-confidence transcripts assembled from direct RNA sequencing (DRS) reads (TSI values close to 1 represent high tissue specificity). Stars indicate significant differences estimated by the Wilcoxon rank sum test with continuity correction (*p*-value < 2.2 <sup>×</sup> <sup>10</sup><sup>−</sup>16). (**B**) Bar plot showing the number of AB lncRNAs with the maximum expression (reads per kilobases per million reads, RPKM) level at a specific stage of wheat development (10, 20, and 30 days post anthesis (dpa)) or tissue (leaves and pistils). (**C**) Heatmap of expression values (log(RPKM)) of AB lncRNAs with >90% of RPKM values based on wheat RNAseq data (10, 20, and 30 dpa; pistils and leaves).

**Figure 3.** Identification, phylogenetic analysis, and expression pattern of RTE-related transcripts. (**A**) Full-length LTR retrotransposons of Tork (TaeST2.19707.1) and Retrofit (TaeST2.45518.1) lineages expressing a short isoform encoding the GAG protein (shGAG). The dark blue rectangles show ORFs encoding all retrotransposon proteins (genome scheme) and the GAG protein (isoform scheme). Orange color indicates untranslated regions (UTRs). The Nanopore direct RNA read alignment on the LTR retrotransposon genome sequence is also shown. (**B**) Neighbor-joining phylogenetic tree built from 15 GAG proteins. Red and gray highlight the GAG proteins expressed from RTEs without one or more encoded proteins and from loci without predicted RTEs, respectively. Red stars indicate GAG proteins with <300 aa. The vertical red line indicates a group of GAGs encoded by truncated retrotransposons possessing only GAG ORFs. (**C**) Multiple alignment of Ty1/Copia GAG proteins. Green and red highlight variable regions between the Tork and Retrofit regions. Gray shows the RNA-binding motif (CX2CX4HX4C, where X is any amino acid). (**D**) Heatmap of the log(RPKM) expression of isoforms encoding GAG proteins in wheat leaves and pistils, and during seed development.

**Table 4.** Number of retrotransposon-related transcripts (RTE-RNAs) in different groups based on completeness of associated LTR retrotransposons (RTEs) and the type of predicted RTE proteins encoded by open reading frames (ORFs) (reverse transcriptase, RT; aspartic protease, AP; RNAse H).

We further focused on GAG-encoding RTE-RNAs (GAG-RNAs) as the most represented group. Oxford Nanopore RNA sequencing provides a unique opportunity to analyze the exon–intron structure of GAG-RNA-encoded loci and predict the deduced full GAG protein sequences. Classification of the deduced GAG proteins showed that 85% (14) and 15% (1) of them belong to Ty1/Copia and Ty3/Gypsy elements, respectively. Based on the information from the transcript assembly, we grouped all GAG-RNAs into three categories based on the number of introns they possessed: (1) a single intron, (2) two introns, and (3) no introns. The vast majority (86%, 13) of GAG-RNAs were found to carry a single intron (Figure 3A), while one GAG-RNA from Ty3/Gypsy had no introns, showing that the exon–intron structure may differ between the two LTR retrotransposon superfamilies. To understand the functional role of splicing in the generation of GAG-RNA transcripts, we also predicted ORFs for unspliced RNA variants corresponding to the regions between two LTRs. We observed that unspliced transcripts expressed from the whole RTEs had significantly longer ORFs, resulting in GAG proteins being fused with other RTE proteins. Additionally, for two GAG-RNA encoding RTEs (TaeST2.19707.1 and TaeST2.45518.1, Figure 3A), only one very long ORF (>4000 bp) was predicted, while other RTEs were found to have two or three ORFs encoding distinct proteins. Thus, we showed that the splicing of the GAG-RNA isoform is critical to ensure the production of the entire GAG protein.

We next performed a comparison and phylogenetic analysis of the amino acid sequences of the 15 GAG proteins. The multiple alignment revealed significant differences between one Ty3/Gypsy GAG (TaeST2.45979, GAG length 514 aa) and the Ty1/Copia GAG proteins. We compared the Ty1/Copia GAG sequences and found that Ty1/Copia GAGs originated from two RTE lineages, Tork (seven sequences) and Retrofit (seven sequences). Phylogenetic tree analysis revealed three and two groups of highly similar (up to 99%) GAG proteins in the Retrofit and Tork lineages, respectively (Figure 3B). In addition, we detected pronounced sequence divergence between the GAG proteins of Tork and Retrofit lineages. The GAG protein of the Retrofit lineage has a ~50 aa-specific C-end which is not found in the Tork GAG sequences (Figure 3C). In turn, Tork GAGs have specific amino acid sequences before the RNA-binding motif (RBM) site. These differences are also reflected in the phylogenetic tree's topology, where the branches corresponding to the Tork and Retrofit Ty1/Copia lineages are readily distinguishable (Figure 3B).

In addition, we noticed that the divergence of groups in a single lineage correlated well with the completeness of the RTEs expressing corresponding GAG-RNAs; the most diverged group of GAG was the one with the Tork lineage and encoded by truncated retrotransposons (Figure 3B) that had a single ORF-encoding GAG.

We then analyzed the GAG protein sequences in more detail. In particular, we estimated the presence of the RNA-binding motif (RBM)(CX2CX4HX4C, where X is any amino acids), a specific part of GAG proteins, which is responsible for GAG–RNA interactions. We found that RBM could be identified in all except three GAG proteins, including one Ty3/Gypsy (TaeST2.45979.1) and two Ty1/Copia elements (TaeST2.14377.1, TaeST2.16660.1) of the Tork lineage. The Ty1/Copia GAG-RNAs without the CX2CX3GHX4X motif (TaeST2\_14377 and TAeST\_16660) are truncated GAG proteins with smaller protein lengths (175 aa and 182 aa vs >300 aa for the full-length Ty1/GAG protein) and probably very limited functionality. Altogether, the results of the phylogenetic analysis suggest that a divergent set of GAG proteins, including truncated GAGs with no RNA-binding capacity, is expressed during triticale seed development.

We then estimated the expression patterns of the GAG-RNAs in several seed developmental stages and in leaf tissues and pistils of wheat. Nine of the 15 GAG-RNA loci had a specific expression pattern with maximum expression levels at early developmental stages or in pistils. For two GAG-RNA loci (TaeST2.14377.1 and TaeST2.44075.1), the expression level was too low, indicating a possible triticale-specific expression pattern (Figure 3D). Thus, the expression data showed that most of the identified triticale GAG-RNAs were also expressed during wheat seed development, and some RTEs expressed genomic RNA (gRNA) as well as a short isoform (shGAG) carrying ORF for GAG protein.

Overall, our results suggested that tens of transcripts encoded by full-length and truncated LTR retrotransposon copies are expressed at early stages of triticale seed development. Three-fourths of these RNAs carry ORFs for encoding a set of GAG proteins of variable length and phylogenetic diversity.

#### *3.5. A Full-Length Ty1*/*Copia LTR Retrotransposon Is Active in Triticale Seeds*

To transpose in the genome, RTEs need to express the full-length genomic RNA (gRNA). Although we did not detect gRNA for the RTEs expressing TaeST2.19707.1 (RTE3B, location in wheat genome: 3B:555,156,557 ... 555,163,131 (6.58 kb)) and TaeST2.45518.1 (named '*MIG*', location in wheat genome: 7B:312,336,869 ... 312,341,902 (5 kb)) shGAG RNAs by Nanopore sequencing (Figure 3A), we performed RT-PCR with primer pairs designed to detect (a) shGAG isoforms and (b) gRNA isoforms. The RT-PCR analysis resulted in detection of the expression levels of shGAG and gRNAs in developing triticale seeds (10 dpa) and flag leaves, although the gRNA expression level was lower (Figure 4A). Next, we assessed whether the RTEs were capable of transposing. To answer this question, we determined the generation of the extrachromosomal circular DNA (eccDNA) by these RTEs using inverted PCR. EccDNAs are byproducts of RTE activity in plants [71]. We first determined that inverted PCR with the designed primers did not produce PCR products with genomic DNA. Unfortunately, the primers on RTE3B produced PCR fragments with genomic DNA; therefore, the activity of RTE3B could not be assessed by inverted PCR. We continued eccDNA detection only for *MIG* (Supplementary File S3). For this, we enriched the eccDNA fraction by exonuclease treatment of total genomic DNA using an enzyme that specifically cut linear DNA while leaving circular DNA molecules intact. The product was then amplified by rolling circle amplification, and inverted PCR was carried out. We enriched eccDNA in genomic DNA isolated from developing grain at 10 dpa, as well as glume and lemma tissue. The specific products were detected only for eccDNA isolated from developing grain, and no products were obtained with eccDNA of glume and lemma tissues. Thus, our results showed that RTE *MIG* expresses both shGAG RNA and gRNA isoforms and has transposition activity.

Here, we provide experimental evidence suggesting that some detected RTE-RNAs originate from autonomous LTR retrotransposons with ongoing transposition activity in triticale at early stages of seed development.

**Figure 4.** Expression and extrachromosomal circular DNA (eccDNA) formation. (**A**) RT-PCR detection of the shGAG (a short isoform carrying ORFs for GAG protein) and gRNA isoforms of RTE3B and *MIG* full-length RTEs. CDC: reference gene (the cell division control protein). (**B**) Inverted PCR with genomic DNA and eccDNA-enriched fractions obtained from developing seeds (10 dpa) and glume and lemma triticale tissues. The positions of the primers on eccDNA are shown in the small representation in the right. PCR with shGAG primers with total and eccDNA-enriched DNA was used as a control.

#### **4. Discussion**

#### *4.1. A Set of Intergenic lncRNAs Detected by Nanopore Sequencing Is Expressed during the Early Stage of Triticale Seed Development*

Wheat seed development is a dynamic multistage process that involves significant changes in the transcriptome landscape. Here, we uncovered the lncRNA transcriptome of triticale seed at 10 days post anthesis (dpa), corresponding to the second stage of grain development, known as embryo differentiation (7–14 dpa) [34]. This stage is characterized by the formation of Type A starch granules and the expression of genes involved in starch biogenesis [36]. In turn, the genes encoding for the main storage proteins (the high molecular weight (HMW) glutenins and gliadin [36]) and the genes involved in storage protein biogenesis (e.g., wbm [75]) are expressed during the mid-development stage (14–21 dpa). In agreement with this, we detected the expression levels of the starch metabolism genes (Table 3), while expression of the storage protein genes (HMW glutenins and gliadins) and the wbm gene (found recently in the analyzed line (L8665) [76]) was not detected by Nanopore read analysis. These results prove that the analyzed triticale transcriptome corresponds to the early stage of grain development.

Our analysis showed that hundreds (798) of lncRNAs were expressed during this stage. Surprisingly, we found that 87% of the A and B lncRNAs were expressed from as yet unannotated regions of Genomes A and B. LncRNAs are often underrepresented in plant genome annotation. For example, 8009 lncRNAs were previously identified in the intergenic regions of barley [77], and 1760

unannotated lncRNAs were identified in foxtail millet [78]. The annotation of the lncRNAs expressed from the intergenic space can be challenging because of the biological properties of lncRNAs, including their high tissue- and stage-specificity. Indeed, we identified a high tissue-specificity index for triticale lncRNAs, suggesting that a large number of lncRNAs are expressed during a narrow time window. Additionally, lncRNAs often possess exons containing transposon-similar multicopy sequences that make the transcript assembly difficult because of unambiguity in short RNAseq read mapping. Here, we found that almost 30% of the triticale lncRNAs possessed exons with similarity to TEs. This is in accordance with previous reports on rice [79], where 73% lncRNAs were found to overlap with different TEs. Furthermore, 9.18% of sunflower lncRNAs have TE-related exons [80]. In these terms, RNAseq-based lncRNA identification can underestimate the number of lncRNAs in a cell or lead to transcript misassembly, because short reads from the repeat regions are often mismapped or discarded from the analysis [25]. However, notably, the study of plant lncRNAs and transposon-derived transcripts has mostly been limited by short read RNAseq data. Here, to escape lncRNA identification biases because of the short length of the RNAseq reads, we applied Nanopore direct RNA sequencing. This approach allowed us to precisely determine the transcribed regions in a complex wheat genome. In addition, Nanopore direct RNA sequencing is strand-specific and can be used to identify natural anti-sense lncRNAs. The application of third-generation sequencing could help to illuminate "the dark side" of developing seed transcriptomes involving lncRNAs and transposon-derived loci, thereby overcoming the obstacles in lncRNA discovery by short RNAseq reads. It will be especially useful for crop plants where lncRNA loci can be included into genome-wide association studies (GWAS) to determine which of them can influence key agronomical traits. Notably, lncRNA loci have not been well used for GWAS analysis so far [81]. We believe that the triticale lncRNAs identified in this study could act as valuable new targets for marker-assisted selection.

#### *4.2. Transcripts Encoding Diverse GAG Proteins Are Expressed during the Early Stage of Seed Development*

Seed development is accompanied by epigenetic relaxation, which may trigger retrotransposon (RTE) activity [21–23]. Here, we identified 20 transcripts with similarities to RTE proteins. Interestingly, we found that 75% of the RTE-related transcripts carry ORFs encoding GAG proteins, the main component of the RTE virus-like particle. It is known that during their lifecycle, RTEs have to produce significantly more GAG proteins that other RTE proteins. To ensure that there are excess GAG proteins compared with other RTE proteins, an *Arabidopsis* RTE called *EVD* encodes a special isoform (shGAG) encoding the GAG protein [82]. Because almost no systematic studies of RTE transcripts have been carried out at the single isoform level in crop plants, it was not previously clear whether shGAG transcript production is a common pattern for plant species. In this study, we detected isoform production for several full-length triticale RTEs and showed their transcription in wheat, implying that shGAG is a very common pattern for plants. Moreover, we showed that one of the full-length RTEs, TaeST2.45518.1, is capable of producing extrachromosomal circular DNAs (eccDNA), which are a byproduct of RTE activity and have been used to isolate transpositionally active RTEs [74]. Whether this RTE is active in the triticale embryo and can produce copies that are transmitted to the next generation or whether it is active in other tissues of developing grains (e.g., endosperm) warrants further investigation.

Based on the knowledge of splicing patterns of GAG-encoding isoforms, we were able to predict the GAG protein sequence and analyze it in more detail. Our results point to the existence of a divergent group of GAG proteins expressing during triticale and wheat seed development. These GAG proteins and the RTEs encoding them have several distinct features: (1) the RTEs encoding these GAG proteins are non-autonomous elements and possess no similarity to POL proteins; (2) the lengths of most of these GAG proteins (169–270 aa) are less than that of the conventional GAG protein (>300 aa); and (3) half of these GAG proteins lack the RNA-binding motif and cannot interact with RTE RNAs. The elements encoding GAG-RNA consist of two LTRs and the internal part is similar to the GAG protein. This structure makes these elements very similar to the previously described TR-GAG elements (terminal repeat with GAG domain) ([83]) found in many angiosperm species. However, the TR-GAG elements described in the current paper are classified as Ty1/Copia, while GAG proteins of previously identified TR-GAG elements are similar to both Ty1/Copia and the Ty3/Gypsy superfamily. While further functional and evolutionary studies are required, we suggest that these GAG loci are intermediate products of RTE diversification or "domestication". The "domestication" of GAG proteins has been documented in animals and insects [84–87]. Because of the short lengths of these GAG proteins, it can be also suggested that they may be involved in control of the activity of functional RTEs via incorporation into their virus-like particles. The mechanism of copy number control of RTEs via virus-like particles (VLP) misassembly, which is caused by a truncated GAG form, known as dominant-negative factor, has been well described for yeast [88,89]. No such examples have been described in plants. In the future, it will be intriguing to check whether TR-GAG proteins are capable of interacting with normal GAG proteins of full-length RTEs, which, as we showed here, are expressed in the same developmental stage.

Previously, Nanopore long-read sequencing of transcriptomes was used to annotate expressed transposons in *Arabidopsis* mutants lacking key systems of TE suppression, resulting in elevated expression of transposons [26]. The authors detected the expression of almost 1300 TEs. However, elucidation of transposon expression a "wild-type" genetic background would provide a unique opportunity to trace natural evolutionary forces shaping plant retrotranscriptome and to more deeply understand the features of host–transposon interactions. Recent studies on maize and sunflower [26,90] and our current results point to the great advantage of Nanopore RNA sequencing to decipher RTE expression in crop plants, even those with complex genomes such as triticale. Together with the growing number of publicly available long-read RNA sequencing datasets, this opens a broad avenue for studies of transposon expression in plants on the isoform-based level.

#### **5. Conclusions**

Here, using Nanopore direct RNA sequencing, we identified hundreds of previously unknown lncRNAs and LTR retrotransposon-derived transcripts that are expressed in the early stages of triticale and wheat seed development. We showed that triticale lncRNAs often possess similar sequences to transposons and their expression has high stage and tissue specificity, with half of the lncRNAs having the highest expression level at 10–30 days post anthesis in wheat. In addition, we found that most of the detected retrotransposon-related RNAs have a single intron, carry ORFs encoding for a divergent set of GAG proteins, and are encoded by potentially autonomous and non-autonomous retrotransposons. Of these, we identified one Ty1/Copia LTR retrotransposon that produces extrachromosomal circular DNA, and we suggest that it has transposition activity in developing triticale seeds. Finally, this study identified a unique set of lncRNAs and LTR retrotransposons expressed in the early stages of seed development, which we believe will be useful for further exploration of their functional potential and the association with phenotypic variation in triticale and wheat.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2223-7747/9/12/1794/s1, Supplementary File S1: Genomic location and overlapped genes for 976 Triticale IncRNAs, Supplementary file S2. Fasta sequences of 796 lncRNAs, Supplementary file S3. Sequence of RTE7B LTR retrotransposon that is expressed and produce eccDNA in triticale developing seed at 10 dpa.

**Author Contributions:** Conceptualization, I.K. and A.S. (Alexander Soloviev); methodology, I.K.; formal analysis, I.K., M.D., P.M., A.S. (Andrey Shingaliev), M.O., E.K., and A.S. (Alexandra Sigaeva); writing—original draft preparation, I.K.; writing—review and editing, I.K., A.S., and G.K.; funding acquisition, A.S. and G.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Education and Science of Russian Federation (Goszadanie No 0431-2019-0005).

**Conflicts of Interest:** The authors declare no conflict of interest.

**Data Availability Statement:** Nanopore data produced for this study are available in Sequence Read Archive (SRA) NCBI under Bioproject Accession PRJNA683988.

#### **References**

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Factors Influencing Genomic Prediction Accuracies of Tropical Maize Resistance to Fall Armyworm and Weevils**

**Arfang Badji 1,\*, Lewis Machida 2, Daniel Bomet Kwemoi 3,\*, Frank Kumi 4, Dennis Okii 1, Natasha Mwila 1, Symphorien Agbahoungba 5, Angele Ibanda 1, Astere Bararyenya 1, Selma Ndapewa Nghituwamhata 1, Thomas Odong 1, Peter Wasswa 1, Michael Otim 3, Mildred Ochwo-Ssemakula 1, Herbert Talwana 1, Godfrey Asea 3, Samuel Kyamanywa <sup>1</sup> and Patrick Rubaihayo <sup>1</sup>**

**Abstract:** Genomic selection (GS) can accelerate variety improvement when training set (TS) size and its relationship with the breeding set (BS) are optimized for prediction accuracies (PAs) of genomic prediction (GP) models. Sixteen GP algorithms were run on phenotypic best linear unbiased predictors (BLUPs) and estimators (BLUEs) of resistance to both fall armyworm (FAW) and maize weevil (MW) in a tropical maize panel. For MW resistance, 37% of the panel was the TS, and the BS was the remainder, whilst for FAW, random-based training sets (RBTS) and pedigree-based training sets (PBTSs) were designed. PAs achieved with BLUPs varied from 0.66 to 0.82 for MW-resistance traits, and for FAW resistance, 0.694 to 0.714 for RBTS of 37%, and 0.843 to 0.844 for RBTS of 85%, and these were at least two-fold those from BLUEs. For PBTS, FAW resistance PAs were generally higher than those for RBTS, except for one dataset. GP models generally showed similar PAs across individual traits whilst the TS designation was determinant, since a positive correlation (R = 0.92\*\*\*) between TS size and PAs was observed for RBTS, and for the PBTS, it was negative (R = 0.44\*\*). This study pioneered the use of GS for maize resistance to insect pests in sub-Saharan Africa.

**Keywords:** prediction accuracy; mixed linear and Bayesian models; machine learning algorithms; training set size and composition; parametric and nonparametric models

#### **1. Introduction**

Insect damage on maize plants and stored grains potentially impedes food security in Africa [1–3]. The fall armyworm (FAW) and stem borers in the field, and the maize weevils (MWs) in storage facilities, are some of the most devastating insect pests on the continent. These insect pests cause yield losses ranging from 10–90% leading to loss of grain marketability, and consumer health concerns due to the possible contamination of the grain with mycotoxins, such as aflatoxins [3–6]. In Africa, tremendous efforts were made during the last two decades to build host plant resistance to insect pests in maize through traditional pedigree (phenotypic)-based selection (PS) with substantial desirable results. Several Africa-adapted maize lines were developed and successfully tested for resistance to MW damage on grains [7–12]. Some of the success stories are from the International Center for Maize and Wheat Improvement (CIMMYT) of Kenya through the Insect Resistant Maize

**Citation:** Badji, A.; Machida, L.; Kwemoi, D.B.; Kumi, F.; Okii, D.; Mwila, N.; Agbahoungba, S.; Ibanda, A.; Bararyenya, A.; Nghituwamhata, S.N.; Odong, T.; et al. Factors Influencing Genomic Prediction Accuracies of Tropical Maize Resistance to Fall Armyworm and Weevils. *Plants* **2021**, *10*, 29. https://dx.doi.org/10.3390/ plants10010029

Received: 13 July 2020 Accepted: 14 September 2020 Published: 24 December 2020

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/ licenses/by/4.0/).

for Africa project (IRMA) that produced several storage pest and stem borer resistant maize lines [7,8,13–15]. On the other hand, the FAW is a new pest on the continent, first reported in 2016 in West and Central African countries [16], from where it spread throughout the African continent [17]. Hence, although efforts to develop FAW resistant varieties are underway at several institutions, including CIMMYT, published reports of FAW resistant varieties are not yet available [18,19].

The complex nature of insect resistance traits makes PS slow and expensive, and thus, difficult to implement, especially for resource-constrained breeding programs [20]. Application of traditional marker-assisted selection (MAS) is hampered by the necessity to first discover resistance-associated genomic regions through genetic linkage and genome-wide association mapping methods, both with several shortcomings, especially for complex traits [21–23]. In addition, genetic linkage and genome-wide association mapping studies have seldom been explored in African germplasm [8,24], which further impedes the application of MAS in the development of insect resistance maize germplasm in Africa. In a previous study, we discovered several quantitative trait nucleotides and genes that are putatively associated with FAW and MW resistance, confirming the quantitative nature of these traits, hence the difficulty in improving these traits through MAS [25]. An alternative to both PS and MAS is genomic selection (GS), which uses whole-genome markers to perform genomic prediction (GP) of breeding values of unphenotyped genotypes, from which one can select superior candidate genotypes for crossing to produce hybrids or to advance to the next generation [26]. GS was reported to achieve up to threefold annual genetic gain in maize improvement when compared to MAS, due to a more efficient accounting of trait-associated quantitative trait loci (QTL), faster selection cycles, and lower phenotyping costs [27–33].

Several statistical and machine learning GP models with various strengths and weaknesses have been developed to adapt to different contexts that are partly influenced by the genetic architecture of traits (number and effect size of QTL, proportions of additive and non-additive genetic effects) and reproductive classes of plants (allogamous vs. autogamous vs. clonally propagated) [34–36]. Therefore, to effectively implement GS in crop improvement programs, it is necessary to employ a holistic approach to determine the best GP strategy for particular breeding targets for given crop species [31,37]. Statistical models for GS vary in their prior assumptions and treatment of marker effects [31]. Parametric models focus on parameter estimates rather than prediction, while nonparametric algorithms give priority to prediction and have fewer assumptions [38]. Some parametric methods assume the SNP effects follow a normal distribution with equal variance for all loci, which seems unrealistic in practice.

Representative parametric methods are ridge regression best linear unbiased predictors (RR-BLUP) [39] and genomic BLUP (GBLUP) [40]. GBLUP was the first GP method to be developed, and it replaced the traditional pedigree-based relationship matrix with a genomic information-based matrix to improve prediction accuracies (PAs) [41]. Parametric methods BayesA [26] and weighted Bayesian shrinkage regression (wBSR) [42], on the other hand, consider a prior distribution of effect with a higher probability of moderate to large effects. Regarding parametric models such as BayesB [43] and BayesCπ [44], assumptions are made that consider the effects of some SNPs to be zero. The Bayesian least absolute shrinkage and selection operator (Bayesian LASSO) assumes that the effects of all markers follow a double exponential distribution [45], whilst the Bayesian sparse linear mixed model (BSLMM), a parametric method developed by Zhou et al. [46], combines the hypotheses of both GBLUP and Bayesian methods and achieves higher PAs than BayesCπ and BayesLASSO. Nonparametric or semi-parametric approaches such as random forest and reproducing kernel Hilbert space (RKHS) [47,48] are better suited for accounting for non-additive genetic effects (37,38), in contrast with parametric genomic prediction models [23,38,47,49]. Several studies compared the performances of GP models under different conditions. In a simulation study, Meuwissen et al. [26] found that while GBLUP achieved PAs of up to 73.2%, BayesA and BayesB comparatively provided additional increases of

around 9% and 16%, respectively. However, when a population is composed of close relatives and the target traits are controlled by several small effect genes, the different methods perform similarly [50–52]. On the contrary, BayesB and BayesCπ are better when dealing with distant relatives and traits affected by a small number of large-effect loci [23]. Kernel methods such as RKHS are robust in predicting non-additive effects and in solving complex multi-environment multi-trait models [53,54]. Compared to the above-mentioned parametric methods, deep learning techniques such as support vector regression (SVR), multilayer perceptron, and convolutional neural networks models performed poorly in some studies [55,56]. However, there are also instances where RKHS outperformed one or several of the parametric methods, for instance, GBLUP, rrBLUP, and Bayesian algorithms, in terms of several traits in several crops including maize [51,57–59]. These results were most likely because nonparametric GS models capture more adequately the non-additive genetic components which are an essential characteristic of complex traits [23,37,38] and hence could be good candidate tools for the prediction of FAW and MW-resistance traits which are controlled by both additive and non-additive gene action [21,23,31,41,60]. Therefore, since GP for maize resistance to insect pests such as FAW and MW is not yet well explored, it is pivotal to compare performances of several available prediction algorithms to inform better future GS programs. Therefore, the Genomic Prediction 0.2.6 plugin of the KDCompute 1.5.2. beta (https://kdcompute.igss-africa.org/kdcompute/home), an online database developed by Diversity Array Technologies (DArT, https://www.diversityarrays.com) for the analysis of DArT marker data, presents great interest for this purpose. It hosts a suite of parametric, semiparametric, and nonparametric GP methods that can be run simultaneously on genotype-phenotype datasets.

Additional factors that influence PAs are the different sizes of the training sets (TSs) and breeding sets (BS) and their genetic relationships, the number of markers used to estimate genomic estimated breeding values (GEBV) of lines, the population structure, and the extent of linkage disequilibrium [21,23,31,41,60]. Since phenotyping is the current bottleneck in plant breeding and one of the disadvantages of GP is the requirement of large TSs for high PAs to be achieved, determination of effective TS composition and size is critical for effective implementation of GS in crop improvement programs [21,61–64]. Additionally, the best TS determination will depend on the genetic architecture and the extent of population structure of the trait targeted for GP [63], two parameters that are substantially variable among plant breeding traits. Another factor that is a determinant of the predictive ability is the kinship between the TS and the BS (63). Several methods are used for TS optimization and these generally fall into two categories—namely, untargeted and targeted approaches. For the untargeted approach, the TS is determined independently of its genomic information, whereas the targeted method considers the genomic relationship between the TS and the BS as a means of maximizing PAs [65]. However, deciding on the best TS selection method is not straightforward and depends on context [66].

Furthermore, in maize, GPs were previously conducted using either genotypic best linear unbiased estimators (BLUEs) [67–69] or best linear unbiased predictors (BLUPs) [31, 41,70] as means of phenotypic correction [70]. BLUEs are obtained by treating the genotypic effect of a mixed linear model as fixed effects and provide an estimated mean for each individual of a population equal to its true value. On the other hand, BLUPs generated by considering the genotypic factors as random and allowing for the shrinkage of the means towards the population mean [71]. Whether to use BLUPs or BLUEs in GP is debatable. Phenotypic BLUEs allow avoiding double penalization which BLUPs suffer from. With phenotypic BLUPs, this double penalization is, however, compensated through maximization of the correlation between predicted and true line values, while phenotypic BLUEs do not rely on this shrinkage [70]. However, the shrinkage in the BLUP procedure accounts better for outliers and environmental variabilities [72], permitting better estimates of individual genetic effects than BLUEs [71], and therefore, it usually yields more accurate predictions of phenotypic performance [70,72,73]. Furthermore, BLUPs are better in handling unbalanced data, wherein, for example, the number of individuals is not the

same in different locations or in the different replications of an experiment [49,70]. On that basis, the current study was conducted to evaluate the efficacies of different parametric, semiparametric, and nonparametric methods from both statistical and machine learning GP models in generating prediction accuracies (PAs) for maize resistance to FAW and MW in a diverse panel using both genotypic BLUEs and BLUPs.

#### **2. Material and Methods**

#### *2.1. Genetic Material and Field Experiments*

The panel used in this study consisted of 358 maize lines with diverse genetic and geographic backgrounds, and they were sourced from the National Crop Resources Research Institute (NaCRRI/Namulonge, Uganda), the International Institute for Tropical Agriculture (IITA/Ibadan, Nigeria), and The International Maize and Wheat Improvement Center (CIMMYT/Nairobi, Kenya). The panel consisted of 71 inbred lines developed for various purposes at NaCRRI; 28 and five stem borer (SB)-resistant inbred lines from CIM-MYT [6,13,14] and IITA, respectively; 19 storage pest (SP)-resistant inbred lines [7,8]; and a doubled haploid (DH) panel of 235 lines developed at CIMMYT using six parents—three of which were stem borer-resistant, one was a storage pest-resistant inbred line (these were also included in the population), and two were CML elite lines (one, CML132 was included in the panel) (Supplementary Materials Table S1).

The panel was planted and evaluated in three environments, at Mubuku Irrigation Experimental Station in Kasese, western Uganda in 2017 (316 lines) during the second rainy season (2017B) and the National Crop Resources Research Institute (NaCRRI), Namulonge, central Uganda in 2018 (92 lines) and 2019 (252 lines) both during the first rainy seasons (2018A and 2019A, respectively). Detailed information on these locations is presented in Table 1.

**Table 1.** Geographical, climatic, and soil characteristics of the planting locations [74].

Each combination of location and season was considered an environment, resulting in a total of three environments. An augmented experimental design was adopted in all three environments using six checks in 2017B, two in 2018A, and four in 2019A replicated in all the blocks. The experiments in 2017B, 2018A, and 2019A consisted of twelve, five, and ten blocks, respectively, containing the replicated checks and unreplicated lines and the experiment in 2018A was replicated twice.

#### *2.2. Genotyping, Quality Control, and Imputation for Genomic Prediction Analyses*

Genotyping of the panel and SNP quality were described in our previous study [25]. In brief, maize leaves at the sixth-leaf stage of development were harvested from 341 of the 358 lines of the panel (5–10 plants per line) in 2017B and 2018A (for lines not captured in 2017B). The leaf samples were oven-dried overnight at 36 degrees Celsius and shipped to the Biosciences east and central Africa (BecA) Laboratory of the International Livestock Research Institute (ILRI, Nairobi, Kenya) for DNA extraction and genotyping. Diversity Array Technology (DArT) genotyping facilities (44) were used to successfully identify 34,509 SNPs from 341 of the 358 lines composing the panel; hence, only these lines were considered for the GP analyses. Duplicate SNPs were first removed using the R package DartR (45), leaving 28,919 unique SNPs (DRSNP) distributed across all the 10 chromosomes of the entire maize genome.

The DRSNP dataset was imputed before GP using KDCompute 1.5.2. beta (https:// kdcompute.igss-africa.org/kdcompute/home), an online database developed by Diversity Array Technologies (DArT, https://www.diversityarrays.com) for the analysis of DArT marker data. KDCompute uses a suit of imputation methods to impute the SNP dataset and scores the imputation results by calculating simple matching coefficient (SMC). The method with the highest SMC is considered as optimal and used to impute the original genotypic dataset.

#### *2.3. FAW and MW Resistance Phenotyping*

After germination, plants were left unprotected to allow sufficient natural pressure of fall armyworm (FAW) population to build up. FAW damage scoring in all the three environments was carried out two months after planting when adequate natural FAW infestation levels had manifested, and scoring was based on a visual assessment using a scale of 1 (no or minor leaf damage) to 9 (all leaves highly damaged) [75], illustrated in Figure S1 [18].

Rearing of and bioassays for MW were performed as described in previous experiments carried out at NaCRRI [76,77]. Weevils were reared prior to the MW bioassay to obtain enough insects aged between 0 to 7 days for infestation. During rearing, standard conditions were provided for weevils to ensure proper acclimatization during the experiment. Rearing was carried out by preparing a weevil-maize grain culture of 300–400 unsexed insects and 1.5 kg of grains contained in 3000 cm<sup>3</sup> plastic jars incubated for 14 days in the laboratory at a temperature of 28 ± 2 ◦C and relative humidity of 70% ± 5%, to enhance oviposition. The lids of the jars were perforated and a gauze-wire mesh with a pore size smaller than one mm was fitted on each lid to allow proper ventilation while preventing the weevils from escaping.

After harvesting and shelling, grains of each line were bulked across environments. Then, samples of 30 g were weighed from each grain bulk, aiming to produce three replicates per line for the MW bioassay experiment. However, due to the lack of an adequate amount of grains for most of the lines of the panel, only 64, 123, and 132 lines could generate three, two, and one replicates, respectively, and were therefore considered for the MW bioassay experiment. Each of these samples was wrapped in polythene bags and kept at −20 ◦C for 14 days to eliminate any weevil infestation prior to the start of the experiment. After this disinfestation process, samples were left to thaw and transferred into 250 cm3 glass jars and infested with 32 unsexed weevils. After 10-days of incubation to allow oviposition, all dead and living adult insects were removed. One month after infestation (MAI), each sample was removed from its jar, and the grains and the flour were isolated and their weights were recorded. The total number of holes inflicted by the weevils on the grains was counted along with the number of grains affected by such damages. Additionally, the numbers of dead and living weevils were recorded. After these measurements were collected, the grains were returned to their respective jars and all the measurements were repeated at two and three MAI. The collected data were used to infer, for each sample, the cumulative grain weight loss (GWL), the cumulative number of emerged adult weevil progenies (AP), and the final number of damage-affected kernels (AK).

#### *2.4. Statistical Analyses of the Phenotypic Data*

Both best linear unbiased estimators (BLUEs) and predictors (BLUPs) were generated using the general linear model with only phenotype option of the software Trait Association through Evolution and Linkage (TASSEL) [78] and the *ranef* function of the R package [79] lme4 [80], respectively. The mixed linear model for generating BLUEs (all factors considered as fixed) and BLUPs (all factors considered as random) for MW traits (GWL, AP, AK, NH, and FP) was as follows:

$$Y = \mu + \text{Replication} + \text{Genotype} + Error$$

The mixed model for generating BLUEs (all factors considered as fixed) and BLUPs (with all factors considered as random) for FAW damage scores across environments model was:

*Y* = *μ* + *Location* + *Block* + *Genotype* + *Location* : *Genotype* + *Error*

where *μ* in the two equations is the intercept.

#### *2.5. Strategies for TS and BS Determination*

#### 2.5.1. MW Resistance Traits

Due to inadequate amount of seeds, only 37% (126 out of 341 that had genotypic data) of lines from the panel had phenotypic data on grain weight loss (GWL), adult progeny emergence (AP), and number of affected kernels (AK). Therefore, to estimate GP accuracies for MW resistance, these 126 lines were used as the TS and the remaining 215 lines with only genotypic data constituted the breeding set (BS).

#### 2.5.2. FAW Damage Resistance

The GP analyses for FAW resistance were carried out on the 341 lines of the panel that were genotyped and phenotyped for FAW damage resistance. To determine TS and BS sizes and compositions for the evaluation of maize resistance to FAW damage, two strategies, namely, random-based TS (RBTS) and pedigree-based TS (PBTS), were used.

#### 2.5.3. Random-Based TS Determination

For the RBTS, 126 (37%) lines used for GPs of MW-resistance traits were used as the TS for FAW to predict the GEBVs of the remaining 215 lines first. To build the second TS for FAW, the 215 (63%) lines used earlier as BS were considered as a TS. Then to determine the third and fourth TSs for FAW, random selections of 75 and 85% of the lines in the entire panel were performed through the Excel formula "*=INDEX(\$A:\$A,RANDBETWEEN(1,COUNTA (\$A:\$A)),1)*" and dragging until the adequate number of lines for each percentage determined above was obtained.

#### 2.5.4. Pedigree-Based TS Determination

The four datasets determined based on the pedigrees of the lines in the panel (PBTS strategy) are presented in Table 2. For the first dataset (FAW.Ped1), the 235 (68.91%) CIMMYT doubled haploid (DH) lines were used as a TS and the remainder (106 lines) as a BS. Regarding the second dataset, the TS and BS were switched to consider the TS in FAW.Ped1 as BS, and BS in FAW.Ped1 as the TS. The third dataset, FAW.Ped3, had a TS composed of the 294 that were neither stem borer (SB) resistant nor storage pest (SP)-resistant lines from CIMMYT, whilst the 28 SB and 19 SP-resistant lines from CIMMYT constituted the BS. In the last dataset, FAW.Ped4, the 235 DH lines, the 28 SB and 19 SPresistant lines from CIMMYT, and the five SB-resistant lines from IITA amounting to 287 (84.16%) genotypes were considered as the TS and the remaining 54 lines from NaCRRI lines were considered as the BS (Table 2).

**Table 2.** Compositions of the pedigree-based test sets (TSs) for fall armyworm (FAW) datasets.

DH = doubled haploid; FAW, fall armyworm; FAW.Ped1 to 4, FAW datasets 1–4 with TS based on pedigree information of the lines in the panel; SB, stem borer; SP, storage pest; TS, training set; CIMMYT, International Center for Maize and Wheat Improvement; IITA, International Institute for Tropical Agriculture.

#### *2.6. Genomic Prediction Algorithms*

The GP analyses were performed using the BLUEs and BLUPs of the phenotypes and the 28,919 DRSNPs. Sixteen algorithms available in 10 GP methods were implemented using the Genomic Prediction 0.2.6 plugin of the KDCompute 1.5.2. beta. The 10 methods were directly translated from functions of five R packages designed for GP analyses:

#### 2.6.1. Bayesian Models

Bayesian models have different prior distributions with a general model that can be as follows: *y* = 1*nμ* + *Zμ* + *ε*, where *y* is the vector of observations, *Z* is the design matrix for random effects, and *μ* is the vector of random effects [31].

The BLR (Bayesian Linear Regression) algorithms from the BLR R Package [81] are used to fit the Bayesian ridge regression. The marker effects are assumed to have a Gaussian prior distribution with mean 0 and variance σ2, where σ<sup>2</sup> is unknown and assumed to have scaled *x*<sup>2</sup> distribution. In the KDCompute genomic prediction 0.2.6 plugin, the Gibbs sampler is run with 4000 iterations and 1000 iterations for burn-in period as default parameters.

The Bayesian Generalized Linear Regression (BGLR) package fits various types of parametric and semi-parametric Bayesian regressions. The parametric Bayesian algorithms used from this package rely on different prior distributions that induce different types of shrinkages of the marker effects [82], including: Gaussian (Bayesian ridge regression, BRR [83]), scaled-t (BayesA [26]), double-exponential (Bayesian LASSO, BL [84]), and two component mixtures with a point of mass at zero and a distribution with a slab that can be either Gaussian (BayesC [44]) or scaled-t (BayesB [43]). In the KDCompute genomic prediction 0.2.6 plugin tool, defaults parameters for running the Gibbs sampler were used: 4000 iterations and 1000 iterations for burn-in period.

Reproducing kernel Hilbert space (RKHS) [47,48] is a semiparametric Bayesian method from the BGLR package implemented on the KDCompute genomic prediction 0.2.6 plugin. The RKHS methods employs a kernel function to convert the molecular markers as a between pairs of observations distances, thereby, generating a square matrix that fits in a linear model. This non-linear regression method is expected to capture dominance and epistasis effects more efficiently. This approach can be modelled as:

$$y = \mathbb{W}\mu + \mathbb{K}\_{\text{fr}}\alpha + \varepsilon\_{\text{r}}$$

where *μ* represents the fixed effects vector and ε is a vector of random residuals. The parameters α and ε are assumed to have independent prior distributions *α* ∼ *N*(0, *Khσ*<sup>2</sup> *α* ) and *ε* ∼ *N*(0, *I σ*2 *ε* ), respectively, and the matrix *Kh* relies on a reproducing kernel function with a smoothing parameter *h*. The parameter *h* measures the genomic distances among genotypes that can be interpreted as a correlation matrix and it controls the rate of decay of the correlation among genotypes [51]. To perform this analysis, the same number of iterations and burn-in parameters as for the other Bayesian methods described above were set on the KDCompute genomic prediction 0.2.6 plugin.

#### 2.6.2. Mixed Models

The Sommer (solving mixed model equations in R) package [85] was used to implement the *mmer* (mixed model equations in R) function on the KDCompute genomic prediction 0.2.6 plugin. The package solves mixed model equations proposed by Henderson [86]. It works incidence matrices and known variance covariance matrices for each random effect using four algorithms: efficient mixed model association (EMMA) [87], average information (AI) [88], expectation maximization (EM) [89], and the default Newton–Raphson (NR) [90].

The model by Sommer can be formulated as [85]: *y* = *Xβ* + *Zμ* + *ε* with variance *V*(*y*) = *V*(*Zμ* + *ε*) = *ZGZ* + *R* Additionally, the mixed model equations for this model are:

$$
\left[\begin{array}{c}X^{\prime}R^{-1}X\,X^{\prime}R^{-1}Z\\Z^{\prime}R^{-1}X\,Z^{\prime}R^{-1}Z+G^{-1}\end{array}\right]^{-1}\left[\begin{array}{c}X^{\prime}R^{-1}y\\Z^{\prime}R^{-1}y\end{array}\right] = \left[\begin{array}{c}\mathcal{S}\\\mu\end{array}\right].
$$

where *G* = *Kσ*<sup>2</sup> *<sup>ω</sup>* is the variance covariance matrix of the random effect μ, from a multivariate normal distribution *<sup>μ</sup>* ∼ *MVN*(0, *<sup>K</sup>σ*<sup>2</sup> *<sup>μ</sup>*), *K* is the additive or genomic relationship matrix (**A** or **Ag**) in the genomics context, *X* and *Z* are incidence matrices for fixed and random effects, respectively, and *R* is the matrix for residuals (here *Iσ*<sup>2</sup> *<sup>e</sup>* ). A mixed model with a single variance component other than the error (*σ*<sup>2</sup> *<sup>e</sup>* ) can be used to estimate the genetic variance (*σ*<sup>2</sup> *<sup>μ</sup>*) along with genotype BLUPs to exploit the genetic relationships between individuals coded in **K**(**A**). The genomic relationship matrix was constructed according to VanRaden where *K* = *ZZ* /2 ∑ *pi*(1 − *pi*) [91].

The ridge regression best linear unbiased predictor (rrBLUP) packages can either estimate marker effects by ridge regression, or alternatively, BLUPs can be calculated based on an additive relationship matrix or a Gaussian kernel. Additionally, using the rrBLUP package, **the mixed model solution (MMS)** that calculates the maximum-likelihood (ML) or restricted-ML (REML) solutions for mixed models to perform GP [92] was fitted in the KDCompute genomic prediction 0.2.6 plugin.

The mixed models fitted by rrBLUP can be formulated as:

$$y = X\beta + Z\mu + \varepsilon\_r$$

where *β* is a vector of fixed effects and *μ* is a vector of random effects with variance *Var*[*μ*] = *Kσ*<sup>2</sup> *μ* . The residual variance is *Var*[*ε*] = *I σ*2 *ε* . This class of mixed models, in which there is a single variance component other than the residual error, has a close relationship with ridge regression (ridge parameter *λ* = *σ*<sup>2</sup> *<sup>ε</sup>* /*σ*<sup>2</sup> *<sup>μ</sup>*) (https://kdcompute.igss-africa.org/ kdcompute/home).

#### 2.6.3. Machine Learning Algorithms

The R package RandomForest that implements Breiman's random forest algorithm for classification and regression [93] was used on the KDCompute genomic prediction 0.2.6 plugin to fit the function missForest. Random forest is a non-linear machine learning algorithm that uses a two-layer randomization process to build decorrelated bootstrapped trees. As a first randomization layer, it builds multiple trees using a bootstrap sample of the marker data in the training. Then, a second randomization process is carried at the novel nodes to grow final trees. The random forest method selects at each node of each tree, a random subset of variables, and only those variables are used as candidates to find the best split for the node [94]. To predict the breeding value of a line in the TS, predictions over trees for which the given observation was not used to build the tree are averaged [51]. On the KDCompute 1.5.2. beta platform, both options for the *mtry*, square root and regression (*sqrt*(*p*) and *p/3*, respectively, where *p* is number of variables in x), for the classification of the number of variables randomly sampled as candidates at each split were implemented in this study. Additionally, the trees to grow (*ntree*) was set to 10, while *node size* (minimum size of terminal nodes) and *max nodes* (maximum number of terminal nodes trees in the forest can have) were set to 5 and 10, respectively. The 16 methods used in this study and their statistical characteristics are presented in Table 3.

**Table 3.** Genomic prediction methods used for the analysis of the different traits and datasets.

#### *2.7. Cross-Validations and PA Estimation*

To calculate the predictive accuracies of each of the 17 methods, a cross-validation approach was performed using the data for the TS with 10 folds and five repetitions amounting to 50 replications. The PAs were estimated as the correlation coefficient (R2) averaged across the 50 cross-validation replications between the observed phenotypic values and the predicted genomic-estimated breeding values (GEBV) (https://kdcompute. igss-africa.org/kdcompute/plugins).

#### **3. Results**

#### *3.1. Higher PAs Achieved for FAW and MW-Resistance Traits with BLUPs when Compared to BLUEs across GP Algorithms*

Both genotypic BLUEs and BLUPs for resistance to FAW and MW traits such as AK, AP, and GWL were used in GPs. In general, BLUPs produced better predictions than BLUEs by at least two orders of magnitude in terms of PAs (Figure 1). The PAs realized with BLUEs (Figure S2) varied from −0.246 for FAW (mms\_ML) to 0.299 for AP (BayesB), while PAs for BLUPs ranged from 0.668 for AP (mmer\_NR) to 0.823 for AP (missForest\_Reg). The differences in terms of accuracies between BLUEs and BLUPs were high, despite the highly significant (*p* < 0.001) correlations between BLUEs and BLUPs for each trait ranging from 0.93 for FAW to close to 1 for AP, AK, and GWL (presented in Figure 1); therefore, only results for BLUPs will be presented hereafter.

**Figure 1.** Boxplot of PAs (prediction accuracies) for best linear unbiased estimators (BLUEs) (in pink) and predictors (BLUPs) (in blue) of maize resistance to FAW and MW across prediction models and correlations (r) between BLUEs and BLUPs for each trait. FAW, fall armyworm; GWL, grain weight loss; AP, adult progeny emergence; AK, number of affected kernels. \*\*\* significant at *p* < 0.001.

#### *3.2. PAs for MW Resistance Traits Using BLUPs*

The PAs were generally high for the tested MW traits, mostly above 0.668 across the 12 GP models that were successfully run on the datasets (Figure 2); however, RKHS failed to work for AK. The highest PAs were achieved for AP with missForest\_reg (0.823), followed by BRR (0.805), and RKHS (0.804), whilst mmer\_NR algorithm had the lowest prediction accuracy of 0.667 (Figure 2). The PAs achieved for GWL ranged from 0.742 for missForest\_Sqt to 0.795 for mmer\_NR, while for AK, they varied from 0.749 for missForest\_sqrt to 0.779 for BRR (Figure 2). In general, Bayesian models predicted better than both mixed model and machine learning methods, although the differences were small (Figure S3).

**Figure 2.** Boxplots of the genomic prediction accuracies of BLUPs for MW-resistance traits: GWL, grain weight loss; AP, adult progeny emergence; AK, number of affected kernels (See Table 3 for GP algorithms).

#### *3.3. PA for FAW Resistance Using BLUPs*

The different maize resistance to FAW datasets showed high predictive abilities with 10 of the 16 GP algorithms used in the study. For the RBTS approach, the PAs were lowest with the dataset that had a TS composed of 37% (lowest size) of the panel and highest with the largest TS (85% of the panel). Even with a TS of 37%, the PAs were still high, ranging from 0.694 to 0.714 for mms\_ML and BLR methods, respectively (Figure 3). However, it should be noted that with equal TS sizes and same composition (37% of the panel), higher PAs were achieved for MW-resistance traits (GWL, AP, and AK) compared to FAWresistance ones (Figure S3). The PA for the RBTS of 63% varied from 0.833 for BL method to 0.838 for the missForest\_Sqt; thus, there was a small variation among different methods. Similarly, there was minimal variation among GP algorithms on the dataset with a 75% TS whose PAs varied from 0.838 for mms\_REML to 0.843 for MissForest\_Reg. The same trend was obtained on the dataset with a RBTS of 85% of the panel, with PAs ranging from 0.843 for the BRR model to 0.847 for the missForest\_Reg method. Furthermore, there was a high and significant (*p* < 2.2.10−16) positive correlation of 0.92 (Figure 4) between the PAs and TS sizes for FAW datasets for the RBTS denoting a steady improvement of the PAs as the TS size increased. However, the PAs for FAW resistance reached a plateau at TS size above 63% of the panel (Figure 5).

**Figure 3.** Boxplot of PAs for maize resistance to the fall armyworm (FAW) datasets with the RBTS approach with random selection of 37, 63, 75, and 87% of the entire panel (see Table 3 for GP algorithms).

**Figure 4.** Pearson correlation between training set (TS) sizes and prediction accuracies (PAs) across the 10 genomic prediction algorithms conducted on RBTS (**A**) and PBTS (**B**) datasets for fall armyworm resistance (FAW) resistance.

**Figure 5.** Prediction accuracies for FAW with RBTS across algorithms and training sets with different sizes in percent of the total panel.

Although the PAs did not vary much among GP algorithms, especially when the analyses involved larger TS sizes equal or bigger than 63% of the panel, the machine learning methods slightly outperformed other GP algorithms for all the traits, except for the TS of 37% where Bayesian methods such as BLR and BayesC showed a slight advantage over the machine learning methods (Figure S4). The PAs for FAW-resistance datasets with PBTS were generally high, mostly above 0.82 (Figure 6). For the first dataset (FAW.Ped1) with a TS of 68.91% of the panel (see Table 2), the PAs varied between 0.828 for BLR to 0.835 for missForest\_Sqt. For FAW.Ped2 (TS = 31.09%), the PAs ranged from 0.862 for BayesC to 0.864 for mms\_REML.

**Figure 6.** Boxplots of PAs for maize resistance to the fall armyworm (FAW) datasets using the PBTS approach (see Table 2 for the PBTS strategy and Table 3 for GP algorithms).

For FAW.Ped4, with a TS of 84.16%, PAs varied between 0.860 to 0.864 for missForest\_Sqt and mms\_ML, respectively. However, for FAW.Ped3 with the largest TS (86.22%), eight of the 10 algorithms achieved low PAs (below 0.20) and only missForest\_Reg and missForest\_Sqt attained PAs of 0.749 and 0.750, respectively. Thus, the Pearson correlation between the sizes of the PBTS datasets and the predictions accuracies for the 10 GP algorithms revealed a significant (*p* > 0.0036) negative relationship of r = −0.45 (Figure 4).

In the FAW datasets, the PAs were more influenced by the composition of the TS and its genetic relationship with the BS (see Table 2). Using the doubled haploid (DH) lines as TS (FAW.Ped1) and vice-versa (FAW.Ped2) or DH and stem borer (SB) and storage pest (SP)-resistant lines as TS (FAW.Ped4) permitted achieving relatively high PAs from all the 10 algorithms, which when considering the CIMMYT SB and SP-resistant lines as BS and the remainder as a TS (FAW.Ped3), only resulted in machine learning algorithms missForest\_reg and missForest\_Sqt achieving relatively high PAs. Furthermore, the composition of the TS and its relationship with the BS determined which GP methods achieved the highest Pas; machine learning algorithms worked best on FAW.Ped1 and FAW.Ped3, linear mixed model approaches outperformed Bayesian and machine learning algorithms on FAW.Ped2 and FAW.Ped4, and Bayesian methods ranked either second or third on all datasets (Figure S5). It should be noted that the PBTS strategy generally achieved better PAs than the RBTS irrespective of the size of the TS, except for the FAW.Ped3 dataset (Figures 3 and 6).

#### **4. Discussion**

Tropical maize germplasm is characterized by rapid linkage disequilibrium (LD) decay with high diversity [95]. These germplasm genetic characteristics make genomic selection (GS) a promising approach to integrate into African breeding programs [96]. However, genomic prediction (GP) models are very diverse and their differential performances depend on crops and trait architectures, besides other parameters such as the size of the training set (TS) and its genetic relationship with the breeding set (BS) [31,37]. Therefore, this study aimed at assessing the feasibility of genomic selection for maize resistance to FAW and MW through estimation of the genomic prediction accuracies achieved by parametric, semiparametric, and nonparametric (machine learning) genomic prediction (GP) algorithms using phenotypic BLUEs and BLUPs, and random and pedigree-based TS determination strategies.

#### *4.1. Higher Pas Were Achieved for BLUPs Compared to BLUEs for Both FAW and MW-resistance Traits*

With a RBTS of 37% of the panel, which was the smallest and expected to give the worst PAs, PAs were higher (at least two-fold) across both FAW and MW-resistance traits and for all GP models when trait BLUPs were used as phenotypes compared to BLUEs, although there were high Pearson correlations between these two categories of phenotypic data for each trait. In general, BLUPs were reported to have higher predictability than BLUEs owing to better accounting for outliers and environmental variabilities permitted by the shrinkage procedure in BLUPs, which results in more accurate estimates of individual genetic effects [70–73]. Furthermore, most of the predictive differences between BLUPs and BLUEs might have stemmed from BLUPs being more suitable than BLUEs in fitting data recorded from unbalanced experiments [49,70] as was the case for both FAW damage scores across environments and MW bioassay in this study. Therefore, for all subsequent analyses with higher RBTS sizes and the PBTS strategy for FAW, only BLUPs were focused at in this study and will be further discussed.

#### *4.2. High PAs Were Achieved for FAW and MW-Resistance Traits Using Moderately Sized Training Sets*

The obtained PAs were high for both MW and FAW-resistance traits even with TS of moderate sizes confirming the potential of genomic selection (GS) in Africa-adapted germplasms [28–30,33]. With a TS of 37% of the entire panel, high PAs (above 0.70) for MWresistance traits, grain weight loss (GWL), adult progeny emergence (AP), the number of

affected kernels (AK), and FAW resistance were achieved in agreement with the moderate to high heritability values for these traits as, reported earlier [21,31,41]. These results are significantly important considering that one of the disadvantages of GS is the requirement of large TS which negatively impacts the reduction of phenotyping cost [62,64].

The PAs increased up to above 0.85 in proportion to the increase in TS (RBTS approach) size for FAW resistance which was the only trait phenotyped for all the lines of the panel. It would be interesting to phenotype other lines of the panel that were not evaluated for MW-resistance traits to establish larger TS which may improve the PAs [31,65,97,98]. Very few reports of GP are available for maize resistance to biotic stresses. High PAs were achieved for maize resistance to chlorotic mottle virus (up to 0.95) and maize lethal necrosis (reaching 0.87) in tropical germplasm [67]. However, lower PAs of up to 0.59 were obtained in a study that assessed the predictability of maize resistance to the European corn borer [99] in temperate germplasm. Additionally, Gowda et al. (69) reported moderate PAs (close to 0.60) for maize resistance to a biotic stress, maize lethal necrosis in tropical maize populations.

#### *4.3. GP Algorithms Performed Differently on FAW and MW Maize Resistance Traits*

In this study, several GP models that included statistical and machine learning algorithms from parametric, semi-parametric, and nonparametric approaches were used to predict FAW and MW-resistance traits. These GP algorithms, as expected, performed differently on the different traits although the predictive variations were generally minimal, especially when large TS were involved, similarly to earlier model benchmarking reports [100,101]. Bayesian models (parametric: BLR and BRR, and semi-parametric: RKHS) performed better on MW traits, GWL, AP, and AK, while nonparametric machine learning algorithms (missForest, here), and to a lesser extent, the linear mixed model (especially in the PBTS approach), achieved the highest PAs on FAW datasets. The differential performances of the different GP algorithms on the insect resistance traits evaluated in this study could be due to differences in the genetic structures (extent of additive vs. non-additive gene action) of the respective traits [23,38,47,49]. Maize resistance to FAW, which was moderately heritable across environments [25], would be expected to be controlled by both additive and non-additive genetic factors, including epistasis [102–104], whereas, MW-resistance traits such as GWL, AP, and AK with heritability values above 90% [25] were most likely characterized by a prevalence of additive gene action [105,106] in the current panel.

This supposed genetic architecture difference between FAW and MW-resistance trait could be the reason for non-linear methods such as random forest performing better at predicting FAW resistance, since these are more capable of integrating epistasis in the statistical modelling [27,51]. However, the RKHS algorithm, also a non-linear GP approach known to efficiently handle epistatic genetic relation [51,59], did not successfully run on FAW dataset, although it was among the best models for predicting MW-resistance traits, except BLUPs for the number of affected kernels (AK), for which the RKHS algorithm did not run successfully. In this study, the reasons for some GP algorithms failing to run either on MW or FAW-resistance datasets are unclear, but this could be related to the BLUPS structure of the datasets that failed to run. It should be noted that all the algorithms ran successfully on phenotypic BLUEs datasets with the smallest TS (37% of the panel) being used to compare PAs between BLUPs and BLUEs in this study. However, the two to three-fold predictive ability gain with BLUPs compared to BLUEs would be an incentive to consider BLUPs in future GS activities for maize resistance to MW and FAW. Overall, future GS efforts for maize resistance to MW and FAW are recommended to focus more on Bayesian and machine learning algorithms such as random forest, BayesA, BayesB, BayesC, BRR, and BLR which outperformed mixed linear models for most datasets considered in the current study.

#### *4.4. Influences of the Sizes and the Compositions of TS and BS on PAs*

Two factors, the relative sizes of the TS and BS (RBTS approach) and their genetic relationship (PBTS approach), influenced the levels of PAs across FAW-resistance datasets, corroborating earlier reports [31,63,65,97,98,107,108]. A net increase in PAs for maize resistance to FAW was realized when the size of the TS was increased from 37% (0.694 to 0.714) to 63% (0.833 to 0.838), similar to earlier reports on wheat yield [109]. This increase was followed by a slight gain in predictability at 75% (0.837 to 0.843) and 85% (0.843 to 0.847), and thus, the PAs plateaued when TS sizes above 63% were considered in this study as reported earlier in other studies [21,64,109–111]. Thus, future GS programs for maize resistance to FAW could be designed around TS composed of a minimum of 60% of the entire breeding germplasm to achieve high genetic gains. These results were further supported by the highly significant (*p* > 2.2.10−16) positive correlation (R = 0.92) between TS size and PAs. Similarly, positive correlations between the number of lines in the TS and the PAs, and plateau for the PAs were also reported by Edwards et al. [109].

The composition of the TS and its relationship with the BS are determinant factors for the genomic predictability of complex traits [63,112–114]. In the current study, using the PBTS approach, these two parameters were more important than the size of the TS since higher PAs were achieved in FAW.Ped2 (0.862 to 0.864) with a TS of 31.09% compared to all other FAW PBTS datasets, including FAW.Ped3 (0.114 to 0.750), with the largest TS of 66.22%. In fact, FAW.Ped3 achieved the lowest PAs among all the PBTS FAW datasets. These results were further illustrated by the significantly (*p* < 0.0036) negative correlation (R = −0.45) between the sizes of the PBTS and the achieved PAs.

However, it is not very clear why the predictions for the BS FAW.Ped3 (47 CIMMYT SB and SP-resistant lines) and the TS (DH, IITA SB, and NaCRRI lines) led to lower PAs for FAW.Ped3. A possible explanation could be that these two sets were distantly related since only two and one CIMMYT SB and SP-resistant lines, respectively, were used as parents to develop the DH lines. Spindel et al. [111] argued that high PAs can be achieved with smallsized TS when lines in the TS and the BS are closely related, since such TS would sample the full genetic diversity of the population. However, the more distantly related the TS and the BS are, the larger the required TS size to reach high PAs [111]. Using the CIMMYT SB and SP-resistant lines as a TS would most likely lead to lower PAs since such a TS would be additionally disadvantaged by its small size (47 lines). The DH lines in the current study are involved as a TS in most of the best performing GP datasets evaluated in the current study (both in the RBTS and PBTS approaches) and as unique lines in the BS of the best performing pedigree-based BS (FAW.Ped2). This DH population could be of interest in future breeding activities targeted at improving insect resistance in maize [23,115–117] and potentially useful for GS of complex traits with low to moderate heritability [118].

#### **5. Conclusions**

This study assessed prediction accuracies of genomic-estimated breeding values for fall armyworm (FAW) and maize weevil (MW)-resistance traits in a diverse Africaadapted maize panel using several parametric, semi-parametric, and non-parametric genomic prediction models. Prediction accuracies for maize resistance to FAW and MW traits were relatively high, even with a moderate training set size. For FAW resistance, although the prediction accuracies were positively correlated with the size of the training set, the composition and the relationship of the training set with the breeding set were more influential in predicting line performance. Additionally, TS determination-related parameters were more important than the type of genomic prediction models in predicting FAW and MW-resistance traits. However, Bayesian models on MW-resistance traits and machine learning models on FAW damage resistance outperformed mixed linear models in almost all the datasets used in this study. Therefore, future genomic selection programs for maize resistance to insect pests such as FAW and MW in Africa should put more effort into designing effective training sets and use selected Bayesian and machine learning GP algorithms to improve genetic gains, shorten breeding cycles, and accelerate variety release. Such programs could greatly benefit from using the genetically diverse maize panel used in this study as a base population, since it consists of lines adapted to several African agro-ecologies.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2223-774 7/10/1/29/s1, Figure S1: Rating of maize plants based on foliar damage by FAW, Figure S2: Boxplot of PA for best linear unbiased estimators (BLUEs) of maize resistance to the fall armyworm (FAW) and maize weevil (MW) with identical training set size (37%) and compositions, Figure S3: Comparisons of genomic prediction accuracies of the three best algorithms for best linear unbiased predictors (BLUPs) of maize weevil resistance traits: number of affected kernels (AK), adult progeny emergence (AP), and grain weight loss (GWL) vs., fall armyworm resistance dataset with identical TS, Figure S4: Genomic prediction accuracies of the three best algorithms for each fall armyworm resistance BLUPs datasets with RBTS of 37, 62, 75, and 85% of the entire dataset, Figure S5: Genomic prediction accuracies of the three best algorithms for each fall armyworm resistance BLUPs datasets with PBTS, Table S1: Descriptions of parents and crosses that constituted the doubled-haploid population.

**Author Contributions:** Conceptualization, A.B. (Arfang Badji), P.R., S.K., M.O., D.B.K., and L.M.; methodology, A.B. (Arfang Badji), D.B.K., and L.M.; investigation, A.B. (Arfang Badji) and D.B.K.; formal analysis, A.B. (Arfang Badji) and L.M.; resources, A.B. (Arfang Badji), G.A., M.O., D.B.K., and L.M.; visualization, A.B. (Arfang Badji); supervision, P.R., S.K., M.O., and L.M.; project administration, P.R. and M.O.; funding acquisition, A.B. (Arfang Badji), M.O., G.A., D.B.K., and P.R.; writing—original draft preparation, A.B. (Arfang Badji); writing—review and editing, all authors (A.B. (Arfang Badji), L.M., D.B.K., F.K., D.O., N.M., S.A., A.I., A.B. (Astere Bararyenya), S.N.N., T.O., P.W., M.O., M.O.-S., H.T., G.A., S.K., P.R.). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the capacity building competitive grant training the next generation of scientists provided by Carnegie Cooperation of New York through the Regional Universities Forum for Capacity Building in Agriculture (RUFORUM: RU/2016/Intra-ACP/RG/001). A. Badji received a Ph.D. scholarship from the Intra- ACP Academic mobility for Crop Scientists for Africa Agriculture (CSAA) project. Genotyping of the lines was carried out through a project of D.B.K. thanks to the Integrated Genotyping Service and Support (IGSS) coordinated by the International Livestock Research Institute (ILRI) and Bioscience east and central Africa (BecA), grant number: IGSS-DL0274. The National Crops Resources Research Institute (NaCRRI) of Namulonge, UGANDA through a grant of the USAID Feed-the-Future Uganda, Agriculture Research Activity/Maize paid the article processing charges. Further, NaCRRI financially and logistically supported field and laboratory activities of this research.