BSR-Seq Analysis Provides Insights into Cold Stress in Actinidia Arguta F1 Populations

The freezing injury, which is one of the important abiotic stresses in horticultural crops, can inuences the growth and development and the production area of kiwifruit (Actinidia Lind1). Actinidia arguta has excellent cold resistance in Actinidia species, but knowledge relevant to molecular mechanisms is still limited so far. Understanding the mechanism of cold resistance in kiwifruit is important for cold-resistant breeding. In our study, a cross of ‘Ruby-3’×‘Kuilv’ male was built for the A. arguta hardiness study, and 20 cold-tolerant and 20 cold-sensitive populations were selected from 492 populations according to LT50. Then, we performed Bulked segregant RNA-seq combined with single-molecule real-time sequencing to identify differentially expressed genes with cold hardiness. We found that the content of soluble sucrose and the activity of β-amylase were higher in the cold tolerant population pool than in the cold sensitive population pool. Upon -30°C low temperature treatment, 126 differentially expressed genes were found, and 59 genes were upregulated and 67 genes were downregulated when comparing the tolerant and sensitive pools, respectively. KEGG pathway analysis showed that the DEGs mainly belonged to starch and sucrose metabolism and amino sugar and nucleotide sugar metabolism. There were main 10 key enzyme encoding genes and two regulatory genes were up regulated in tolerant pool, regulated genes of CBF pathway were found to be different expressed, especially, 14-3-3 gene was down regulated and EBF gene was up regulated. To validate the BSR-Seq results, 24 DEGs were determined by qRT-PCR, and the results were consistent with BSR-Seq.


Abstract Background
The freezing injury, which is one of the important abiotic stresses in horticultural crops, can in uences the growth and development and the production area of kiwifruit (Actinidia Lind1). Actinidia arguta has excellent cold resistance in Actinidia species, but knowledge relevant to molecular mechanisms is still limited so far. Understanding the mechanism of cold resistance in kiwifruit is important for cold-resistant breeding.

Results
In our study, a cross of 'Ruby-3'×'Kuilv' male was built for the A. arguta hardiness study, and 20 coldtolerant and 20 cold-sensitive populations were selected from 492 populations according to LT50. Then, we performed Bulked segregant RNA-seq combined with single-molecule real-time sequencing to identify differentially expressed genes with cold hardiness. We found that the content of soluble sucrose and the activity of β-amylase were higher in the cold tolerant population pool than in the cold sensitive population pool. Upon -30°C low temperature treatment, 126 differentially expressed genes were found, and 59 genes were upregulated and 67 genes were downregulated when comparing the tolerant and sensitive pools, respectively. KEGG pathway analysis showed that the DEGs mainly belonged to starch and sucrose metabolism and amino sugar and nucleotide sugar metabolism. There were main 10 key enzyme encoding genes and two regulatory genes were up regulated in tolerant pool, regulated genes of CBF pathway were found to be different expressed, especially, 14-3-3 gene was down regulated and EBF gene was up regulated. To validate the BSR-Seq results, 24 DEGs were determined by qRT-PCR, and the results were consistent with BSR-Seq.

Conclusion
Our research provides valuable insights into the mechanism related to cold resistance in Actinidia and identi ed potential genes that are important for cold resistance in kiwifruit.

Background
Low temperature drastically in uences plant development, productivity and geographic distribution. In recent years, extremely low temperatures have occurred frequently. The kiwifruit industry suffers from an array of threats from low temperature stress [1]. Therefore, it is important to enhance cold resistance to minimize the economic loss from low temperature injury. Kiwifruit has been domesticated only in the past 100 years, and it has abundant wild resources, which contain excellent cold resistance traits, such as Actinidia arguta, which was found to withstand − 38℃ in our previous report (Sun et al., 2020). However, the lack of a comprehensive low temperature transcriptome, unexplored cold resistance genes and low temperature signaling hinder our full understanding of cold resistance in kiwifruit. Therefore, identifying cold resistance genes in A. arguta is a method for cold resistance breeding and improving kiwifruit cold resistance.
High plants have evolved elaborate mechanisms against cold stress. Cold stress is one of the most important mechanisms against low temperature stress in winter. The response of cold stress to low temperature is a highly complex process involving multiple levels of regulation [2]. A series of physiological and biochemical changes occurred during midwinter in plants [3]. These changes are involved in various pathways and ultimately increase freezing tolerance. The CBF/DREB1 pathway is a well-studied cold regulatory pathway that plays an important role in cold acclimation in Arabidopsis, an adaptive response where plants increase freezing tolerance after exposure to low nonfreezing temperatures [4][5][6]. Recent studies revealed that the CBF-dependent cold response involves transcriptional, posttranscriptional and posttranslational changes, expanding our knowledge of cold stress regulatory pathways [7]. The role of the starch metabolism pathway in plants under low temperature stress has been widely studied, and sugars accumulate rapidly in plants under low temperature [8]. The source of soluble sugar is generally thought to be from the metabolism of starch in plants [9].
The diploid 'Hongyang' was sequencing sequenced [10], the Kiwifruit Genome Database was built [11]; however, A. arguta is a tetraploid species, its polyploid nature and the incompleteness of its genome sequences and annotation limited the transcriptome analysis. Currently, single-molecular sequencing technologies provide an opportunity to thoroughly investigate the molecular mechanisms of the kiwifruit response to low temperature. Single-molecule real-time (SMRT) long-read sequencing technology from Paci c Biosciences (PacBio) is the most popular means of sequencing full-length (FL) cDNA molecules and has been used for whole-transcriptome pro ling [12,13]. FL transcript sequences that eliminate the need for assembly could provide direct information on the transcript isoforms of each gene [14]. SMRT sequencing has also been widely used to predict and validate gene models related to some unique traits.
However, the SMRT methodology cannot be directly used to quantify the expression level of transcripts, which might be corrected with next-generation sequencing (NGS) reads [15]. Bulked segregant analysis (BSA) can be used to identify markers linked to any speci c gene or genomic region using two bulked DNA pools. Each pool, or bulk, contains individuals who are identical in a particular trait or genomic region but arbitrary at all unlinked regions [16]. Bulked segregant RNA-seq (BSR-Seq) possesses the advantage of BSA and RNA-seq together, which has the full capability to identify differentially expressed genes (DEGs) and the ability to identify SNPs between different pools [17]. This method does not require genome information. The BSR-Seq method has been extensively applied to identify major genes in plants such as maize, ginkgoaceae, wheat and cabbage [18][19][20][21][22][23].
A. arguta possesses the strongest cold resistance, which might be involved in its genetic mechanism [24]. However, studies on this species are scarce. Therefore, studies of the genetic mechanism underlying the freezing tolerance trait of this species are still needed. In this study, we used PacBio Sequel and BSR-Seq to identify the DEGs in response to cold stress between tolerant and sensitive pools in A. arguta F1 populations. We identi ed some key genes in starch and sucrose metabolism and related regulated genes to this pathway. The development of these candidate genes will be the focus of future research, and these results will facilitate the study of the molecular mechanism of freezing tolerance in kiwifruit.

Low Temperature Treatment and Evaluation of Cold Resistance
Through the cross of 'Ruby-3' × 'Kuilv' male, a total of 492 populations were obtained, and all the shoots of populations were well planted(Additional le 6: Figure S1). When the dormancy shoots were treated at -30°C for 8 h, the REL of populations showed an approximately normal distribution (Fig. 1), and the REL ranged from 42% to 97% (Additional le 1: Table S1). Fifty populations with lower RELs (with an average REL of 49%) and 50 populations with higher RELs (with an average REL of 80%) were selected to calculate the LT50. The cold-resistant trait in the populations showed the phenomenon of the superparent. The detailed LT50 is shown in Table 1. Finally, 20 populations with the highest LT50 (tolerant pool) and 20 populations with the lowest LT50 (sensitive pool) were chosen for BSR-Seq analysis. The average LT50 of the 20 higher and lower cold-resistant populations were -30.52°C and -13.97°C, the highest LT50 was -36.90°C, and the lower LT50 was -7.51°C (Additional le 2: Table S2). β-amylase Activity and Total Soluble Sugar Content β-amylase activity and total soluble sugar content were measured in F1 populations in the tolerant pool and sensitive pool. The β-amylase activity was higher in the tolerant pool than in the sensitive pool, and the average β-amylase activity in the sensitive resistance and  Functional annotation of unigenes and analysis of DEGs A total of 28,496 unigenes were obtained for the following analysis. GO classi cation showed that most unigenes were associated with metabolic process, cellular process, single-organism process and biological regulation with the molecular functions of binding and catalytic activity. KEGG showed that the top clusters involved unigenes were signal transduction, carbohydrate metabolism, folding, sorting and degradation. Annotation against the NR database showed that unigenes in the PacBio transcriptome were identical to Vitis vinifera (25.6%), followed by Sesamum indicum (7.1%) and Juglans regia (7.0%), while the unigenes in the Illumina transcriptome were identical to Vitis vinifera (36.6%), followed by Sesamum indicum (6.8%) and Theobroma cacao (6.0%) (Additional le 7: Figure S2).
The differentially expressed genes (DEGs) between tolerant and sensitive pools were also determined.
After low temperature treatment, 126 genes displayed signi cantly differential expression between tolerant and sensitive pools, with 59 genes upregulated and 67 genes downregulated (Fig. 4).

GO and KEGG Pathway Enrichment Analysis
BLAST analysis and GO term annotation were performed to improve our understanding of the functions of these speci cally regulated genes. There were 17 GO terms were enriched in biological process, and 4 molecular function (starch synthase activity, transferase activity, transferring hexosy groups, glucosyltransferase activity and transferase activity)(Additional le 3: Table S3). These DEGs were signi cantly involved in KEGG pathways, including starch and sucrose metabolism and amino sugar and nucleotide sugar metabolism (Fig. 5). Seven unigenes related to starch and sucrose metabolism were upregulated by low temperature treatment. There were AGPase, granule-bound starch synthase, sucrose synthase (SUS), 1,4-alpha-glucan-branch enzyme (GBE), alpha-1,4 glucan phosphorylase, beta-amylase (BAM), glucan water dikinase (GWD), and neutral-alpha-glucosidase and disproportionating enzyme 2 (Additional le 4: Table S4).

Con rmation of Differentially Expressed Genes by qRT-PCR
To verify the reliability of the cold responsive gene expression pro les for DEGs, 27 DEGs that contained 18 upregulated and 6 downregulated were analyzed by quantitative real-time PCR(Additional le 5: Table  S5). The tolerant pool, sensitive pool and the three populations for random selection were used as templates. As shown in Fig. 6, the fold change values obtained by qRT-PCR were highly consistent with those based on BSR-Seq data for all of the selected cold responsive genes, despite the difference in the absolute fold change between the two methods. Therefore, some alleles originating from the tolerant pool were preferentially induced to be expressed under low temperature.

Discussion
A. arguta is a deciduous fruit tree, and it is a specie that has the higher cold resistance than other Actinidia species. A study on the cold resistance of A. arguta was signi cant for understanding the mechanism of cold resistance in Actinidia. Transcriptome analysis has been widely used in studies of kiwifruit, including investigations of fruit development and ripening [25,26], fruit color [27], and biotic and abiotic stresses such as waterlogging stress [28] and psa [29]. However, the materials in most studies were cultivars; in this study, F1 populations were used as the materials for the rst time. We performed the combined transcriptome analysis of NGS and SMRT sequencing and investigated the mechanism in response to low temperature.
Proline-rich proteins (PRPs) was found to be upregulated in the tolerant pool of populations, and some evidence suggests that PRPs is responsible for cell wall structure, such as GhHyPRP4, which may be involved in the plant response to cold stress in cotton [30]; Brassica BnPRP genes could be induced by cold [31]; and the Arabidopsis HyPRP gene protects the cells during freezing stress [32]. In this study, the higher expression of PRP genes, leading to proline accumulation, may be because the increase in proline added mechanical strength to the cell wall and stabilized the structure of organs under low temperature stress [33,34].
The pathways associated with starch and sugar metabolism were signi cantly enriched. Cold treatment also seemed to trigger enzymes responsible for the production of amylose, starch, maltose, and dextrin [35]. A few key changes in gene expression suggested that these pathways were utilized differently between F 1 populations with different cold tolerances. Starch is synthesized by starch synthase (SS) and 1,4-alpha-glucan-branch enzyme (GBE); then, starch is degraded to maltose by β-amylase (BAM), maltose transfers from plastid to cytoplasm, and maltose is catalyzed by disproportionating enzyme (DPE2), in which one glucose residue is released to the cytoplasm [36]. In our study, SS, GBE, BAM and DPE2 were all upregulated (Fig. 7). In particular, β-amylase plays an important role in abiotic stress [37]. The β-amylase gene family is induced by low temperature, and BAM1 and BAM 3 can degrade starch in Arabidopsis [38]. In trifoliate orange, PtrBAM1 is induced by cold but repressed by maltose; it is a member of the CBF regulon and plays an important role in cold tolerance by modulating the levels of soluble sugars acting as osmolytes or antioxidants [39]. In our results, the activity of β-amylase was higher in the tolerant pool, and the bam gene expression in populations was largely consistent with the activation of the enzyme, suggesting that both the expression and activity of β-amylase are upregulated by low temperature, which was in agreement with an earlier report [40].
SUS was upregulated in this study. Sucrose transforms the end product D-glucose. Soluble sugars such as sucrose, glucose, fructose, and ra nose have been shown to increase in concentration when plant tissues are in cold stress. These sugars are thought to provide freeze protection through stabilization of membranes, scavenging of reactive oxygen species, acting as signaling molecules, and decreasing the freezing point as compatible osmolytes [41]. We detected soluble sugars in F 1 populations. The sugar content was higher in the tolerance pool than in the sensitive pool, implying that tolerant populations produced more soluble sugars to resist cold and sugar reserves between nonsoluble and soluble sugar forms when faced with cold.
In our study, DnaJ and HSP70 were upregulated in the tolerant pool. Tomato DnaJ protein LeCDJ1 acted as an essential molecular chaperone in protein homeostasis and protein complex stabilization under stress conditions. Heat-shock protein 70 was identi ed as the partner of LeCDJ1 [42]. This indicates that DnaJ and HSP70 have essential functions in A. arguta against cold.
Plants have evolved sophisticated strategies in which hormone and cold signaling pathways are coordinated to better adapt to freezing stress. Ethylene insensitive 3 (EIN3) is a key transcription factor involved in ethylene signaling that negatively regulates freezing tolerance in Arabidopsis, partially by affecting the expression of CBFs by binding to the EBS motifs in their promoters [35,43]. Two F-box proteins, EIN3-binding F-box 1/2 (EBF1/2), positively regulate CBF expression by mediating the degradation of EIN3 and PIF3 via the 26S proteasome pathway [44,45]. In this study, EBF1/2 was upregulated in the tolerance pool, implying that A. arguta could increase cold resistance by increasing the expression of EBF to regulate the CBF gene. Cold perception seemed to differ in the sensitive pool and tolerance pool. 14-3-3 was upregulated in the sensitive pool. 14-3-3 family proteins are extensively involved in plant freezing tolerance. Under cold stress, CRPK1 phosphorylates 14-3-3 proteins and triggers their translocation from the cytoplasm to the nucleus; 14-3-3 proteins interact with and destabilize CBF1/3 proteins, thus attenuating the CBF pathway and preventing an excessive cold response [46]. Another study showed that Arabidopsis RCI1A, a 14-3-3 protein, negatively regulates freezing tolerance [18]. In our study, in the sensitive pool, the higher expression of 14-3-3 and the phenotype of sensitive cold resistance led us to speculate that 14-3-3 negatively regulates freezing tolerance in A. arguta, possibly by destabilizing the CBF pathway to stop COR genes (Fig. 8). In addition, we found that the transcription factor CHY was upregulated in the sensitive pool, CHY could respond to low temperature in pineapple, and A. arguta may try to resist cold in other pathways when the CBF pathway is suppressed.
Although some regulators of CBF were detected, we also failed to detect CBF transcription factors, just as the same results were found in Poncirus trifoliata [47]. First, most of the genes coding for the regulatory proteins, such as the CBFs, exhibited a transient and temporary change under abiotic stresses. At our low temperature (with 8 h), the mRNA abundance of the regulatory genes may be indistinguishable. Second, woody plants differ from herbal plants (such as Arabidopsis) when faced with low stress due to their natural variation in the length of the growth cycle [48]. Third, there were nonsigni cant differences in the tolerance pool and sensitive pool under a long time low temperature treatment. The identi cation of DEGs may be in uenced by several factors, such as stress intensity, duration of stress imposition, experimentation strategy and plant species/genotypes [49], and our results also support these conclusions.

Conclusions
In summary, the biparental F1 populations were constructed. Then, we performed the BSR-Seq combined with SMRT sequencing in F1 populations. DEGs induced by -30℃ low temperature treatment were identi ed. Low temperature treatment increased the content of total soluble sugar and the activity of βamylase in cold tolerant pool. 10 key enzyme encoding genes and 2 regulatory genes were upregulated, encoding genes including AGPase, granule-bound starch synthase, sucrose synthase (SUS), 1,4-alphaglucan-branch enzyme (GBE), alpha-1,4 glucan phosphorylase, beta-amylase (BAM), glucan water dikinase (GWD), neutral-alpha-glucosidase, disproportionating enzyme 2 and Proline rich protein (PRP). Regulatory genes including EIN-bingding F box (EBF) and 14-3-3. Most of the DEGs were enriched in starch and sugar metabolism and the regulated pathway. We concluded that starch degradation and soluble sugar synthesis is important for kiwifruit to resistance to cold.

Measurement of Electrolyte Leakage
After low temperature treatment, the shoots without buds were cut into 1 mm-2 mm thick slices. Then, 0.2 g of the slices were incubated in 30 ml of double-distilled water for 2 h with shaking at 200 rpm at room temperature. The initial electrical conductivity (C1) was measured using a digital conductivity meter (DDS-307, Rex, China). The samples were heated to boiling for 30 min and then cooled down to room temperature with shaking for 30 mins. The second electrical conductivity (C2) was taken. The REL was calculated as indicated by Eq. 1: The FT was expressed as LT50 (half lethal temperature at which REL reaches 50%) by tting the response curve obtained by the REL with a logistic sigmoid function (Eq. 2): y=k/(1+ae -bx ) (2) where x is the treatment temperature, y is the REL value, k indicates the extreme value when x approaches in nity, and a and b are the equation parameters. If the correlation coe cient r is close to 1, the equation is used to calculate LT50 (Sun et al., 2020).

Measurement of β-amylase Activity and Total Soluble Sugar Content
For measurement of β-amylase activity and total soluble sugar contents, the assay was repeated at least three times, with three biological replicates for each sample, and β-amylase activity and total soluble sugar content were determined using the relevant kits (Nanjing Jiancheng Bioengineering Institute) according to the manufacturer's instructions.

RNA-seq and Bulked Segregant Analysis
Tolerance pools and sensitive pools were generated for BSR-Seq analysis. Two extreme pools were constructed using equal amounts of tissues from the same position of plants at the same growing stage.
Shoots collected from the most cold-tolerant individuals and the most cold-sensitive individuals were mixed as tolerant and sensitive pools, respectively. Transcriptome sequencing was performed on the Illumina HiSeq 2500 platform (China Novogene, Beijing). The Illumina library was prepared according to the protocol described in the next-generation transcriptome article [50].

Gene Functional Annotation and DEG Analysis
Seven databases (NR, Swiss-Prot, GO, NT, KOG, Pfam, KEGG) were selected to map the nonredundant transcript sequences and obtain the annotation information of the transcript with e-values of 1e -5 against a total of seven databases [51][52][53][54][55][56].
RNA sequencing reads were aligned to the abovementioned SMRT transcripts obtained by Trinity using [57], and expression levels were estimated using RSEM [58]. The expression levels of the unigenes were expressed as fragments per kilobase of transcript per million mapped reads (FPKM) values to eliminate the in uences of gene length and sequencing quality difference on the estimated gene expression [59].
Differential gene expression analysis was carried out between the tolerance pool and sensitive pool using DEGSeq [60]. For the samples without biological repetitions, TMM was used to standardize readcount data [61], and then DEGSeq was used for differential gene expression analysis. In this study, genes with adjusted q values less than 0.005 and fold changes greater than or equal to 1 were identi ed as DEGs. Gene ontology (GO) term annotation was conducted for functional classi cation of DEGs [62], and further metabolic pathway enrichment analysis was carried out using KOBAS (2.0) against the KEGG (Kyoto Encyclopedia of Genes and Genomes) database [63,64].

qRT-PCR Analysis
Shoot tissues of F1 populations were analyzed in three biological replicates. Samples were ground as described in the Materials. RNA samples were converted to single-stranded cDNA using the cDNA Synthesis Kit (TOYOBO). Candidate genes were selected from previous data. Primer Premier 5 design primers were used (Table S1). SYBR Green-based RT-qPCR analyses were performed in a LightCycler 480 (Roche480) on a 96-well plate. The conditions for the PCR ampli cations were as follows: 95 ℃ for 5 min, followed by 45 cycles of 10 s at 95 ℃, 20 s at 60 ℃, and 20 s at 72 ℃. At the end of each experiment, a melt-curve analysis was carried out using the default parameters (5 s, 95 ℃ and 1 min, 65 ℃). β-actin in the kiwifruit was considered the control gene for normalization. All analyses were repeated three times using biological replicates. The relative expression levels were calculated using the 2 -ΔCt method. RT-qPCR primers of target genes were designed using Primer Premier 5 software. Availabiilty of data and materials The raw reads of the BSR-Seq data in this study have been deposited in the NCBI SRA database under accession number of PRJNA663731. The activity of beta-amylase and the content of soluble sugar in shoots of populations. A: the activity of beta-amylase, B: the content of soluble sugar