The plants used in this study were provided by Dr. Fred Eickmeyer, ESKUSA GmbH, and originate from a wild-type population of highly heterozygous individuals which is due to the fact that T. koksaghyz is a facultative outcrossing, sexually reproducing species. Seeds were collected from a large set of T. koksaghyz plants grown under the greenhouse conditions described below after open pollination. Seeds were sown in rows in large buckets containing a mixture of 70% topsoil (Botanical Garden Münster, Germany) and 30% sand. Seeds started to germinate after 3 days and only seedlings that had emerged on days three or four were used further. The plants were grown in the greenhouse at 22–25°C (day) and 14–18°C (night) with a 16-h photoperiod (20 klx light intensity). A first set of plants was grown for transcriptomic analysis. We harvested the plants after 12 weeks and assigned them to three pools of 21–36 plants based on root morphology: thick roots of medium length (pool A), short and thin roots (pool B), or long and thin roots (pool C) (detailed classification in Supplementary Table S1). The roots were measured lengthwise using a ruler and the diameter was determined using digital calipers. Each root was then divided into two parts: proximal (upper root including lateral roots) and distal (root tip up to the first emerging lateral root). The parts were separately frozen in liquid nitrogen, lyophilized and ground to powder for transcriptome analysis.

We then cultivated a second set of wild-type T. koksaghyz plants as described above, but the pots were placed outside so that the plants were exposed to natural weather and day length. This second set of plants was used to validate first possible candidate genes gained from the transcriptomic analysis of the first set of wild-type plants. Quantitative real-time PCR (qRT-PCR) was used for the gene expression analysis of individual plants harvested from this second set. The roots were harvested after 12 weeks and measured as described above. We then processed whole roots without division into proximal and distal segments, using the same process of freezing, lyophilization and grinding as described above. The whole-root material was used to provide sufficient amounts of RNA for purification and cDNA synthesis followed by qRT-PCR.

For spatial gene expression analysis, T. koksaghyz plants were cultivated for 12 weeks as described above in the greenhouse. The plants were harvested individually and subsequently separated into flower, peduncle, leaf, root, and latex. The tissue material was then processed as described above.

Massive analysis of cDNA ends (MACE) libraries (one library per pool) were prepared and sequenced by GenXPro (Frankfurt, Germany) as previously described (Zawada et al., 2014; Bojahr et al., 2016). The sequencing data generated by MACE are deposited in the SRA under the bioproject accession PRJNA783530.

The natural rubber content of the root material used to prepare MACE libraries was quantified by accelerated solvent extraction as previously described (Benninghaus et al., 2020). The natural rubber content was determined gravimetrically by weighing the empty collection vials and those containing dried hexane extracts, then calculating the weight percentage as a proportion of root dry weight using the exact quantity of root material.

To increase mapping accuracy, low-quality bases (Phred score <20) were trimmed from the 3′ end of all reads using cutadapt v1.16 (Martin, 2011). Given the nature of the MACE technique, reads originate from 50 to 800 bp upstream of the 3′ ends of the transcripts, which often results in reads containing poly(A) tail sequences. Therefore, following the previous step, a de novo script was used to detect and remove poly(A) tail sequences, which would otherwise interfere with mapping. Finally, reads <25 bp were discarded using cutadapt to prevent nonspecific mapping. Quality-trimmed reads were mapped to the reference genome (Lin et al., 2017) using segemehl (Hoffmann et al., 2009) with split reads included.

Raw read counts were generated using bamutils v0.5.9 (Breese and Liu, 2013) with partially counted multiple hits, which represented the basis for differential expression analysis using the MA plot method with random sampling (MARS) in the DEGSeq R package v1.32.0 (Wang et al., 2009). Transcripts were annotated by using the NCBI BLAST suite v2.7.1 (Altschul et al., 1990) to screen the NCBI non-redundant proteins (O’Leary et al., 2016) and UniProtKB/SWISS-PROT (Bateman et al., 2017) databases. Transcripts with a normalized one-fold or greater differential expression (log-2-fold change ≥ 1 or ≤ −1, respectively) combined with a p-value less than 0.001 (p < 0.001) were considered significant and the corresponding genes were defined as differentially expressed genes (DEGs). Candidate genes for root length and diameter were selected according to following conditions: transcripts that were differentially expressed in pools A vs B and B vs C but not A vs C were considered as candidates for root length, whereas those differentially expressed in pools A vs B and A vs C but not B vs C were considered as candidates for root diameter.

Multi-dimensional scaling (MDS) plot analysis was performed using limma R package v3.46.0 (Ritchie et al., 2015). Venn diagrams were prepared using Vennplex (https://www.nia.nih.gov/research/labs/vennplex; Cai et al., 2013). Heatmaps showing differential gene expression were generated using heatmapper (http://www.heatmapper.ca; Babicki et al., 2016). Gene Ontology (GO) enrichment was analyzed using WEGO (http://wego.genomics.org.cn; Ye et al., 2018, 2006) with all transcripts in the corresponding root part used as the background. Glycosylphosphatidylinositol (GPI) anchors, protein domains and signal peptides were predicted using PredGPI (Pierleoni et al., 2008), CD search (Lu et al., 2020) and SignalP v5.0 (Almagro Armenteros et al., 2019), respectively. Phylogenetic trees were constructed using MEGA X (Kumar et al., 2018) and graphically rendered using iTOL v6.3.1 (Letunic and Bork, 2021).

Total RNA was extracted from T. koksaghyz root material using the innuPREP RNA Mini Kit (Analytik Jena, Jena, Germany) according to the manufacturer’s instructions, and cDNA was synthesized using PrimeScript RT Master Mix (TAKARA, Clontech, Sain-Germain-en-Laye, France) according to the manufacturer’s instructions.

Quantitative real-time PCR was carried out as previously described (Laibach et al., 2015) with slight modifications. The elongation factor 1α (TkEF1α) and ribosomal protein L27 (TkRP) reference genes were used for the normalization of gene expression as described before (Niephaus et al., 2019; Benninghaus et al., 2020). The qRT-PCR data were analyzed, and normalized expression levels were calculated using Bio-Rad CFX Manager v3.1 (Bio-Rad Laboratories, Hercules, CA, United States). All oligonucleotides used for qRT-PCR analysis are listed in Supplementary Table S2 along with corresponding primer efficiencies. For determination of primer efficiencies, a standard curve using a serial dilution of a representative cDNA sample was generated and subsequently analyzed using Bio-Rad CFX Manager v3.1. Efficiencies are calculated by the software as described by Pfaffl (2001) and Vandesompele et al. (2002), converted to percentage efficiency.

Data analysis using the two-tailed t-test, Kolmogorov-Smirnov normal distribution analysis, Grubbs’s test for the identification of outliers, and correlation studies (linear fit with calculation of Pearson’s correlation coefficient) was carried out using OriginPro 2019b (OriginLab Corporation, Northampton, MA, United States).

To identify genes that potentially affect T. koksaghyz root morphology, we generated a set of wild-type plants and subdivided it into three pools of plants with different root characteristics. The root morphology of the pools differed significantly in terms of primary root length, diameter directly below the leaf rosette, and the volume of the upper 5 cm, which was calculated by determining the root volume at two defined points (Figure 1, Supplementary Table S1). However, all three pools produced comparable amounts of natural rubber (Supplementary Table S1). Pool A comprised plants with medium to long roots that were also thick (large volume), representing the most desirable root phenotype for biomass accumulation as well as harvesting and handling. Pool B comprised plants with short and thin roots, whereas pool C comprised plants with long and thin roots.

We compared gene expression profiles between pools in two segments of the root. Each complete root was divided into the proximal segment (upper root including lateral roots) and distal segment (root tip up to the first emerging lateral root), resulting in a total of six pooled samples. T. koksaghyz is a sexually reproducing and self-incompatible species, so wild-type plants show a high degree of genetic diversity. By pooling plants with similar root morphology prior to transcriptome analysis, we aimed to reduce the background noise and emphasize transcripts that are differentially expressed between pools, thus capturing those representing genes related to the development of specific root characteristics.

We prepared and sequenced cDNA libraries from all six samples using MACE technology (Kahl et al., 2012). The six libraries contained between 4.65 × 10⁶ and 6.68 × 10⁶ single-end raw reads of 16–145 bp in length, which were further processed using a custom Python script to remove poly(A) tails. Following the removal of adapters and the trimming of low-quality bases, the libraries contained between 4.55 × 10⁶ and 6.48 × 10⁶ clean reads, 86.3–87.3% of which mapped to the reference genome (Lin et al., 2017), and 74.3–78.5% of which were uniquely matched (Table 1). We obtained 30,955 unique transcripts from the six libraries, which were then annotated using the NCBI and UniProt databases. To confirm sufficient reduction of background noise by pooling, we constructed an MDS plot, which showed the clear separation of all samples (Supplementary Figure S1).

Next, we identified the genes with significant differential expression between pools, defined as those with a log2 fold change (l2fc) ≥ 1 or ≤ –1 (p < 0.001) in any pairwise comparison. We then constructed a Venn diagram and used this to find genes related to either root volume or root length (Figure 2A).

Genes with significant differential expression when comparing pool A vs B and pool A vs C were considered as candidates influencing root volume (RV in Figure 2A) because the volume of the roots in pool A was significantly higher compared to both pool B and pool C. The initial list of candidates was reduced by excluding those also showing differential expression when comparing pool B vs C because the roots in these pools were similar in volume. By applying this last criterion, we arrived at a final set of RV candidate genes (schematic description in Supplementary Figure S2). Genes with significant differential expression when comparing pool A vs B and pool B vs C were considered as candidates influencing root length (RL in Figure 2A) because the roots in pools A and C were significantly longer than the roots in pool B. As above, we narrowed the initial list of candidates by excluding those also showing differential expression between pools A and C because the roots in these pools were similar in length. By applying this last criterion, we arrived at a final set of RL candidate genes (schematic description in Supplementary Figure S2). The Venn diagrams showed that DEGs defining the RV intersection were either upregulated or downregulated in both comparisons (A vs B and A vs C), but no counter-regulated genes were observed. In contrast, all DEGs in the RL intersection were counter-regulated (either upregulated in A vs B but downregulated in B vs C or vice versa).

We identified 573 DEGs for the upper root segment. We extracted 40 RV candidate genes from the 43 DEGs in the RV intersection by manually excluding genes showing differential expression when comparing pool B vs C with p > 0.001, because we considered these l2fc values as not statistically robust (Figure 2B). The corresponding heat map revealed that 33% of the RV candidate genes were expressed more strongly in roots with a higher volume (pool A) compared to those with a lower volume (pools B and C) as indicated by the positive l2fc values (Figure 3A). We also extracted 55 RL candidate genes from the 58 DEGs in the RL intersection (Figure 2B), 67% of which were expressed at higher levels in long roots (pools A and C) compared to short roots (pool B) (Figure 3C).

We identified 1,118 DEGs for the root tip, including 129 RV candidate genes from 146 DEGs assigned to the RV intersection and 56 RL candidate genes from 62 DEGs assigned to the RL intersection (Figure 2C). We found that 58% of the RV candidate genes were expressed at higher levels in roots with a larger volume (Figure 3B), and 55% of the RL candidate genes were expressed at higher levels in longer roots (Figure 3D). All DEGs and candidate genes are listed in Supplementary Table S3.

GO term enrichment analysis was carried out to determine the putative functions of the RV and RL candidate genes. Overall, we observed a similar distribution of second-level GO terms for all sets of candidate genes in the three main GO categories: cellular component, molecular function and biological process (Supplementary Figure S3). The most abundant GO terms assigned in the cellular component category were cell, cell part, organelle, organelle part, and membrane. For the molecular function category, the most represented GO terms were catalytic activity and binding. For the biological process category, the most represented GO terms were cellular process, metabolic process, response to stimulus, biological regulation, and regulation of biological process. Several GO terms were significantly enriched or underrepresented at level 2 (Supplementary Figure S3). Among the RV candidate genes in the root tip dataset, the cellular component extracellular region was enriched compared to all transcripts, whereas multiple terms (cell, cell part, membrane, and membrane part in the cellular component category; cellular process, metabolic process, biological regulation, and cellular component organization or biogenesis in the biological process category) were underrepresented compared to all transcripts. In addition to this general characterization, more specific GO terms at deeper levels were also evaluated by enrichment analysis (Table 2). Compared to all transcripts, hydrolase activity was enriched in the RV candidate genes of the upper root. Several GO terms describing response to different stimuli (chemical, hormone, oxygen-containing compound) were enriched in the RV candidate gene set in the root volume. Among others, intrinsic component of membrane, integral component of membrane, regulation of metabolic process, cellular metabolic process, and cellular macromolecule metabolic process were underrepresented in this candidate gene set. For the GO term distribution in the RL candidate genes of the upper root we found components not associated with membrane-bound organelles to be enriched (non-membrane-bounded organelle, intracellular non-membrane-bounded organelle) and vice versa (membrane-bounded organelle, intracellular membrane-bounded organelle). Cell periphery was enriched in the RL candidate genes of the root tip.

To characterize our transcriptome dataset further, we searched for typical factors affecting root growth and development like transcription factors, hormone pathway genes and cell wall/glucan biosynthesis/degradation genes based on the UniProt annotation (Supplementary Table S4–S6).

We identified 10 differentially expressed transcription factors which mostly originate from the volume and length trait groups in the root tip. Members of this set include two basic helix-loop-helix transcription factors (Transcription factor bHLH62, GWHPAAAA007503; Transcription factor bHLH145, GWHPAAAA042185), a transcription factor of the TCP family (Transcription factor TCP8, GWHPAAAA044352), an AP2 domain containing or RAP2 (related to AP2) family protein (Ethylene-responsive transcription factor RAP2-4, GWHPAAAA002086) and Protein LIGHT-DEPENDENT SHORT HYPOCOTYLS 10 (GWHPAAAA017340) of the ALOG family. Other members of the ALOG family were proposed to be developmental regulators because they promote cell growth in Arabidopsis thaliana (Zhao et al., 2004; Cho and Zambryski, 2011). This might indicate that the member we identified in our dataset plays a role during root development in T. koksagyhz.

We also found 16 genes associated with or induced by certain hormone pathways, including brassinosteroids (Protein EXORDIUM, GWHPAAAA001556 and GWHPAAAA007090; Shaggy-related protein kinase eta, GWHPAAAA008535), auxin (Auxin response factor 2A, GWHPAAAA014478; Remorin, GWHPAAAA014754) and jasmonates (Protein TIFY 10A, GWHPAAAA022842; Protein SGT1 homolog B, GWHPAAAA001300; Allene oxide synthase 3, GWHPAAAA036026). Most genes were found in the volume trait groups of the upper root and root tip. Protein EXORDIUM was shown to promote growth in the root and shoot when overexpressed in A. thaliana (Coll-Garcia et al., 2004) and Auxin response factor 2A downregulation increased root branching in tomato (Hao et al., 2015).

The last group comprised 11 genes encoding cell wall/glucan biosynthesis/degradation factors. Members of this group included not only glucan endo-1,3-β-D-glucosidases (β-1,3-glucanases; GWHPAAAA003778, GWHPAAAA042727, GWHPAAAA038049, GWHPAAAA027039) and xyloglucan endotransglucosylases (GWHPAAAA046431, GWHPAAAA046431, GWHPAAAA013192, GWHPAAAA013193), but also a cellulose synthase A catalytic subunit 6 (UDP-forming)-like gene (GWHPAAAA015656). The latter was previously found to play a crucial role in root and hypocotyl elongation in A. thaliana (Hauser et al., 1995; Fagard et al., 2000) and is also associated with the root length trait of the root tip in our dataset.

To confirm the suitability of our approach for root transcriptome analysis, we generated an independent set of wild-type T. koksaghyz plants (Supplementary Table S7). These plants were grown outside under natural climate conditions to approximate field growth. The roots were harvested individually after 12 weeks. Subsequently, qRT-PCR analysis was performed using primers specific for a small set of genes from the previous step, which were not obviously associated with root growth, as our goal was to identify novel factors influencing root morphology. They were also selected to represent various types of developmental and cellular mechanisms.

GWHPAAAA042727 was differentially expressed in the root tip and was annotated as a glucan endo-1,3-β-d-glucosidase (β-1,3-glucanase; UniProt ID Q94G86). This gene was strongly expressed in thin roots (pools B and C) but only weakly expressed in thick roots (pool A) and therefore represented a putative negative regulator of root volume. This was confirmed by qRT-PCR (Figure 4A). Additionally, spatial gene expression analysis showed the highest expression in roots (Supplementary Figure S4A). Protein domain analysis revealed a glycoside hydrolase family 17 (GH-17) domain spanning amino acids 23–343 and an X8/CBM43 superfamily domain spanning amino acids 372–456. The first part of the N-terminus corresponded to a signal peptide with a cleavage site between residues 21 and 22, but no GPI anchor site was predicted, suggesting the protein is not attached to the membrane but probably binds to carbohydrates through its X8/CBM43 domain (Barral et al., 2005). Phylogenetic analysis (Supplementary Figure S5) revealed that this β-1,3-glucanase belongs to the α-clade characterized by a diverse set of protein domain structures and expression patterns (Doxey et al., 2007). Like approximately one third of the β-1,3-glucanases in this clade, the T. koksaghyz protein represents protein domain architecture type II, consisting of a signal peptide, a GH-17 domain and an X8/CBM43 domain.

GWHPAAAA036026 was differentially expressed in the proximal root and was annotated as allene oxide synthase 3 (TkAOS3; UniProt ID Q8GZP5). This gene was also expressed at higher levels in pools B and C compared to pool A and thus represented another putative negative regulator of root volume. Amino acid sequence analysis showed that it shares the same characteristics as other CYP74 family P450 cytochromes, lacking the threonine residue required for the binding and activation of oxygen in the I-helix of P450 monooxygenases (Chapple, 1998). Instead, the small hydrophobic residue valine is found at this position, which is common for CYP74 proteins (Itoh et al., 2002). TkAOS3 also fulfils the CYP74 consensus sequence for cysteinyl heme ligand residues near the C-terminus with the sequence NKQC (A/P)(G/A)K(D/N)XV (Itoh et al., 2002). An unrooted maximum-likelihood tree (Supplementary Figure S6) containing several members of the CYP74 family from different species was constructed to assign TkAOS3 to one of the four CYP74 subfamilies. TkAOS3 clustered with other members of CYP74 subfamily C, which show high affinity for 9-hydroperoxides of linolenic or linoleic acid (9-HPOT/9-HPOD) derived from 9-lipoxygenase (9-LOX) activity, yielding the corresponding allene oxides (9,10-EOD/9,10-EOT) (Itoh et al., 2002; Stumpe et al., 2006). These unstable intermediates are likely to undergo spontaneous cyclization (forming 10-OPDA/10-OPEA) or hydrolysis (forming α-ketol/γ-ketol). However, the tomato AOS3 gene (LeAOS3) is thought to encode a multifunctional enzyme that can also catalyze the cyclization or hydrolysis of the intermediate to form rac-cis-10-OPEA or 9R-α-ketol, respectively (Grechkin et al., 2008). The expression of AOS3 is restricted to the roots in tomato (Itoh et al., 2002) and to the sprouting eyes, stolons, tubers and roots in potato (Stumpe et al., 2006). Our qRT-PCR analysis did not find any correlation because the expression of GWHPAAAA036026 did not negatively correlate with root volume (Figure 4B) as initially suggested by the MACE dataset. Among the different plant tissues, we observed the highest gene expression in latex, which is highly abundant in root tissue. In all other examined tissues, GWHPAAAA036026 is only weakly expressed in comparison to latex (Supplementary Figure S4B).

GWHPAAAA022842 was differentially expressed in the root tip and was annotated as a homolog of TIFY10A/JAZ1 (UniProt ID Q9LMA8). This gene was expressed at significantly higher levels in pools A and C compared to pool B, and was thus considered a candidate for a positive regulator of root length. Our qRT-PCR analysis in single plants revealed a positive correlation with the primary root length (Figure 4C), consistent with the MACE dataset. GWHPAAAA022842 is expressed in almost all tissues except flowers at a moderate level, and the highest expression was observed in the peduncle (Supplementary Figure S4C). The protein contains a TIFY domain (spanning residues 106–138) with a TIFLNG motif differing slightly from the consensus sequence [TIF (F/Y)XG] (Vanholme et al., 2007). It also features a Jas/CCT domain (spanning residues 182–207). A maximum likelihood tree was constructed from A. thaliana TIFY/JAZ proteins, revealing that our protein is closely related to AtTIFY10A and AtTIFY10B (Supplementary Figure S7). AtTIFY10A/JAZ1 belongs to the JAZ family of repressors that target transcription factors such as MYC2, MYC3 and MYC4 (Chini et al., 2007; Fernández-Calvo et al., 2011). These transcription factors promote the expression of several early jasmonate response genes and are involved in defense responses as well as developmental processes such as lateral root formation (Yadav et al., 2005; Grunewald et al., 2009; Kazan and Manners, 2013). The degradation of JAZ proteins depends on jasmonyl-isoleucine (JA-Ile)-mediated interactions with COI1 through its Jas domain, which is part of the SCF(COI1) E3 ubiquitin ligase complex, and subsequent processing by the 26S proteasome (Thines et al., 2007). The TIFY domain of JAZ1 interacts with NINJA/AFPH2, which recruits TOPLESS to repress MYC2 (Pauwels et al., 2010). In turn, the expression of multiple JAZ genes (including JAZ1) is promoted by MYC2 as well as the closely related transcription factors MYC3 and MYC4, resulting in a negative feedback loop to attenuate the jasmonate signal (Chini et al., 2007; Niu et al., 2011).