Significant Parent-of-Origin Effects for Seed, Cotyledon, and Early Plant Growth Traits in Cucumber

Parent-of-origin effects have long been recognized and exploited in plant breeding and genetic studies. These effects can be conferred by preferential expression of an allele from one parent, organellar effects, or altered organellar-nuclear interaction. The goal of this work was to evaluate parent-of-origin effects on seed, cotyledon, and early growth traits in cucumber using a full eight-by-eight diallel from crossing two doubled haploids (DHs) extracted from each of four cucumber populations. Significant general combining ability (GCA), specific combining ability (SCA), and reciprocal effects were observed for all traits, and direction and magnitude of effects were DH rather than population specific. Transcriptome analyses of reciprocal hybrids with and without significant reciprocal effects for early plant growth revealed that different pathways were associated with the significant reciprocal differences. These findings are consistent with the DH-specific nature of combining abilities and reciprocal effects across cucumber populations. Because reciprocal effects were DH and hybrid-combination specific, cucumber breeders should generate and evaluate both hybrids from reciprocal crossing for improved hybrid development.


Introduction
Parent-of-origin effects have long been recognized in plant breeding and genetic studies [1][2][3][4]. Parent-of-origin effects occur when a particular phenotype or allele is asymmetrically expressed in the progeny based on its inheritance from one parent over the other. These effects can be conferred by preferential expression or silencing of an allele from one parent, specific organellar effects, or altered organellar-nuclear interaction [4][5][6]. Allelic parent-of-origin effects are primarily caused by genomic imprinting where epigenetic regulation of gene expression occurs in a parent-dependent manner [7]. Examples of genomic imprinting are more widespread in mammals, and imprinting in plants has primarily been reported for endosperm and developmental traits [8][9][10][11][12]. Examples of genomic imprinting in plant embryos have been reported in maize, rice, and Arabidopsis, but these effects dissipate during early embryo or seedling development [13][14][15][16]. For example, only the maternal allele of the maternally expressed in embryo 1 (mee1) gene was expressed in maize embryos and endosperm, associated with transient differential methylation of the alleles, and neither allele was expressed in the zygote [13]. Raissig et al. [15] found monoallelic gene expression at nine loci in the embryo of Arabidopsis, of which eight were maternal and one paternal, all of which had either biallelic or no expression by the early seedling stage.
Plants contain two organellar genomes, plastid and mitochondrial, in addition to their nuclear genome. Organelle-nuclear interaction and coordination are crucial for normal plant development, growth, and productivity, as the majority of essential genes for organelle function are encoded by the nuclear genome [17][18][19][20]. Organellar effects are one type of parent-of-origin effect and are commonly referred to as cytoplasmic effects. One of the most economically valuable cytoplasmic traits is cytoplasmic male sterility (CMS), which has been identified and exploited for hybrid development across species [21][22][23]. Various patterns of leaf variegation are commonly identified as cytoplasmic variants, the majority resulting from sorting of organellar variants [24,25]. Organelle-associated herbicide resistance has also been reported across species [26,27]. In addition to discrete traits, Dobler et al. [4] found that cytoplasmic effects are widely underestimated and underreported, with evidence from a large meta-analysis reporting prevalent cytoplasmic effects across taxa and traits. Joseph et al. [5] and Flood et al. [28] reported significant cytoplasmic effects and cytoplasmic-nuclear interactions across various traits related to plant productivity and growth in Arabidopsis. Joseph et al. [5] found that cytoplasmic effects on natural variation in the metabolome of Arabidopsis not only affected a higher number of metabolites than identified nuclear QTL, but also had higher average effect sizes.
Identification of beneficial parent-of-origin effects is valuable for plant breeders for production of better performing hybrids. Evaluation of a full diallel in cucumber allows for elucidation of parent-of-origin effects based on specific combinations of the organellar genomes and/or unique organellar-nuclear interactions. Shen et al. [43] generated a diallel mating scheme with reciprocal crosses among cucumber doubled haploid (DH) lines, producing reciprocal hybrids with identical nuclear genotypes but differing for their chloroplast and mitochondrial genomes. Shen et al. [43] reported significant differences in performance between reciprocal hybrids of DH lines. For example, significant reciprocal effects were reported for DH line TMG1, resulting in more vigorous progeny when this DH was used as the female parent [43]. These results suggest positive effects of the chloroplast genome from DH TMG1, which supports use of this DH line as a female parent in hybrid development. However, significant maternal effects on plant growth could also be affected by seed size, which has been reported in bean [44], pea [45,46], maize [47], and squash [48]. Shen et al. [43] reported significant paternal effects for DH line GY14, potentially due to a positive effect of the mitochondrial genome and supporting the use of this DH as a male parent. Because Shen et al. [43] used a single DH line from each of four cucumber populations for diallel crossing, it is not certain if significant reciprocal cross differences are an attribute of the DH or population from which the DH was extracted. Evaluation of multiple DH lines in a full diallel scheme would reveal population versus DHspecific effects and provide further insight into these interactions and their potential use in breeding programs. Additionally, gene expression analysis of significant reciprocal-cross effects could reveal potential candidate genes and pathways conferring beneficial nuclearorganelle interaction for enhanced performance. If hybrids with significant reciprocal cross effects share gene expression patterns, this could validate resource allocation into further analyses of the genetic basis of these parent-of-origin effects and reveal selection targets for beneficial nuclear-cytoplasmic interactions to potentially increase performance of hybrids. In addition to evaluating reciprocal effects, which indicate variation across reciprocal hybrids, parents and hybrids can be evaluated based on general combining ability (GCA) and specific combining ability (SCA), respectively. GCA refers to the relative contribution of a single parent to a trait across multiple crosses, while SCA is the relative effect size of a specific combination of parents to a trait [49]. The goals of this work were to use a diallel mating scheme using DH lines of cucumber to (1) identify GCA, SCA, and reciprocal effects across seed, cotyledon, and early plant growth traits to identify beneficial parent-of-origin effects on these traits and (2) identify if there are common patterns of nuclear gene expression in hybrids with significant reciprocal effects.

Materials and Methods
Two DH lines were extracted from each of four cucumber populations: North American pickling cucumber 'GY14' (GY14-9 and GY14-15), Asian slicing cucumbers '9930' (9930-3 and 9930-5) and 'TMG1' (TMG1-4 and TMG1-5), and North American slicing cucumber 'Straight 8' (ST8-2 and ST8-4) by culturing of immature female flowers [43]. Homozygosities of all DH lines were confirmed using 50 simple sequence repeats (SSRs) distributed across the nuclear genome [39]. All DH lines were grown in the greenhouse and crossed in a full diallel mating scheme without self-pollination to produce 56 hybrids. Seeds were harvested from two fruits from each cross. Individual seeds were extracted from each fruit and placed into 48-well microtiter plates (18 seeds per fruit × 2 fruit per cross = 36 seeds per cross). Individual seed length, width, perimeter, and area were determined by image analysis using ImageJ as previously described [50].
The identity of each imaged seed was tracked through planting and subsequent plant measurements in order to estimate correlations among seed size, cotyledon size, and early plant growth. Seeds from 112 fruits (56 hybrids × 2 fruits per cross) were sown in late June through early July into a soilless substrate (PRO-MIX HP Mycorrhizae; Premier Tech Horticulture, Quakertown, PA, USA) in 48-cell trays and covered with vermiculite. Trays were placed onto heating pads at 28 • C for 4 days for uniform germination. Six days after sowing, cotyledon length and width were measured using a digital caliper. All plants were transplanted into 11.4 cm diameter pots one week after sowing. Plants were then grown in greenhouses at 28 • C day and 24 • C night temperatures with a supplemental lighting by sodium vapor lamps for 16-h days. At 20 days after sowing, plants were destructively harvested by cutting just above the cotyledons. Fresh weight was immediately measured for each sample using a digital balance, and samples were then dried for at least 4 days at 60 • C. Dry weights were then measured with a digital balance.
The experiment was a randomized complete block design (RCBD) with three blocks (benches in a single greenhouse) and the experiment was repeated three times (three different greenhouses) in greenhouses on the campus of the University of Wisconsin-Madison. Six seeds from each fruit from each of the 56 hybrids were sown and two plants from each fruit were randomly placed into each of three blocks for each replicate.
Pearson's correlation coefficients were calculated for all pairwise comparisons of traits using the corrplot function of the corrplot package in R [51]. Significance was tested and p-values were adjusted using the Bonferroni correction using the cor.test function in the stats package in R (https://www.R-project.org (accessed on 15 April 2020)). Means were calculated for all traits from the four plants (2 plants per fruit × 2 fruits) per hybrid within each block. Analysis of variance (ANOVA) was performed based on Method 3, Model 1 from Griffing [52], including effects of experiment, experiment by genotype interaction, and block nested within experiment. The statistical model used based on Griffing [52] was Y ijkl = µ + H ij + E k + (H × E) ijk + B(E) kl + ε ijkl with H ij = G i + G j + S ij + R ij , where µ is the grand mean for each trait (Y); H ij is the effect of the hybrid genotype from the i × j cross; E k is the effect of the kth experiment; (H × E) ijk is the interaction of the ijth hybrid genotype and kth experiment; B(E) kl is the effect of the lth block within the kth experiment; ε ijkl is the experimental error for the ijklth observation; G i and G j are the GCA effects for the i and j parents, respectively; S ij is the SCA effect for the i × j hybrid; R ij is the reciprocal effect for the i × j hybrid; such that the sum of all G = 0, S ij = S ji , and R ij = −R ji . The GCA, SCA, and reciprocal effects of each hybrid were calculated for each trait and all analyses were performed using an SAS (version 9.3; SAS Institute, Cary, NC, USA) program for diallel analysis [43]. Statistical significance of GCA, SCA, and reciprocal effects were determined Agronomy 2021, 11,1908 4 of 14 using the t-test and variance calculations as described by Griffing [52]. GCA and SCA were tested based on deviations from zero. Significance of reciprocal effects was tested based on deviation between positive and negative effects for each reciprocal cross.
Two sets of reciprocal hybrids with significant reciprocal effects (hybrids of ST8-4 with 9930-3 and hybrids of ST8-2 with GY14-15) and two sets with non-significant reciprocal effects (hybrids of ST8-4 with GY14-15 and hybrids of ST8-4 with TMG1-5) were selected for RNA-seq analysis. Nine seeds of each hybrid were planted and grown as described above. At 20 days after planting, all plants were destructively harvested. The apical meristem and youngest expanding leaf from each of three plants from each hybrid were combined and immediately frozen in liquid nitrogen and then ground in liquid nitrogen using a mortar and pestle. Tissues from three plants from each hybrid were combined to create a biological replication and three biological replicates were generated for each hybrid yielding 24 samples (8 reciprocal hybrids × 3 biological reps) for sequencing.
Total RNA was extracted from each sample according to the manufacturer's protocol using the Zymo RNAeasy kit (Zymo Research, Irvine, CA, USA). Total RNA sample quality and concentration were determined by electrophoresis through 1% agarose gels and using the NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), respectively. Samples were stored at −80 • C until shipment on dry ice to Novogene (Sacramento, CA, USA) for library construction and mRNA sequencing. At Novogene, total RNA samples underwent a second set of quality control procedures to test purity using NanoPhotometer spectrophotometer (IMPLEN, München, Germany), RNA degradation and contamination using 1% agarose gel electrophoresis, and RNA integrity and using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). After quality control, mRNA was enriched using oligo(dT) beads and a cDNA library was constructed with 250-300 bp inserts. Libraries were then sequenced on Illumina NovaSeq platforms with paired-end 150 bp (PE 150) sequencing.
Raw fastq files were trimmed to remove adapter sequences using Skewer [53]. After trimming, reads were filtered for genes with zero or low-abundance and normalized by the method of trimmed mean of M-values (TMM) [54]. Paired-end reads were aligned to the 9930 cucumber nuclear genome reference (ASM407v2) using STAR (Spliced Transcripts Alignment to a Reference) [55]. Mapped reads were quantified using RSEM (RNASeq by Expectation Maximization), implementing TPM (Transcripts per million mapped reads) for accurate comparison across libraries [56]. Analysis of differentially expressed genes (DEGs) between pairs of reciprocal hybrids were performed using the glm function in the edgeR package [57]. Significance was assessed based on the adjusted p-value using a Benjamini-Hochberg correction to control for the false discovery rate (FDR) [58].
DEG lists were compared using Venny 2.1 online analysis tool [59]. DEGs were classified based on the reference Ensembl annotation and the PANTHER (Protein Analysis Through Evolutionary Relationships, http://pantherdb.org (accessed on 10 November 2020) classification system using the PANTHER tools Gene List Analysis [60][61][62]. Gene Ontology (GO) enrichment analysis was performed for sets of differentially expressed genes using the PANTHER tools statistical overrepresentation test and significance was tested based on Fisher's exact test with a correction for FDR [60][61][62]. Overrepresentation in gene sets was compared to a custom reference gene list which included all genes detected with at least two TTM-adjusted read counts across the experiment.

Trait Correlations, Combining Abilities, and Reciprocal Effects
Traits measured across 56 hybrids were seed area, length, width, and perimeter, cotyledon length and width, and plant fresh and dry weights 20 days after seed sowing. The identity of each seed was tracked to determine correlations among seed, cotyledon, and plant growth traits. Seed area had significant (p < 0.001) correlations with all other seed traits with Pearson coefficients ranging from 0.81 to 0.92 (Table 1), as expected, because larger seeds would have greater lengths, widths, and perimeters. Correlations of seed traits with cotyledon sizes were significant (p < 0.001), but weaker, ranging from 0.18 to 0.38 (Table 1). This indicates that larger seeds produce plants with larger cotyledons, but the strength of this correlation was lower than among seed traits. Cotyledon length and width were significantly (p < 0.001) correlated at 0.88, indicating that longer cotyledons were wider. Fresh and dry weights of young plants were significantly (p < 0.001) correlated at 0.88 (Table 1). Seed and cotyledon sizes, except for seed width, had significant (p < 0.001) but relatively low correlations (0.09 to 0.13 and 0.15 to 0.22, respectively) with fresh and dry weights of plants (Table 1). Table 1. Pearson's correlation coefficients between seed size (area, length, width, and perimeter), cotyledon (length and width), and early growth (fresh and dry weight 20 days after planting) traits. Because of the significant correlations among seed traits, cotyledon sizes, and plant weights, results from statistical analyses are presented for seed area, cotyledon length, and plant fresh weights; analyses of other traits are available from Olberg [63]. Genotype and experiment were significant (p < 0.001) sources of variation for these three traits ( Table 2). The genotype by experiment interaction was significant for cotyledon lengths (p < 0.001) and plant fresh weights (p < 0.001), but not for seed area. Although this interaction was significant, the genotype variance was greater than variance of the genotype-byexperiment interaction across all traits, indicating consistency in the ranking of genotypes across environments. The block within experiment (Block(E)) interaction was significant (p < 0.001) for seed area, cotyledon length, and fresh weights ( Table 2). Variation across experiments (greenhouses) may have been caused by different light levels or temperatures across the three greenhouse plantings which were temporally separated by one week. Variation in microclimates within greenhouses may have contributed to the significant block effects.
GCA, SCA, and reciprocal effects were all highly significant (p < 0.001) sources of variation for seed area, cotyledon length, and fresh weight (Table 2), indicating both additive and non-additive effects on these traits. GCA for seed area was significant for all DHs (Table 3). However, direction of deviation from zero was not consistent for DHs extracted from the same population, indicating that overall parental influence on seed area is DH specific. Both DHs from 9930 had negative GCAs for seed area, whereas DHs from TMG1 had positive GCAs. These indicate a general reduction and increase in seed areas for progeny of 9930 and TMG1, respectively. For all DHs except ST8-2, GCAs for cotyledon lengths were significant (p < 0.05) ( Table 3). For all DHs, GCA was not significant for fresh weights, indicating that no individual DH had a significant overall influence on early growth of progeny (Table 3). Shen et al. [43] reported greater GCA deviations for fresh and dry weight measurements across populations, and this may be due to a longer period of growth (22-30 days before sampling) as compared to the 20 days in this study.  SCA was significant for specific hybrids across all traits (Table 4). More hybrids had significant SCA for seed areas and cotyledon lengths compared to plant fresh weights. Hybrids of DHs extracted from 9930 and GY14 had highly significant positive SCA for fresh weights, indicating beneficial combining abilities of these populations. The hybrid of GY14-9 and ST8-4 had a significant negative SCA for fresh weight, while other GY14 and ST8 hybrids did not show this effect, indicating a DH-specific interaction. Hybrids produced from crossing of the two DHs extracted from the same source population had highly significant negative SCAs (Table 4), likely due to higher homozygosity and inbreeding depression for within-population hybrids. Although low levels of inbreeding depression and heterosis have been reported for yield components in cucumber, these findings are consistent with moderate levels of heterosis previously reported for growth and days to flowering in cucumber [43,64,65]. Table 4. Specific combining ability (SCA) for seed area, cotyledon length, and plant fresh weights expressed as deviation from overall mean of zero for progeny of two separately derived DHs from cucumber populations of 9930 (9930-3; 9930-5), GY14 (GY14-15; GY14-9), Straight 8 (ST8-2; ST8-4), and TMG1 (TMG1-4; TMG1-5) crossed in a full 8 × 8 diallel mating scheme. All DHs showed highly significant (p < 0.01) reciprocal effects for seed areas and cotyledon lengths (Table 5). Each DH displayed significant (p < 0.05) reciprocal effects for plant fresh weight in at least one hybrid combination ( Table 5). Direction of effect was consistent across DHs extracted from populations 9930 and ST8, which indicate specific maternal or paternal effects. For example, both DHs from 9930 had positive paternal effects for fresh weight in hybrids with ST8 and GY14 and may indicate a beneficial mitochondrial effect on early growth conferred by 9930. In agreement with results by Shen et al. [43], significant reciprocal effects were observed for fresh weight for hybrid TMG1xGY14, with greater fresh weights when TMG1 was used as the maternal parent. This combination may warrant further investigation into this specific positive interaction. In contrast, while Shen et al. [43] also found positive maternal effects on early growth for TMG1 in a cross with ST8, we found negative maternal effects for this hybrid combination (Table 5). These results indicate that while reciprocal effects are often significant contributors to early plant growth across hybrids, direction of effects are DH as opposed to population specific. The significances of SCA and reciprocal effects for fresh weights reveal the importance of nonadditive genetic effects on early plant growth. This is in contrast to findings that GCA and additive effects were the main contributor to early plant vigor in maize, though significant SCA and reciprocal effects also contributed [47].

RNA-Seq of Reciprocal Hybrids
RNA-seq was completed using two sets of reciprocal hybrids with significant reciprocal effects (hybrids of ST8-4 with 9930-3 and hybrids of ST8-2 with GY14-15) and two sets with non-significant reciprocal effects (hybrids of ST8-4 with GY14-15 and hybrids of ST8-4 with TMG1-5). Reads were aligned to the 9930 reference (ASM407v2) and across all samples there was an average raw count of 44.5 M reads, primary read alignment >90%, and an average of only 2.7% unaligned reads. Transcriptome analyses revealed only two differentially expressed genes (DEGs), Csa7G368160 and Csa5G623470, in common between the hybrid pairs with significant reciprocal effects (Figure 1), consistent with hybrid-specific reciprocal effects (Table 5). Gene Csa7G368160 was upregulated in both higher growth hybrids with a log2(fold change) of 4.51 and 6.86 between ST8-4x9930-3 versus 9930-3xST8-4 and ST8-2xGY14-15 versus GY14-15xST8-2, respectively. This was the highest fold change observed across both sets of reciprocal hybrids. Gene Csa7G368160 encoded an uncharacterized protein in cucumber and partial sequence homology with Arabidopsis thaliana revealed similarity to an mRNA for auxin-induced root culture-like protein (AIR12). AIR12 orthologs were identified as plasma membrane ascorbate-reducible b-type cytochromes in bean and soybean [66]. Further investigation into the function of this gene could reveal its potential role in reciprocal effects on early growth in cucumber. The other common DEG (Csa5G623470) encodes an MLO-like protein and was down-regulated in both higher growth hybrids. Down-regulation of this gene has also been associated with powdery mildew tolerance in cucumber [67].
Many of the DEGs were annotated as uncharacterized proteins, which may be due to relatively poor annotation of the cucumber reference sequence. Annotated gene descriptions in combination with PANTHER gene classification [62] of the DEG sets between reciprocal hybrids with significant reciprocal effects revealed a wide array of predicted functional categories (full list of DEG annotations and classifications available in Table A1 of Olberg [63]). Both sets of DEGs include chloroplast-and mitochondrial-related genes, indicating that organelle-nuclear interactions play a role in the significant reciprocal effects. Investigation into Gene Ontology (GO) overrepresentation in the DEG sets between hybrid pairs with significant reciprocal effects revealed no commonalities. DEGs between ST8-2xGY14-15 and GY14-15xST8-2 were significantly overrepresented (FDR < 0.05) for genes relating to enzyme regulation, and specifically protease and hydrolase activity (Table S1). Eliminating the genes that were differentially expressed between the control reciprocal hybrid pairs from this set of DEGs nullified all previously identified overrepresented GO terms for ST8-2xGY14-15 and GY14-15xST8-2 (Table S2). Therefore, it is unclear whether these DEGs are related to the observed reciprocal differences in early plant growth. The genes differentially expressed between hybrids ST8-4x9930-3 and 9930-3xST8-4 were significantly enriched (FDR < 0.05) for GO terms related to DNA-binding and transcriptional regulation, and were associated with the nuclear chromosomes and chromatin (Table S1). Genes related to these GO terms, within the molecular function and cellular compartment domains, remained significantly overrepresented when all genes differentially expressed between the control hybrids were removed from the DEG set (Table S2). This indicates that differences in DNA-binding transcription factor activity may be involved in mediating the reciprocal effects on early plant growth observed for ST8-4x9930-3 and 9930-3xST8-4. These differences in gene expression patterns between the two sets of reciprocal hybrids with significant reciprocal effects suggest separate pathways are involved in mediating the observed differences in early growth for reciprocal hybrids of ST8-2 with GY14-15 and ST8-4 with 9930-3.
Among genes differentially expressed between control hybrids ST8-4xTMG1-5 and TMG1-5xST8-4, genes related to protein folding were significantly (FDR < 0.05) overrepresented, while no GO terms were found to be significantly enriched between GY14-15xST8-4 and ST8-4xGY14-15 (Table S1). These gene expression patterns are not associated with differences in early growth, but indicate potential expression differences between reciprocal hybrids which could contribute to differences in other phenotypes across these reciprocal hybrids.
Kollipara et al. [68] reported many DEGs across various functional categories for maize reciprocal hybrids with a divergent response to cold germination and desiccation stress. It is possible that differential allelic contribution could influence early plant growth, though there is little evidence for varied allelic expression in reciprocal hybrids beyond the embryo stage. Springer and Stupar [69] observed significant maternal effects in embryo allelic expression in Arabidopsis reciprocal hybrids, but found no evidence of parent-oforigin effects on allelic expression in seedlings. Guo et al. [70] reported reduced yield and heterosis in hybrids with paternally biased allelic expression in maize, but these differences were not attributed to reciprocal differences. Variation in DNA methylation can contribute to asymmetric allelic contribution as well as variation in overall gene expression [71][72][73]. He et al. [69] reported differential gene expression patterns between rice reciprocal hybrids associated with variation in DNA methylation, but observed no differences in allelic contribution across reciprocal hybrids. Epigenetic variation between reciprocal hybrids may contribute to variation in gene expression and early growth, warranting future investigation.

Conclusions
The specificity of gene expression differences within the two reciprocal hybrid pairs of cucumber provides evidence for hybrid-specific parent-of-origin effects rather than common pathways across hybrids mediating significant reciprocal effects on early plant growth. This specificity is consistent with the significant SCA observed for the phenotypic traits. Further investigation into more sets of reciprocal hybrids with significant reciprocal effects could give a broader picture of the various pathways involved in these reciprocal differences. Understanding the mechanisms of varied interactions could reveal potential targets of selection for beneficial parent-of-origin effects. Future investigation of the methylome in these reciprocal hybrids could reveal altered DNA methylation profiles associated with changes in gene expression patterns. The hybrid combinations with significant reciprocal differences for early plant growth provide a unique system for studying epigenetic variation and its potential effects on heterosis and growth in cucumber.
Given the observed DH-specific nature of reciprocal effects, it is not surprising that varied pathways are likely involved in mediating differences in early plant growth across hybrid combinations. While selection of cucumber hybrids with larger seeds or cotyledons could confer some increase in early plant growth, the benefit of this selection is likely to be minor. Given the DH-specificity and variation in effect size and direction of reciprocal effects across populations and hybrid combinations, we recommend that breeders generate and evaluate both reciprocal hybrids to identify better performing hybrids.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/agronomy11101908/s1, Table S1: Gene Ontology (GO) terms enriched within the sets of significantly (FDR < 0.05) differentially expressed genes (DEGs) between reciprocal hybrid pairs with significant reciprocal effects (ST8-2/GY14-15 and ST8-4x/9930-3) and no reciprocal effects (controls; ST8-4/TMG1-5 and GY14-15/ST8-4) compared to number of genes represented in reference set (Ref) of all detected genes across samples in the experiment. Significance was tested based on Fisher's exact test with the false discovery rate (FDR) correction. Table S2: Gene Ontology (GO) terms overrepresented within sets of significantly (FDR < 0.05) differentially expressed genes (DEGs) that are exclusive to each reciprocal hybrid pair with significant reciprocal effects compared to number of genes represented in reference set (Ref) of all detected genes across all samples in the experiment. Significance was tested based on Fisher's exact test with the false discovery rate (FDR) correction.  Data Availability Statement: All data are available from the corresponding author; all sequences are available from https://www.ncbi.nlm.nih.gov/genbank as BioProject PRJNA765045.

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