Homozygosity mapping identified loci and candidate genes responsible for freezing tolerance in Camelina sativa

Homozygosity mapping is an effective tool for detecting genomic regions responsible for a given trait when the phenotype is controlled by a limited number of dominant or co‐dominant loci. Freezing tolerance is a major attribute in agricultural crops such as camelina. Previous studies indicated that freezing tolerance differences between a tolerant (Joelle) and susceptible (CO46) variety of camelina were controlled by a small number of dominant or co‐dominant genes. We performed whole genome homozygosity mapping to identify markers and candidate genes responsible for freezing tolerance difference between these two genotypes. A total of 28 F3 RILs were sequenced to ∼30× coverage, and parental lines were sequenced to >30–40× coverage with Pacific Biosciences high fidelity technology and 60× coverage using Illumina whole genome sequencing. Overall, about 126k homozygous single nucleotide polymorphism markers were identified that differentiate both parents. Moreover, 617 markers were also homozygous in F3 families fixed for freezing tolerance/susceptibility. All these markers mapped to two contigs forming a contiguous stretch of chromosome 11. The homozygosity mapping detected 9 homozygous blocks among the selected markers and 22 candidate genes with strong similarity to regions in or near the homozygous blocks. Two such genes were differentially expressed during cold acclimation in camelina. The largest block contained a cold‐regulated plant thionin and a putative rotamase cyclophilin 2 gene previously associated with freezing resistance in arabidopsis (Arabidopsis thaliana). The second largest block contains several cysteine‐rich RLK genes and a cold‐regulated receptor serine/threonine kinase gene. We hypothesize that one or more of these genes may be primarily responsible for freezing tolerance differences in camelina varieties.


INTRODUCTION
Camelina (Camelina sativa), also referred to as false flax, gold of pleasure, or linseed dodder is an under-used oilseed crop from the Brassicaceae family (Berti et al., 2016;Vollmann et al., 2007;Kagale et al., 2014). Camelina is an emerging biofuel crop (Moore et al., 2017) that has both winter and spring biotypes (Anderson et al., 2018). In more moderate regions, it can be planted as a winter annual, but the spring variety, which does not require vernalization, is generally preferred by growers in the Northern Great Plains because winter survival is inconsistent for most available commercial varieties (Gesch & Archer, 2013;Berti et al., 2016). Camelina has reportedly been cultivated in Europe for over 5000 years (Bouby, 1998) but has gained notoriety recently as a potential source of omega-3 fatty acids and for producing niche oil and animal feed stocks (Malik et al., 2018;McVay & Lamb, 2008;Luo et al., 2020). The major use of camelina oil is for a hydroprocessed renewable jet fuel (Berti et al., 2016;Moore et al., 2017). The US military revealed that camelina oil met all jet engine performance standards and reduced greenhouse gas emissions significantly (Shonnard et al., 2010). Like other winter Brassica oilseed crops, camelina sets seed early in the growing season and produces good yields (743-2300 kg ha −1 in west central Minnesota) under harsh dry and cold environments with low inputs-even in saline soil (Gesch, 2014). Combined, these attributes make it an excellent choice for cropping in the Northern Great Plains region (Gesch & Cermak, 2011;McVay & Lamb, 2008). In addition, camelina shows potential as a weed suppressing cover crop (Saucke & Ackermann, 2006) that reduces soil nutrient losses (Berti et al., 2016) and offers an early season nectar source for pollinators (Eberle et al., 2015). Camelina also has the advantage of a short growth cycle (85−100 days) crop (Abdullah et al., 2016;McVay & Lamb, 2008) and thus, has promise for dual cropping opportunities for North Dakota (ND) growers (Berti et al., 2016;Gesch & Archer, 2013). Despite having potential, relative to canola, there has been little research to improve biotic and abiotic stress tolerance (Berti et al., 2016). Camelina is allohexaploid (2n = 40) in nature and shows close syntenic relation to arabidopsis (Arabidopsis thaliana). A draft reference sequence of the camelina genome was first released in 2014 and it appears to have an ∼782 Mb genome (Kagale et al., 2014). Anderson et al. (2018) also released whole genome and transcriptome sequence data from a winter biotype "Joelle'" and summer biotype "CO46" at ∼25× coverage. The camelina chromosome scale reference genome reported a total of 89,418 annotated protein-coding genes (Kagale et al., 2014) and a transcriptome atlas with a coverage of 88% of the annotated genes (Kagale et al., 2016). This available reference genome, together with camelina's ability to be transformed using agrobacterium with a floral dip pro-

Core Ideas
• We refined a homozygosity mapping procedure for freezing tolerance loci in camelina. • We identified several regions on chromosome 11 of camelina associated with freezing tolerance. • Candidate genes include cysteine-rich receptorlike protein kinases and receptor serine/threonine kinases.
cedure (Liu et al., 2012;Lu & Kang, 2008) has made camelina a unique crop species for genomic studies and has opened the opportunity to introduce useful agronomic traits through genetic engineering (Abdullah et al., 2016;Anderson et al., 2019;Berti et al., 2016). The camelina genome is hypothesized to be composed of an allotetraploid sub-genome with seven chromosomes each and a diploid sub-genome with six chromosomes (Kagale et al., 2014). Although having reasonable synteny with arabidopsis, the presence of three copies of homologous genes (allohexaploid nature), highly undifferentiated polyploidy, and small fractionation bias within its genome has made camelina more complicated for studies compared to arabidopsis (Anderson et al., 2018;Horvath et al., 2019;Kagale et al., 2014). Camelina poses a low level of genetic diversity, and a majority of the domesticated cultivars are spring type, whereas most of its wild relatives are winter type (Berti et al., 2016;Chaudhary et al., 2020). Freezing stress significantly impairs plant growth and productivity and is often fatal in extreme conditions (Thomashow, 1990(Thomashow, , 1999. Freezing damage is initiated by extracellular ice crystal formation, which later causes irreversible damage to plants through extreme cellular dehydration and cellular dysfunction (Ding et al., 2019;Fiebelkorn et al., 2018;Pearce, 2001;Thomashow, 1990). The ability of plants to survive freezing is advantageous for winter cover crops (Wilke & Snapp, 2008) and unraveling the freezing mechanisms both at genetic and physiological level is a well-researched phenomenon (Soorni et al., 2017). Plants can reduce freezing damage if first subjected to a period of time at low but non-freezing temperatures (Thomashow, 1990) often referred to as cold acclimation. Camelina can survive exceptionally low temperatures (−23˚C for 8 h) following acclimation (Berti et al., 2016). Cold acclimation results in increased levels of protective proteins such as late embryogenesis abundant (LEA) proteins, anti-freeze proteins (AFPs), cold shock proteins (CSPs), and synthesizing hydrophilic compounds such as proline and soluble sugars (Ding et al., 2019;Guy et al., 1994;Thomashow, 1990). In arabidopsis, cold acclimation also induces a family of transcription factors known as C-REPEAT BINDING FACTORs (CBFs), which in turn, regulates a set of COLD-REGULATED (COR) genes (Thomashow, 1990). This signaling process is often referred to as the CBF regulon (Ding et al., 2019;Jaglo-Ottosen et al., 1998). However, a transcriptomic study on two contrasting camelina varieties Joelle (a freezing tolerant winter biotype) and CO46 (a freezing sensitive spring biotype) revealed that the activation of CBF regulon was significantly induced in both varieties after 8 weeks of cold acclimation, although there were differences in the level of induction (Anderson et al., 2022;Horvath et al., 2019;Wang et al., 2021). These studies also indicated several sets of transcription factors and other potential regulatory factors besides CBFs were associated with differential survival following freezing in camelina. Elucidating the freezing survival mechanisms in camelina may lead to increased agricultural production of camelina and potentially enhance winter hardiness for other significant brassica species such as canola through genetic engineering approaches with camelina genes (Berti et al., 2016;Soorni et al., 2017).
Homozygosity mapping, also known as autozygosity mapping, is considered an efficient gene mapping procedure to identify recessive disorders/traits in inbreed families (Gholipoorfeshkecheh et al., 2020;Seelow et al., 2012). The principal mechanism behind the homozygosity mapping is to trace candidate genes by identifying homozygous regions in the affected individuals and identify causal mutation for a disorder through genome wide screening using a high-density single nucleotide polymorphism (SNP) based platform (Gholipoorfeshkecheh et al., 2020). A schematic for how this works is shown in Figure 1. Although homozygosity mapping is similar to bulk segregant analysis Michelmore et al., 1991;Wang et al., 2021), it differs in that rather than just selecting two extremes of the population with no regard to genotype, only individuals known to be homozygous for the traits are selected and compared to individuals known to be heterozygous. Also, rather than just looking for statistical associations between markers and the trait of interest, the goal of homozygosity mapping is to identify contiguous stretches of DNA that are homozygous in all individuals that are fixed for the trait, and which are heterozygous in all individuals where the trait is segregating. Homozygosity mapping has been successfully applied in humans and livestock research, but there are limited applications of this technique in plant species. Studies on several human recessive diseases had revealed that the chances of identifying the candidate genes within a homozygous block is over 90% (Soorni et al., 2017). Since homozygosity mapping is quite dependent on the availability of a reference genome, current cost-effective genotyping technology has opened the opportunity to explore this technique in several plant species (Kumar et al., 2020).
The apparently limited genetic loci controlling the major differences in freezing tolerance between Joelle and CO46, and the availability of a reference genome for this species makes it an ideal target to attempt homozygosity mapping to identify the primary causal genes underlying the freezing tolerance differences. In this study, we attempted to identify a reliable set of molecular markers through whole genome sequencing of winter biotype Joelle and summer biotype CO46 and the limited number of derived F3 RIL families that appeared to be homozygous for the loci controlling freezing tolerance or susceptibility, respectively. Our objective was to identify homozygous blocks and to look for causal candidate genes within and near these loci. Resulting candidate gene(s) should indicate potential targets to improve freezing tolerance in camelina and its related species such as canola.

Plant materials and phenotyping of F3 RIL
A reciprocal cross between two parental cultivars with contrasting freezing tolerance: Joelle (freezing tolerant) and CO46 (freezing susceptible), generated and 255 F2 seeds. Thus, F3 RIL families were generated by selfing the F2 seeds. Seeds from all 255 F3 segregating RILs were grown in 500 mL pots containing Sun-shine mix #1 (Fisons Horticulture Inc., Bellevue, WA, USA) under greenhouse conditions (16 h photoperiod with supplemental halogen lighting at 22-25˚C) with daily watering and weekly fertilizer application to the 6-8 leaf stage. Additionally, 254 F7 RIL families derived by single seed decent from these F2 seeds were also similarly tested to confirm the homozygosity of the phenotypes from the fully tolerant and susceptible F3 RIL families (one line was lost during generation of the F7 RIL populations).
Plants were then cold acclimated for 6 weeks in temperature-controlled walk-in growth chambers at 5˚C under 12 h photoperiods augmented with full spectrum LED lighting with an average intensity of 4500 lux at plant level prior to freezing. Acclimated plants were then subjected to a freezing treatment consisting of a linear ramp up to 15˚C at noon followed by a linear ramp down in temperature to 5˚C at 6:00 p.m. with lights on. Lights were turned off and plants were ramped down to −15˚C at 4:00 a.m. and held at that temperature for 4 h. Lights were then turned on, and plants were rapidly ramped up to 0˚C over the next hour and then back to 15˚C by noon. Such conditions are designed to mimic a hard early morning frost in the fall or early spring. Plants were then transferred to the greenhouse and scored for survival after 2 weeks using a visual survival scale of 0-3 (0 = dead, 1 = more than 70% leaf damage but surviving meristematic tissue, 2 = 30%-70% leaf damage and surviving meristematic tissue, 3 = 0%-30% leaf damage and surviving meristematic tissue) (Figure 2). Based on the above-mentioned phenotypic scale, freezing tolerant families that were fixed for freezing tolerance/susceptibility and families that appeared to

Resistant parent SuscepƟble parent
F I G U R E 1 Schematic workflow for homozygosity mapping, homozygous block identification, and candidate gene identification for freezing tolerance. The red lines represent the chromosome of parent that was homozygous for the trait of interest (freezing tolerant) and the blue lines represent the chromosome of the other parent which was homozygous for the wild type trait (freezing sensitive). Chromosomes with red and blue segments represent crossovers that occurred in the generation of the F2 individuals. The gene(s) of interest must be within the region that was both homozygous in the families of the F3 recombinant inbred lines (RIL) which were fixed for that trait of interest, and which were heterozygous in the F3 RIL families that were segregating for the trait of interest (the region depicted between the two horizontal lines).

F I G U R E 2
Photos of freezing damage with scale used in phenotyping.
still be segregating for freezing tolerance/susceptibility were identified.

Experimental design
A complete randomized design (CRD) was used for each experimental run with 12 technical replicates for each of the 255 F3 families, and 6 technical replicates for the two parental lines as positive and negative controls. After completing one round of phenotyping, we selected all the promising lines that were seemingly tolerant, susceptible, and segregating. We performed phenotyping one additional time for those selected lines. Based on the phenotyping data of the F3 RIL population, we selected 6 families with high ratios of freezing The Plant Genome tolerance, 6 families with low ratios of freezing tolerance, and 16 families that appeared to be segregating for freezing tolerance for whole genome sequencing in the F3 families.

DNA extraction and library preparation
To perform whole genome sequencing, young leaf tissues from 12 individual plants of each of the 28 selected families, along with leaf tissue from 6 plants of each parental variety, were collected and separately pooled. All samples were frozen immediately in liquid nitrogen and kept at −80˚C prior to DNA extraction. Pooled samples were homogenized, and DNA was extracted from each pool using a standard CTAB extraction (Doyle & Doyle, 1990). DNA concentrations were quantified using a NanoDrop 2000 (Thermo Fisher Scientific, USA) and the DNA integrity was evaluated by electrophoresis on 1% agarose gels. The high-quality DNA was used to produce genomic libraries using the NEBNext Ultra II FS DNA Library prep Kit for sheared DNA (200-450 bp) Illumina sequencing (catalogue #E7805L) according to manufacturer's protocols with the following minor modifications. Ampure beads were used for initial size selection followed by size selection using Pippin Prep gels (CDF1510), 2% Dye-Free, range 100-600 bp, and the peak between 200 and 300 bp was selected. Additionally, Pacific Biosciences high fidelity (PacBio HiFi) was used to sequence the parental lines (∼40× coverage for CO46 and ∼28× coverage for Joelle). The resulting sequences from both parents were assembled independently using the program Hifiasm with default parameters. The resulting ∼5000 contigs from CO46 along with an additional ∼1000 contigs from Joelle that did not appear to have any matches to the CO46 contigs were combined to serve as a new reference genome set (supplemental file S1) for mapping of short read sequences. The program MUMmer4 (Marçais et al., 2018) was used with default parameters to align CO46 contigs to the published camelina reference genome and to identify those contigs unique to Joelle. Additionally, BlastN was used to map previously published camelina gene sequences to this new reference genome set.

Whole genome sequence read mapping and SNP calling
Tagged genomic cDNA libraries were pooled and sent to Novogene Genomic Service & Solutions, China for whole genome sequencing at ∼30× coverage with 150 base paired end reads using Illumina NovaSeq technology. Raw whole genome sequenced FastQC data was quality trimmed and contaminating primers were removed using BBMAP (Bushnell, 2014) sub-programs (supplemental file S2). Briefly, trimming was done using BBDuk maintaining a parameter of minlen = 70, qtrim = rl, trimq = 20, ktrim = r, k = 28, mink = 12, hdist = 1, ref = /bbmap/resources/adapters.fa. Trimmed and filtered sequenced reads were then aligned to the assembled camelina contigs (supplemental file S1). In addition to these 30 libraries, sequence data from a previous whole genome sequencing project that provided additional files with greater than 25× coverage of both parental lines (Anderson et al., 2018) were added to the downstream analyses. Mapping of fragments was done using Bowtie2 with default parameters (supplemental file S3). This reference-guided approach was used to generate a vcf file using BCFtools software. All the procedures were completed using python bash scripts on SSH Linux platform. The resulting vcf file was then filtered for calling high-quality SNPs using VCFtools software, where read counts smaller than 40 for either major or minor allele were discarded and max-missing = 0.5, maf = 0.05, and minimum DP of 10 were maintained. The raw reads have been deposited in the SRA database (accession number PRJNA880815).

Identifying markers and candidate genes
From the above-mentioned quality filtered SNP vcf file, we identified markers that were homozygous but divergent in at least three of the four parental lines in both the current sequencing run and those from a previous whole genome sequencing of the two parental lines (Anderson et al., 2018). This resulted in 125,899 markers from ∼2 million seemingly high-quality raw SNP markers from the 28 F3 families and both sets of sequence from parental lines. These markers were then extracted from the vcf files of the F3 RIL families. Markers that were homozygous in both independent sequencing runs for the parental lines, and which were homozygous in the corresponding tolerant/susceptible F3 families, and which were also heterozygous in the F3 families segregating for freezing tolerance were identified by a simple sorting algorithm in excel. Briefly, from our identified set of markers, we screened for homozygous blocks using block length with a threshold of 50 kb distance apart. We defined a sliding window size of 50 kb distance and quantified the number of markers within each window. By dividing the number of markers per window by the total number of markers in our study, we obtained the relative frequency of markers per window value. Those relative frequency of markers per window values were plotted to visualize peaks indicating homozygous block. We did not consider a peak as homozygous block if the peak forming windows had less than 10 markers or a relative frequency of markers designated as homozygosity score of less than 0.016. Genes present within or adjacent to the homozygous blocks were identified using the BlastN-created Note: Ratios shown are the number of freezing tolerant individuals (with survival scores greater than 2) to freezing susceptible individuals (survival scores less than 2) from each F3 RIL family.
map described above. We also examined our homozygous blocks with the syntenic regions of other brassica species for freezing tolerance/susceptibility by mapping sequences within our blocks to loci in other species by sequence similarities and gene order and location from various sources and genome browsers. The workflow of homozygosity mapping is depicted in Figure 1.

Phenotyping analysis to select F3 RILs for whole genome sequencing
Phenotyping all F3 RIL families for freezing tolerance indicated that there was considerable variation in the response to freezing conditions between experiments. Thus, multiple experimental runs were needed to identify F3 families that consistently showed limited segregation for freezing tolerant/susceptible phenotypes. This analysis eventually identified three non-segregating F3 families with freezing tolerance similar to the Joelle parent, and 54 families with weak freezing tolerance similar to the CO46 parent. There were 198 families that appeared to still be segregating for freezing tolerance (Table 1). We selected 28 F3 families for whole genome sequencing based on survival score for freezing tolerance. Among those selected F3 families, three had >90% freezing tolerant individuals and three had >75% freezing tolerant individuals in both experimental runs, six appeared to have a high percentage of freezing sensitive individuals (100% in both experimental runs), and the remainder appeared to be segregating for tolerance/susceptibility (Figure 3). A list of selected accessions is shown in supplemental file S4. However, analysis of the F7 RIL mapping population indicated that only two of the six selected lines were truly homozygous for freezing tolerance and only two maintained homozygosity for freezing susceptibility (data not shown), suggesting that only a subset of the 12 selected F3 families were likely fixed for freezing tolerance/susceptibility.

Analysis of sequenced data and discovery of homozygous markers
We obtained ∼2 million seemingly high-quality SNP markers from the 28 F3 families and both sets of sequence from parental lines. Among them 125,899 markers were homozygous in both parents and consistently differentiated the two parents in the earlier study by Anderson et al. (2018). Two families did not have sufficient sequencing depth and were dropped from the analysis. Additionally, several families were dropped because their phenotype was divergent from what was expected in the F7 generation or were lost during advancement to the F7 generation. This left five freezing tolerant lines (three were fully tolerant with a survival score of 3 and the other two were highly freezing tolerant, showing maximum survival score of close to 3), four freezing susceptible lines, and six lines that were likely still segregating for freezing tolerance. Additional filtering for high-quality SNPs using VCFtools software resulted in 57,251k SNP markers which were homozygous in both independent sequencing runs of the two parents (the current run and the previous run described by Anderson et al., 2018). Of these, only 617 markers had limited or no missing data and were homozygous in all families fixed for freezing tolerance/resistance and heterozygous in the families that were segregating for freezing tolerance. Intriguingly, these 617 markers mapped to just two contigs (contig 2 and contig 25) (supplemental file S5). Out of the 617 markers, contig 2 and contig 25 had 365 and 252 markers, respectively. Unexpectedly, many of the other markers on these contigs, although homozygous in the parents, still showed segregation in the four non-segregating progeny. Thus, we looked for blocks of homozygosity where closely linked markers (within a block length of 50,000 bases) were found in consecutive stretches of 10 or more markers (see yellow highlighted regions in supplemental file S5). However, a total of seven and two closely linked homozygous blocks were detected on contigs 2 and 25, respectively. These resolved into four distinct peaks on contig 2 and one peak on contig 25 (Figure 4). It should be noted that both contigs map contiguously on Chromosome 11 in the previously published camelina reference genome (Table 2).

Discovery of probable candidate genes for freezing tolerance or susceptibility
We identified 22 camelina genes that had homology to sequences within or very near (within 5 kb upstream or downstream) of our selected homozygous blocks (Table 2). Additionally, we identified two genes within or near our homozygous blocks that were differentially expressed (Table 2green highlighted) between Joelle and CO46 after cold T A B L E 2 Camelina genes that map to regions of contigs 2 and 25 which contain extended stretches of markers that are homozygous in both parent and offspring that are fixed for the trait, and which were heterozygous in offspring that appeared to still be segregating for the trait ID

Note:
Bolded genes indicate those that were shown to be differentially expressed between the freezing tolerant and freezing susceptible genotypes following cold acclimations in a previous study (Anderson et al., 2022).

F I G U R E 3
Boxplot showing phenotypic distribution and resistance level of selected 28 F3 RIL families and their parents (based on average survival score from two experiments) that were advanced for whole genome sequencing.

F I G U R E 4
Peak indicating homozygous block based of relative frequency of markers per window of 50 kb distance apart (generated from 617 markers).
acclimation in previous studies (Anderson et al., 2022). The largest peak was on contig 2 between the position of 9,518,950 and 9,698,140. Three contiguous blocks of 50 kb distances fall within this peak. Surprisingly, of the four peaks observed on contig 2, most were not enriched with genes previously reported to be involved in freezing tolerance. However, on contig 2, one block within the large peak contained four genes located between position 6,844,201 and 6,905,478 with two genes of interest. One of these genes is similar to a rotamase cyclophilin 2 that has been implicated in freezing tolerance in arabidopsis (Weng et al., 2020), and the other encodes a plant thionin gene that was differentially expressed following cold acclimation in camelina (Anderson et al., 2022). The peak on contig 25 contained two blocks with six genes (Table 2). Fortuitously, these six genes form two clusters of three cysteine-rich RLK (RECEPTOR-like protein kinase) genes and three cold-regulated receptor serine/threonine kinase gene each. However, we cannot rule out that other loci within or outside of these two contigs might also be playing a role in the differences in freezing tolerance observed in the parental lines.

DISCUSSION
Extreme winters and short summer growing seasons are a common phenomenon in Northern Great Plain regions where a relatively short growing season and good winter hardiness of winter camelina have made it an excellent choice for growing as a cover crop (Berti et al., 2017a;Gesch & Archer, 2013;Gesch & Cermak, 2011;Gesch et al., 2014Gesch et al., , 2018. However, in the Northern Great Plains, summer biotypes of camelina are predominantly grown because winter kill is still occasionally problematic (Berti et al., 2016). To obtain the maximum economic benefits, camelina of a double-or relay-cropping system, more consistent winter survival, and early flowering is required . This is only possible if the genotype has strong winter survival potential in Northern Great Plain regions. The winter genotype Joelle has shown excellent winter survival and overwintering at extremely low temperature in field conditions (Gesch & Cermak, 2011;Gesch et al., 2018;Walia et al., 2018). Thus, our study was designed to elucidate the freezing tolerance mechanism that differentiates it from a highly freezing sensitive spring genotype C046. Identifying the loci regulating these differences will provide a valuable target for future breeding to enhance freezing tolerance in camelina and potentially other related brassica crops such as canola.

Efficacy of the phenotyping in our study
We repeatedly phenotyped all the F3 RIL families to select homozygous and segregating progeny for freezing tolerance. Additionally, we repeatedly phenotyped 254 of the original 255 F7 families derived by single seed decent from these F3 families. Thus, we are confident that the four F3 families chosen for the homozygosity mapping are homozygous for the freezing tolerance/susceptible trait. Previous studies suggested that freezing tolerance between the parents of these families was controlled by only one or two dominant genes . This is clearly not the case. The discrepancy is most likely due the limited number of F2 progeny examined in the earlier study and the observation that there is considerable phenotypic plasticity in the response of lines that are segregating for freezing tolerance. The data here, with 9 lines out of 254 showing full parent-like phenotypes consistently from the F3 through the F7 generation suggest that there may be as few as three co-dominant genes that must be homozygous to give a parental phenotypic response to freezing stress.

Homozygosity mapping in other systems
Homozygosity mapping technology has previously been limited to human and animal studies. Several researchers used different approaches to find homozygous blocks and candidate genes with variable sized homozygous blocks/windows. In humans, homozygosity mapping identified regions and candidate genes for congenital heart disease (Gholipoorfeshkecheh et al., 2020), recessive genes for 13 different diseases (Hildebrandt et al. 2009), and candidate genes for the trait of hereditary hair loss (Zeinali et al., 2021). Although the size of each identified homozygous block varied significantly, in 93% of the cases the disease-causing candidate genes were found within the identified homozygous blocks. Limited applications of homozygosity mapping have been reported in plants. Kumar et al. (2020) applied homozygosity mapping in pear plants to reveal their population history and trait architecture using SNP markers. In another study, Machado et al. (2016) used microsatellite markers to detect homozygosity and genetic variability in an F4 castor bean population without performing a genome-wise scan.  used a modified homozygosity mapping approach to identify the gene underlying the barley variegation mutant luteostrians. However, to our knowledge, the application of homozygosity mapping to identify candidate genes for a nonmutant trait from a bi-parental segregating population has not been reported in any plants. We are the first attempting to use this relatively new technology to identify freezing tolerant genes in camelina.

High throughput sequencing ensures good quality markers for freezing tolerance or susceptibility
According to several studies, freezing tolerance in plants is a complex trait (Augustyniak et al., 2020;Huang et al., 2018;Wrucke et al., 2020). Phenotyping plants for complex traits such as freezing tolerance is difficult, highly labor intensive, and time consuming because one must screen multiple plants at a large scale (Huang et al., 2018). Thus, developing good quality molecular markers for freezing tolerance would be beneficial to breeders. A set of reliable molecular markers linked to freezing tolerance not only can reduce breeding time and space but also increase the efficacy of reliable selection (Huang et al., 2018). Our data identified multiple markers associated with several loci that appear to be linked to freezing tolerance in two camelina varieties. Additionally, we have identified over 126K markers that consistently differentiate the two parents used in our study. It is likely that these markers could also be used to identify other genes of interest controlling agronomically important traits that differ between these two varieties including flowering time, spring versus winter biotype, yield, seed size, and oil content/quality. However, the small size of the homozygous blocks in the F3 families that were fixed for freezing tolerance suggests that a fair number of the markers that differentiate the two parents may not be segregating as expected and should be further screened to identify those that show a probable Mendelian segregation pattern in a segregating population.

4.4
Homozygosity mapping reveals candidate region for detecting causal gene for freezing tolerance In our study, we identified 617 markers that were homozygous in the parental varieties and that were already fixed for freezing tolerance in the F3 families. The majority of these markers mapped to two contigs that represent a large portion of chromosome 11 based on homology to the published reference genome. The contigs in question are large (18,617,245 bp for contig 2 and 7,057,017 bp for contig 25) and appear to be adjacent and contiguous with less than a 2000 base gap between them as based on the published reference sequence. However, not all markers within these contigs showed homozygosity in both parents and the F3 families. It is unlikely that the median density of homozygous markers (1155 bases between markers) was too sparse to clearly delineate homozygous blocks as crossover rates on average should be roughly 200 kb per every 1% crossovers (1 centimorgan) (Goodman et al., 1995). It is more likely that the error rate in calling homozygous markers was too high to accurately generate the likely existence of large homozygous blocks. That said, we did identify multiple stretches of closely linked homozygous markers on these two contigs that ranged in size from a few hundred bases to more than 20 kb in size. It is noteworthy that among the 22 genes with strong homology to sequences within or near these homozygous blocks, each of the two largest peaks contained a gene that was differentially expressed between the two parents following 8 weeks of cold acclimation or were differentially expressed in one or both parents as they cold acclimated over time (Anderson et al., 2022). Thus, this observation would be consistent with these homozygous regions having a higher-than-average role in freezing tolerance.
Genes that make up the CBF regulon are often key targets when studying freezing tolerance (Ding et al., 2019;Thomashow, 1990). Indeed, several association mapping studies for freezing tolerance in arabidopsis have identified a region that include the three CBF genes on chromosome 4 of arabidopsis (Gehan et al., 2015;Horton et al., 2016;Oakley et al., 2014). Intriguingly, two CBF genes are located on chromosome 11 in camelina. Although these genes were not within or adjacent to any long stretches of closely linked homozygous markers, the fact that these critical genes are colocalized within the two contigs that host the majority of the homozygous markers means we cannot rule out the possibility that selection of strong freezing tolerance/susceptibility observed in the parental lines is not the result of genetic differences within these genes. However, although there is differential expression of the CBF genes between acclimated CO46 and Joelle, CBF is induced in both varieties by cold treatments (Anderson et al. 2022;Horvath et al., 2019).
We also examined other genes within or near homozygous blocks on chromosome 11 that were previously identified as having been associated with syntenic regions in other brassica species. Several quantitative trait loci (QTL) associated with freezing tolerance have been identified in arabidopsis on chromosomes 1, 4, and 5 (Horton et al. 2016). Intriguingly, the majority of the candidate genes identified in the arabidopsis study from both chromosome 4 and 5 mapped to syntenic regions on camelina chromosome 11 including CBF. However, the syntenic regions, although mapping relatively close to the homozygous blocks identified in our study, were not overlapping with them (within the 50 kb zones). Few QTL studies on freezing tolerance have been done in brassica crops where the QTL regions have been mapped to the reference genome. We found only one such study which identified a large region of Chromosome A08 of Brassica napus that was associated with differential freezing tolerance (Huang et al., 2018). Surprisingly, there are no large syntenic blocks between B. napus and camelina in chromosome A08. However, over 300 genes from chromosome 11 of camelina had their best matches to chromosome A08 of B. napus based on a BlastN analysis (data not shown), and there are genes within the large QTL on chromosome A08 that map to chromosome 11 of camelina. However, there were no obvious syntenic regions within the homozygous blocks on chromosome 11 that might suggest synteny of freezing tolerance genes between camelina and B. napus.

Possible candidate genes responsible for freezing tolerance/susceptibility
Our homozygosity mapping study detected 22 candidate genes associated with freezing tolerance mechanisms that were similar to sequences located within 5 kb upstream and downstream of homozygous blocks. Intriguingly, BlastN indicated multiple transcripts from chromosome 11 and 12 annotated as cysteine-rich RLK had their top match to the second highest peak that was located on contig 25. Additionally, BlastN indicated three receptor serine/threonine kinase-like protein transcripts that had their top match to sequences within this peak. Of the likely transcripts that mapped to this locus, one of the receptor serine/threonine kinase-like transcripts was differentially expressed between the two parents. Also, several closely related cysteine-rich RLK genes that mapped near this locus but not directly in it were also differentially expressed in previous studies (Anderson et al., 2022). Furthermore, this locus had the largest stretch of homozygous markers. Considering the above-mentioned facts, we hypothesize that this locus may be one of the primary loci responsible for the differences in freezing tolerance between the two parental varieties.
The only gene located in our targeted locus on contig 25 that was clearly differentially expressed between the two parents, encodes a receptor serine/threonine kinase-like protein and annotated as Csa11g022770. The most similar gene in arabidopsis is AT4G18250. According to the expression browser from The Arabidopsis Information Resource website (https://www.arabidopsis.org/), AT4G18250 is induced in the roots in response to salt stress and the leaves in response to light stress as well as following various biotic stress events. However, it was not shown as being induced by cold stress. Surprisingly, this gene, like the surrounding cysteine-rich RLK genes, was induced in response to cold acclimating conditions in CO46, but not in Joelle. Interestingly, the arabidopsis orthologue of the receptor serine/threonine kinase genes were previously found to play a role in ABA signaling and plant abiotic stress response (Lim et al., 2015;Osakabe et al., 2010). The overexpression of Serine/threonine kinase-based genes were reported in camelina (Anderson et al. 2022), in pepper (Tang et al., 2022), in arabidopsis (Mao et al., 2010), and in grapevine (Sawicki et al., 2019) during cold exposure/treatment. Intriguingly, modulating serine/threonine kinase-based genes have shown to enhance freezing tolerance in arabidopsis (Ding et al., 2015;Mao et al., 2010) and in pepper (Lim et al., 2020). These strengthen our finding and hypothesis that this serine/threonine kinase-based gene is probably involved in freezing tolerance and ABA signaling in camelina.
The other genes located within this peak encode cysteinerich RLKs. It should be noted though that the ∼50 kb homozygous region is too small to encode all of the genes that show high homology to this stretch of DNA. Indeed, an analysis of the sequence suggests it contains only one full-length copy of a cysteine-rich RLK gene that is highly similar to but not identical to the reference camelina transcript Csa11g022780 and several portions of other likely cysteine-rich RLK sequences which could be multiply spliced. The fact that previously predicted cysteine-rich RLK genes have their best matches to this region-many of which were preferentially expressed in the freezing sensitive CO46 parent, suggests that there may be some structural anomalies in the predicted location of these genes in the published camelina reference genome and possibly some expansion of the gene family resulting from mischaracterized splice variants of these genes. However, being a part of large RLKs family, cysteine-rich RLKs were reported to have influence in growth and development in arabidopsis and can modulate abscisic acid (ABA) response and drought tolerance (Lu et al., 2016;Pelagio-Flores et al., 2020;Zhang et al., 2013aZhang et al., , 2013b. ABA influences both drought and freezing tolerance through plant stomatal movement (Lu et al., 2016) and it has been reported that around 10% of genes that were expressed during drought induction were also induced during cold stress (Shinozaki et al., 2003). In addition, cysteine-rich RLK1 (CRK1) was shown to play a positive role in protecting plants from cold stress (Ye et al., 2017) and were also differentially expressed during cold treatment in tomato (Chen et al., 2015) and in Medicago sativa (Song et al., 2016). Expression of cysteine-rich RLK proteins such as CRK4, CRK5, CRK19, and CRK20 were found to cause a hypersensitive responseassociated programmed cell death in arabidopsis (Burdiak et al., 2015;Lu et al., 2016;Wrzaczek et al., 2010). Thus, the presence of CRK gene may enhance freezing tolerance by increasing ABA production or may cause death of plants by negatively regulating the ABA pathway.
Even though the cysteine-rich RLKs and receptor serine/threonine kinase genes from contig 25 seems to be the most reasonable candidate for freezing tolerance based on our finding and cited literature reviews, we cannot rule out other possible genes, particularly from the longest stretches of homozygous markers located in contig 2. Interestingly, majority of the seven homozygous blocks from contig 2 contained either a limited number of genes or genes with limited information or involvement with freezing tolerance trait (for instance, "Csa11g045260" and "Csa11g045270" genes appeared in one of our homozygous blocks but had no matched genes in arabidopsis). However, one block from contig 2 positioned between 9,648,098 and 9,698,140 contained a gene encoding a rotamase cyclophilin 2 (Csa11g052140). The rotamase cyclophilin 2 gene is analogous to arabidopsis AT3G56070 which contains a peptidyl-prolyl cis-trans isomerase domain. This gene was previously reported to increase cold tolerance in arabidopsis and regulate plant growth and development (Weng et al., 2020). Another interesting gene from the largest homozygous peak of contig 2, is Csa11g052230, which is analogous to AT1G66100 and encodes a plant thionin function. Although we did not find any literature that shows it is involvement in cold/freezing tolerance, it is involved in defense responses and showed differential expression between the two parents. Taking these together, our results suggest two genes from contig 25 with cysteine-rich RLKs and receptor serine/threonine kinase function and one gene from contig 2 with cyclophilin function may play a significant role in freezing tolerance in camelina.
Further QTL and fine mapping study on the freezing tolerance trait in camelina may provide validation of the involvement of our identified genes in freezing tolerance mechanism in camelina.
Although there is high similarity between the gene sequences from the reference genome and the assemblies used for this study, there are significant differences as well. For example, multiple annotated genes mapping to two different chromosomes in the previously published reference sequence had their best matches to the limited number of cysteine-rich RLK sequences and the receptor serine/threonine kinase gene located in the peak on contig 25. Such anomalies indicate the limitations of the existing available reference genome. Thus, our findings support the need for an improved camelina draft reference genome with maximum coverage that can deal with such situations efficiently. We expect to develop a new reference for the spring type CO46 and winter variety Joelle from our assemblies. The molecular markers developed for this study should prove invaluable for ordering the various contigs into chromosome-level scaffolds.

CONCLUSION
Application of relatively new but effective homozygosity mapping technology can provide new insights for genetic studies in plants. Using homozygosity mapping derived markers highlighted regions along chromosome 11 of camelina that appear to differentiate the freezing tolerance differences between the spring biotype CO46 and the winter biotype Joelle. Candidate genes associated with these markers identified homozygous blocks containing multiple cysteine-rich RLK genes and a receptor serine/threonine kinase gene that may be responsible for freezing tolerance differences in our F3 camelina population. Likewise, the differences in the rotamase cyclophilin 2 gene could also be responsible for the variation in freezing tolerance between Joelle and CO46. However, the causative nature of these genes still needs to be confirmed by QTL and fine mapping studies and engineering lines camelina that over-or under-express these genes. Additionally, we cannot rule out the possibility that other loci on portions of camelina chromosome 11 may also play a role in freezing tolerance/susceptibility. Thus, further studies are needed to test the hypothesis that these and other loci/genes might also be important in determining the relative freezing survival between these two divergent camelina varieties. Once we are certain about the roles these candidate genes play in freezing tolerance, it should be possible to manipulate those genes by transformation or CRISPR modification and increase freezing tolerance of other spring varieties of camelina. These candidate genes may also provide valuable resource to enhance freezing tolerance in other brassica species such as canola.

A C K N O W L E D G M E N T S
Funding for this work was provided by AFRI grant 2019-67014-29302 and USDA-ARS CRIS project # 3060-21220-033-000-D.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.

D I S C L A I M E R
Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture. USDA is an equal opportunity provider and employer.