text:
Metagenomic analysis of cucumber rhizospheric microbiome provides insights into NaHCO3 stress mitigation by Rheinheimera pacifica NYJ

 Background: NaHCO3 in soil leads to alkalinization and causes stress in plants, thus hindering agricultural sustainability. A strain NYJ was identified as Rheinheimera pacifica and applied into cucumber-planted soil to mitigate NaHCO3 stress in this study. We aimed to detect the taxonomy and functional potential of rhizospheric microbiome affected by NaHCO3 stress and NYJ application, obtaining NaHCO3-tolerant microbes and exploring one mechanism underlining NaHCO3 stress mitigation by NYJ application.Results: The strain NYJ grew well under 0.01-0.06 g L-1 NaHCO3 and had plant growth-promoting capacity. Application of NYJ into NaHCO3-contaminated and cucumber-planted soil improved plant growth, increased relative water contents and osmotic potential of leaves, and decreased malonaldehyde and hydrogen peroxide contents in seedlings, mitigating NaHCO3 stress in cucumber. In cucumber rhizospheric soil, Proteobacteria, Actinobacteria, and Bacteroidetes of microbes were dominant phyla with 69, 17, and 6% of total abundance. To respond to NaHCO3 stress, 669 species and 208 genera of microbes were enriched in rhizospheric soil, while 190 sodium ion-related genes, 1242 osmotic stress-related genes, and 121 other stress-related genes from microbial genera including Pseudomonas showed the increased abundances. Under NaHCO3 stress, 26 microbial genera were distributed in connectors, the interaction network of rhizospheric microbes was changed, and soil enzyme activities were decreased. When NYJ being applied into the NaHCO3-contaminated soil, 1327 species and 487 genera of microbes were enriched than those under NaHCO3 stress alone. Meantime, soil enzyme activities were increased, while 2238 symbiont- and immune-related genes and 2434 plant growth-promoting genes including 16 genes encoding chitinase had the elevated abundances in microbial genera. Moreover, NYJ negatively interacted with one bacterial genus and positively interacted with 3 microbial genera under NaHCO3 stress. As a result, 25 bacterial genera were distributed in connectors, and microbial interaction network in rhizospheric soil was different from the one caused by NaHCO3 stress alone.Conclusions: 669 NaHCO3-tolerant species and 208 NaHCO3-tolerant genera of microbes were identified. NYJ application affected microbial community structure and interaction network in rhizospheric soil, increased the abundances of genes related to symbiont, immune, and plant growth-promoting traits in microbes, and activated soil enzymes, thereby mitigating NaHCO3 stress in cucumber.1 These authors contributed equally to this work.

4 rhizospheric soil, increased the abundances of genes related to symbiont, immune, and plant growthpromoting traits in microbes, and activated soil enzymes, thereby mitigating NaHCO3 stress in cucumber. 1 These authors contributed equally to this work.

Background
Microbes in plant rhizosphere interact with each other and also with plants [1]. These interactions play important roles in many vital processes such as nutrient cycling [2] and have impacts on plant growth [3]. Therefore, rhizospheric microbiome has received an increasing attention in recent years [4]. In the rhizospheric microbiome, plant growth-promoting rhizobacteria (PGPR) are beneficial to plants [5]. When PGPR are isolated from plant rhizosphere and then applied into planted soil, they not only change the levels of phytohormones, antioxidants, and epoxypolysaccharide [6], but also influence the rhizospheric microbial communities [7] by having competitive skills over the existing rhizosphere microbes [8], thus improving plant resistance to environmental stresses. For example, inoculation with Bacillus amyloliquefaciens NBRISN13 stimulates the population of betaine-, sucrose-, trehalose-, and glutamine-utilizing bacteria in rhizosphere and induces salt tolerance in rice [9].
Application of Streptomyces canus GLY-P2 has positive impact on microbial community composition in rhizosphere and mitigates stresses of ferulic and p-hydroxybenzoic acids in cucumber [10].
Soil salinization and alkalinization hinder plant growth and lead to desertification [11]. Alkaline soils, which are mainly caused by NaHCO 3 or Na 2 CO 3 , have an area of 434 million hectares in the world [12]. So NaHCO 3 stress has become one of the major constraints to crop yield and agricultural sustainability [13]. Under the stress of NaHCO 3

, inoculation of Enterobacter aerogenes LJL-5 and
Pseudomonas aeruginosa LJL-13 increases alfalfa yield and quality [14]. Application of Trichoderma asperellum decreases the infection incidence of pathogenic fungi and alleviates NaHCO 3 stress in poplar [15]. To the best of our knowledge, it has not been reported that PGPR would mitigate the stress of NaHCO 3 in plants by changing the microbial community composition and structure in plant rhizosphere. Moreover, uncovering the functional information of microbiome might contribute to 5 further understanding the PGPR-mediated stress mitigation [6], however, the function potentials of genes have not been reported in the rhizospheric microbes when PGRR being applied under NaHCO 3 stress. In addition, salt affects soil microorganisms, and many salt tolerance bacteria have been isolated [16]. Up to now, effects of NaHCO 3 stress on rhizospheric microbiome have not been analyzed, and NaHCO 3 tolerance microbes have not been identified.
In this study, a strain NYJ was isolated from rhizospheric soil, identified as Rheinheimera pacifica, and then applied into cucumber-planted soil under NaHCO 3 stress. It was hypothesized that the strain NYJ might belong to PGPR and mitigate NaHCO 3 stress in cucumber by affecting rhizospheric microbiome.
We aimed to determine the changes in taxonomy and functional potential of rhizospheric microbiome caused by NaHCO 3 stress and NYJ application, thus obtaining NaHCO 3 tolerance microbes and exploring one mechanism underlining the NaHCO 3 stress mitigation by NYJ application.

Materials And Methods
Isolation, identification, and properties of the strain NYJ Using a spread plate technique [17], the strain NYJ was isolated from the rhizosphere of wheat seedlings that have been cultured in saline-alkaline soils for 6 months. Each dilution was added to 20 mL of beef extract-tryptone agar medium (containing 10 g L − 1 tryptone, 5 g L − 1 NaCl, 5 g L − 1 beef extract, and 1.5% agar), in which the pH value was adjusted to 10 with NaOH. To identify the strain NYJ, its 16S rRNA gene sequence was amplified with the primers 1492R (5'-GGTTACCTTGTTACGACTT-3') and 27F (5'-AGAGTTTGATCCTGGCTCAG-3') and used to perform sequence alignment and phylogenetic analysis [18]. In order to analyze the ability of NYJ to tolerate NaHCO 3 , 1 mL of strain cultures grown exponentially in Luria-Bertani (LB) liquid medium were inoculated into 100 mL of beef extract-tryptone liquid medium supplemented with 0.01, 0.03, 0.04, 0.06, 0.08, or 0.1 g L − 1 NaHCO 3 . Subsequently, they were incubated at 37 °C and 180 rpm. The alteration of the optical density was observed at 600 nm [19]. To verify the NaHCO 3 tolerance ability of NYJ, 1 mL of strain cultures grown exponentially in LB medium were inoculated into 100 mL of the beef extract-tryptone liquid medium, in which the pH value was adjusted to 6, 6.5, 7, 8, 9, 10, 10.5, or 6 11. To determine the plant growth-promoting properties of NYJ and assess whether the strain would belong to PGPR, its ability of indole acetic acid (IAA) production [20], phosphate solubilization [21], and siderophore production [22] were detected, while the activity of 1-aminocyclopropane-1carboxylic acid (ACC) deaminase of the strain were analyzed [23].

Nahco Stress Mitigation In Cucumber By The Strain Nyj
According to the procedures of Zhang et al. [18], each cucumber seedling was cultured in one 10-cm plastic pot filled with 400 g of soil. The soil contained 1.37 g kg − 1 total N, 29.95 mg kg − 1 NO 3 -N, 13.87 mg kg − 1 NH 4 -N, 20.96 mg kg − 1 Olsen-P, 86.76 mg kg − 1 rapidly-available K, and 0.75% organic matter, and its pH value was 7.34. At the two-leaf stage, 32 cucumber seedlings were chosen and divided into four groups (8 plants per group) to impose the following treatments: I (control group), watered with 10 mL of autoclaved water; II (NaHCO 3 treatment group), watered with 10 mL of NaHCO 3 (0.12 g mL − 1 ) to make the concentration of NaHCO 3 reach 3 mg g − 1 dry soil, which induced stress in cucumber but did not cause plant wilting; III (NYJ treatment group), watered with 10 mL of NYJ (10 9 cfu mL − 1 ); IV (NYJ + NaHCO 3 treatment group), watered with 10 mL of solutions containing NYJ (10 9 cfu mL − 1 ) and NaHCO 3 (0.12 g mL − 1 ). Then, all of the seedlings were cultured according to Wu et al. [10]. After 12 days of treatment, the second leaves and rhizospheric soil of cucumber were collected.
In order to confirm the ability of NYJ to mitigate NaHCO 3 stress in cucumber seedlings, we measured the effects of the strain on cucumber seed germination under NaHCO 3 conditions. When being sterilized with 10% sodium hypochlorite, 60 cucumber seeds were put into the cultures of NYJ grown exponentially in LB liquid medium. Subsequently, the seeds were incubated at 28 °C and 180 rpm for 12 h. Meantime, another 60 sterilized seeds were incubated in LB liquid medium and were not treated with NYJ. After that, 30 NYJ-treated and 30 NYJ-untreated seeds were placed on sterile filter paper moistened with 3% NaHCO 3 in petri dishes. They were named as the "NYJ + NaHCO 3 treatment" and "NaHCO 3 treatment", respectively. Meantime, the remaining 30 NYJ-treated and 30 NYJ-untreated seeds were put on sterile filter paper moistened with sterilized distilled water, and they were named as the "NYJ treatment" and "control", respectively. Then, all of the filter paper was kept moistened with sterilized distilled water everyday. At 7 d, cucumber roots were collected.
To assess the effects of NYJ application on NaHCO 3 -stressed cucumber, the levels of plant growth parameters such as the fresh and dry weights of shoots and roots were determined [24]. Meantime, the activities of soil enzymes including catalase (CAT) [25], alkaline phosphatase [26], urease [27], and sucrase [28] were assayed in rhizospheric soil. The levels of relative water content [29], osmotic potential [30], and malonaldehyde (MDA) [31] in leaves were also measured. Since hydrogen peroxide (H 2 O 2 ) will accumulate in stressed plants and result in the over-production of MDA, the content of H 2 O 2 was assayed in leaves [32] to confirm the MDA level. In order to verify the H 2 O 2 content, the activities of H 2 O 2 -related antioxidant enzymes such as CAT [33], ascorbate peroxidase (APX) [34], guaiacol peroxidase (GPX) [35], dehydroascorbate reductase (DHAR) [36], monodehydroascorbate reductase (MDHAR) [37], glutathione reductase (GR) [38], and superoxide dismutase (SOD) [39] were determined in leaves. Since reduced glutathione (GSH) is regenerated by GR, its content were detected [40] to verify the activity of GR.
16s Rrna Amplicon Sequencing Of Microorganisms In Rhizospheric Soil Microorganism DNA in cucumber rhizospheric soil was extracted by using an EZNA® Soil DNA Kit (Omega Bio-tek, Norcross, GA, USA). Then, 16S rRNA amplicon sequencing and data processing were performed at Personal Biotechnology Co., Ltd. (Shanghai, China) [41].

Metagenome Sequencing Of Rhizospheric Microbes
When being extracted from rhizospheric soil, microbial genomic DNA was used to construct metagenome shotgun sequencing libraries by using the Illumina TruSeq Nano DNA LT Library Preparation Kit. Subsequently, each library was sequenced, and all data were processed at Personal Biotechnology Co., Ltd. (Shanghai, China) [42].

Results
NaHCO 3 tolerance ability and plant growth-promoting properties of the strain NYJ The 16S rRNA gene sequence of the strain NYJ showed 100% homology to Rheinheimera pacifica 8 (Additional file 1: Fig. S1), so the strain was identified as Rheinheimera pacifica NYJ. At 10 h, the cell numbers of NYJ under 0.01-0.06 g L − 1 NaHCO 3 were the same as the ones under 0 g L − 1 NaHCO 3 and were higher than those under 0.08-0.1 g L − 1 NaHCO 3 (Additional file 1: Fig. S2a). These were verified by the results that the strain NYJ had the same numbers of cells under pH 6.5-10.5 at 20 h and showed low numbers of cells under pH 6 and 11 (Additional file 1: Fig. S2b). Therefore, NYJ could grow under NaHCO 3 environment. On the other hand, NYJ had the capacity of IAA production, phosphate solubilization, and siderophore production and showed the activity of ACC deaminase, belonging to PGPR (Additional file 1: Fig. S3).

Nahco Stress Mitigation In Cucumber By Nyj Application
In comparison to the control, the shoot fresh and dry weights, plant height, areas of the third leaves, relative water contents and osmotic potential of the second leaves were significantly increased in cucumber seedlings of the NYJ treatment and were significantly decreased in the NaHCO 3 treatment, while the contents of MDA and H 2 O 2 in the second leaves were significantly decreased in the NYJ treatment and were significantly increased in the NaHCO 3 treatment ( Fig. 1a and Additional file 1: Fig.   S4). When compared the NYJ + NaHCO 3 treatment to the NaHCO 3 treatment, the shoot fresh and dry weights, plant height, areas of the third leaves, relative water contents and osmotic potential of the second leaves were significantly increased, and the contents of MDA and H 2 O 2 in the second leaves were significantly decreased. These were consistent with the results that the activities of APX, CAT, GPX, MDHAR, DHAR, GR, and SOD and the content of GSH in the second leaves were significantly increased in cucumber seedlings of the NYJ + NaHCO 3 treatment in comparison to the NaHCO 3 treatment (Additional file 1: Fig. S5). So inoculation with NYJ into soil mitigated NaHCO 3 stress in cucumber seedlings. Similarly, NYJ application increased root length, decreased the contents of MDA and H 2 O 2 , and elevated the activities of APX, CAT, DHAR, MDHAR, and SOD in roots, alleviating NaHCO 3 stress in germinating cucumber seeds (Additional file 1: Fig. S6).

Effects of NYJ application and NaHCO 3 stress on soil enzyme activity
Compared to the control, the activities of urease (Fig. 1b), phosphatase (Fig. 1c), and sucrase ( Fig. 1d) were significantly increased in the rhizospheric soil of the NYJ treatment and were significantly decreased in the NaHCO 3 treatment, while the activity of CAT (Fig. 1e) did not change in the NYJ treatment and were significantly reduced in the NaHCO 3 treatment. In the NYJ + NaHCO 3 treatment in comparison to the NaHCO 3 treatment, the activities of urease, phosphatase, sucrase, and CAT were significantly increased in the rhizospheric soil. Therefore, NYJ application increased soil enzyme activity under NaHCO 3 stress.
Effects of NYJ application and NaHCO 3 stress on microbial richness and diversity in rhizospheric soil Whether using the shotgun metagenomic sequencing or the 16S rRNA amplicon sequencing, the operational taxonomic units (OTUs) of observed species in the treatments of control, NYJ, NaHCO 3 , and NYJ + NaHCO 3 tended to be flat with the increasing sequence numbers (Additional file 1: Fig. S7).
Thus, the sequencing depth was reasonable to detect the majority of species in each treatment. The microbes in the rhizospheric soils of all treatments were identified ( Fig. 2 and Additional file 1: Fig. S8-S12). Based on the shotgun metagenomic sequencing, there were 7880 species shared in all treatments (Additional file 1: Fig. S13). Considerable overlap was observed in the NYJ-treated rhizospheric soils, with the NYJ and NYJ + NaHCO 3 treatments having 8404 shared species. Meantime, the rhizospheric soils treated with NaHCO 3 had considerable overlap, with the NaHCO 3 and NYJ + NaHCO 3 treatments having 8268 shared species. Moreover, in the control, NYJ, NaHCO 3 , and NYJ + NaHCO 3 treatments, the numbers of species were 8749, 8843, 8679, and 8772, respectively, and the distinct species were 64, 108, 97, and 74, respectively. Similarly, using the 16S rRNA amplicon sequencing, we not only observed the shared OTUs of microbial communities among four treatments, but also found the distinct ones in each treatment (Additional file 1: Fig. S14). According to the principal component analysis (PCA) based on the abundances of species (Additional file 1: Fig. S15) and genera (Additional file 1: Fig. S16), the microbial communities of the NaHCO 3 treatment separated from those of the control or NYJ + NaHCO 3 treatment. So soil microbial richness and diversity were changed by NaHCO 3 treatment and affected by NYJ application under NaHCO 3 stresses. 10 Effects of NYJ application and NaHCO 3 stress on microbial community composition and structure in rhizospheric soil In rhizospheric soil of the control, NYJ, NaHCO 3 , and NYJ + NaHCO 3 treatments, the majority of microbes were bacteria (Additional file 2: Table S1). Based on the shotgun metagenomic sequencing, Proteobacteria, Actinobacteria, and Bacteroidetes in rhizospheric soil were the dominant phyla with 69, 17, and 6% of total abundance, respectively ( Fig. 2a and Additional file 1: Fig. S17). The three phyla were also dominant in all treatments according to the 16S rRNA amplicon sequencing ( Fig. 2b and Additional file 1: Fig. S18).  (14), Acinetobacter (13), and Sphingomonas (13). Using the 16S rRNA amplicon sequencing, the genera Pseudomonas, Bacillus, Rhizobium, Microbacterium, and Acinetobacter were also found to respond to NYJ application (Additional file 1: Fig. S20 and Additional file 3: Table S3). Therefore, NYJ application affected microbial community composition and structure in rhizospheric soil.
To respond to NaHCO 3 stress, 834 OTUs were enriched and 1489 OTUs were depleted according to the shotgun metagenomic sequencing (Additional file 3: Table S4) Sphingomonas (14). Based on the 16S rRNA amplicon sequencing, the genera Pseudomonas, Rhizobium, Bradyrhizobium, Mycobacterium, and Sphingomonas were still in response to NaHCO 3 tress (Additional file 4: Table S5). Moreover, 208 genera were enriched and 421 genera were depleted in the NaHCO 3 treatment in comparison to the control by using the shotgun metagenomic sequencing (Additional file 4: Table S6), indicating that NaHCO 3 stress changed microbial community composition and structure.
Using the shotgun metagenomic sequencing, 1539 enriched and 612 depleted OTUs were found when compared the NYJ + NaHCO 3 treatment to the NaHCO 3 treatment, among which 1327 enriched and 489 depleted species were identified (Additional file 5: Table S7) Xanthomonas (12). Using the 16S rRNA amplicon sequencing, the genera Pseudomonas, Paenibacillus, Sphingomonas, Bradyrhizobium, Bacillus, and Legionella were also in response to NYJ application under NaHCO 3 stress (Additional file 6: Table S8). As shown in Additional file 6: Table S9, 487 genera were enriched and 137 genera were depleted in the NYJ + NaHCO 3 treatment when compared to the NaHCO 3 treatment. There were 257 depleted and 20 enriched genera in the NaHCO 3 treatment in comparison to the control, and these genera were enriched when compared the NYJ + NaHCO 3 treatment to the NaHCO 3 treatment (Additional file 7: Table S10). Meantime, 6 depleted and 78 enriched genera were found in the NaHCO 3 treatment when compared to the control, and they were depleted in the NYJ + NaHCO 3 treatment in comparison to the NaHCO 3 treatment. So NYJ application affected microbial community composition and structure under NaHCO 3 stress.

Effects of NYJ application and NaHCO 3 stress on the interactions of microbial communities
The interaction networks for rhizospheric microbes of the control, NYJ, NaHCO 3 , and NYJ + NaHCO 3 treatments were different from each other whether using the shotgun metagenomic sequencing or the 16S rRNA amplicon sequencing (Figs. 4 and 5, and Additional file 1: Fig. S21-S26). Based on the shotgun metagenomic sequencing, the strain NYJ negatively interacted with 10 species in the rhizospheric soil of the NYJ treatment (Additional file 1: Fig. S24), and it negatively interacted with 1 species and positively interacted with 11 species in the NYJ + NaHCO 3 treatment (Fig. 5). Meantime, using the 16S rRNA gene sequences, one OTU was identified as Rheinheimera pacifica, and its relative abundance was significantly increased in the NYJ treatment in comparison to the control or in the NYJ + NaHCO 3 treatment when compared to the NaHCO 3 treatment (Additional file 1: Fig. S27). The OTU of Rheinheimera negatively interacted with one genus in the NYJ treatment (Additional file 1: Fig.   S25), while it negatively interacted with one genus and positively interacted with 3 genera in the NYJ + NaHCO 3 treatment (Additional file 1: Fig. S26). So NYJ application and NaHCO 3 stress affected the interaction of microbes in rhizospheric soil.
The majority of the OTUs in the control, NYJ, NaHCO 3 , and NYJ + NaHCO 3 treatments were peripherals with most of their links being inside modules (Additional file 1: Fig. S28-S31). One network hub, which highly linked to many species in its own module and highly connected to several modules, was assigned to Acidovorax whether in the treatments of control, NYJ, NaHCO 3 , or NYJ + NaHCO 3 . There were three module hubs, which highly linked to many species in their own modules and respectively belonged to Pseudoxanthomonas, Flavobacterium, and Glutamicibacter in all treatments. In the connectors that highly connected to several modules, there were 18, 23, 26, and 25 identified genera in the control, NYJ, NaHCO 3 , and NYJ + NaHCO 3 treatments, respectively.

Effects of NYJ application and NaHCO 3 stress on the functional profiles of microbial communities
Based on the shotgun metagenomic sequencing, microbial communities in rhizospheric soil responded to abiotic and biotic stimulus and showed plant growth-promoting traits including nitrogen utilization (Additional file 1: Fig. S32). The genes related to osmotic stress were found in microbial communities (Fig. 6). When compared the NaHCO 3 treatment to the control, the relative abundances of osmotic stress-related genes were significantly decreased in GO:0006970, GO:0009651, GO:0071470, and GO:0030104 and did not change in GO:0061392, GO:0047484, GO:0006972, and GO:0071474 (Additional file 8: Table S11). However, 1242 osmotic stress-related genes showed significantly increased relative abundances in the NaHCO 3 treatment in comparison to the control. In 2 major functional genera for GO:0047484, 15 major functional genera for GO:0006970, 4 major functional genera for GO:0009651, 4 major functional genera for GO:0071470, 1 major functional genus for GO:0006972, 2 major genera for GO:0030104, and 1 major functional genus for GO:0071474, these abundance-increased genes were 22, 805, 165, 128, 35, 75, and 12, respectively. Among these, the osmotic stress-related genes of Pseudomonas had significantly enhanced relative abundances to respond to NaHCO 3 Table S12). There were 190 sodium ion-related genes, which showed significantly enhanced relative abundances in the NaHCO 3 treatment in comparison to the control (Additional file 8: Table S13). In 3 major functional genera (Pseudomonas, Bacillus, and Pseudoxanthomonas) for GO:0005343, 1 major functional genus (Pseudomonas) for GO:0005295, and 2 major functional genera (Pseudomonas and Pseudoxanthomonas) for GO:0005283, these abundance-increased genes were 96, 21, and 73, respectively.
In the microbial communities of rhizospheric soil, the genes implicated in other stresses except for the osmotic stress were found (Additional file 8: Table S14). When compared the NaHCO 3 treatment to the control, the relative abundances of other stress-related genes were significantly increased in GO:0043618, GO:0032101, GO:0032103, GO:0043619, GO:0034644, and GO:0009750. We observed 121 stress-related genes, which had significantly increased relative abundances in the NaHCO 3 treatment in comparison to the control (Additional file 8: Table S15). In 1 functional genus (Pseudomonas) for GO:0032101 and GO:0032103 and 2 major functional genera (Pseudomonas and Xanthomonas) for GO:0034644, these abundance-increased genes were 54, 54, and 13, respectively.

Discussion
Salt changes soil microbial community [16], and salt-tolerant microbes are screened by using a metagenomic library of saline soil [45]. In this study, 834 OTUs of microbes were enriched to respond to NaHCO 3 stress by using the shotgun metagenomic sequencing, among which 669 enriched species were identified. Meantime, there were 208 enriched genera by using the shotgun metagenomic sequencing when compared the NaHCO 3 treatment to the control in the current experiment. Although the 16S rRNA amplicon sequencing will show several artifacts [46], enriched genera were also found to respond to NaHCO 3 stress by using this approach. These results indicated that NaHCO 3 stress altered microbial community structure in rhizospheric soil, and 669 NaHCO 3 -tolerant species and 208 NaHCO 3 -tolerant genera were obtained for the first time in this paper. On the other hand, salt tolerance genes are mined from the metagenomic profiling of soil microbes [45]. In this study, 190 sodium ion-related genes, 1242 osmotic stress-related genes, and 121 other stress-related genes from the functional genera including Pseudomonas showed the increased relative abundances in the NaHCO 3 treatment in comparison to the control, suggesting that these sodium ion-, osmotic stress-, and other stress-related genes might play roles in tolerating NaHCO 3 stress. This could be verified by the result that the relative abundance of Pseudomonas was increased when compared the NaHCO 3 treatment to the control in this paper and would be confirmed by the report that Pseudomonas mendocina is tolerant to osmotic stress [47]. Moreover, elevated atmospheric concentrations of CO 2 change network of soil microbial communities [48]. In rhizospheric soil of this study, the interaction network of microbial communities and the OTUs of connectors had difference between the NaHCO 3 treatment and the control, indicating that NaHCO 3 stress altered the soil microbial interaction network and the connectors might contribute to the alteration.
Hou et al. [40] isolated a drought-tolerant strain of Bacillus, and they found that application of the strain into cucumber-planted soil improves plant growth, increases the relative water contents and osmotic potential of leaves, and decreases the MDA and H 2 O 2 contents of seedlings, thus mitigating drought stress in cucumber. In this study, inoculation with NaHCO 3 -tolerant NYJ into NaHCO 3supplemented and cucumber-planted soil had the same effects on plant growth, relative water contents, osmotic potential, MDA contents, and H 2 O 2 levels of leaves and thereby alleviated NaHCO 3 stress in cucumber seedlings. Application of Acinetobacter calcoaceticus CSY-P13 into soil mitigates the stresses of ferulic and p-hydroxybenzoic acids by altering soil bacterial community [7]. When compared the NYJ + NaHCO 3 treatment to the NaHCO 3 treatment in this paper, 1539 enriched OTUs of microbes were found, and the interaction networks were different. Moreover, NYJ negatively interacted with 1 species and positively interacted with 11 species in the NYJ + NaHCO 3 treatment of this study. Thus, application of NYJ affected microbial community structure and interaction network under NaHCO 3 stress. These results were consistent with the improved plant growth and the decreased contents of MDA and H 2 O 2 in the NYJ + NaHCO 3 treatment in comparison to the NaHCO 3 treatment in this paper. So NYJ application might mitigate NaHCO 3 stress by affecting microbial community in rhizospheric soil.
Symbiotic relationship between plant and microbes affects plant growth [49], and soil microbes modulate plant immune system [50]. Therefore, genes implicated in the symbiont and immune response were investigated in this study to understand the NaHCO 3 stress mitigation by NYJ application. In the NYJ + NaHCO 3 treatment in comparison to the NaHCO 3 treatment of this paper, the relative abundances of symbiont-or immune-related genes in GO:0002764, GO:0002757, GO:0002253, GO:0045087, GO:0045088, GO:0045089, GO:0002758, GO:0002218, GO:0044409, GO:0046718, and GO:0051828 were increased. Meantime, there were 2238 symbiont-and immunerelated genes, which were from major functional genera including Massilia and had the significantly increased relative abundances in this study. These results indicated that NYJ application might affect symbiotic relationship between plant and microbes and make some soil microbes to modulate plant immune system under NaHCO 3 stress.
Proteobacteria, Actinobacteria, and Bacteroidetes were the dominant phyla in rhizospheric soil of this paper whether using the shotgun metagenomic sequencing or the 16S rRNA amplicon sequencing.
Proteobacteria and Bacteroidetes are important contributors to N cycling in soil [51] and will induce phosphatase [52]. Actinobacteria plays roles in soil C mineralization [53]. Thus, N-related urease, phosphatase, and C-related sucrase were investigated in rhizospheric soil of this study to explore the NaHCO 3 stress mitigation by NYJ application. CAT in soil, an indicator of aerobic microbial activity [54], was also detected. Application of Pseudomonas putida CSY-P1 into soil increased soil enzyme activity under stress of p-hydroxybenzoic acid [17]. In this paper, inoculating NYJ into NaHCO 3supplemented and cucumber-planted soil elevated the activities of urease, phosphatase, sucrase, and CAT in rhizospheric soil, which were consistent with the improved plant growth and the decreased contents of MDA and H 2 O 2 in seedlings of the NYJ + NaHCO 3 treatment in comparison to the NaHCO 3 treatment. These results suggested that NYJ application might mitigate NaHCO 3 stress by inducing soil enzymes.
Plant growth-promoting traits of rhizospheric microbial communities contribute to root colonization of microbes and act beneficially on plant growth [55]. The strain NYJ in this study belonged to PGPR but did not produce chitinase. Compared the NYJ + NaHCO 3 treatment to the NaHCO 3 treatment, the relative abundances of chitinase-encoding genes were increased, and 16 chitinase-encoding genes of Massilia showed the enhanced relative abundances. These results indicated that NYJ application might induce microbes in rhizospheric soil to produce chitinase, which would play roles in NaHCO 3 stress mitigation by NYJ application. Similarly, the chitinase isolated from Pseudoalteromonas is salt-tolerant [56]. We also found some major functional genera, in which NYJ application induced genes implicated in nitrogen metabolism, IAA production, general plant growth promoting activity, phosphate solubilization, siderophore, polyamine, and extracellular polysaccharide under NaHCO 3 stress. These genera might contribute to plant growth [55] when NYJ being applied to mitigate NaHCO 3 stress in cucumber.

Conclusions
NaHCO 3 stress altered rhizospheric microbial community structure and interaction network, while 669 NaHCO 3 -tolerant species and 208 NaHCO 3 -tolerant genera were obtained from rhizospheric soil. NaHCO 3 -tolerant NYJ had plant growth-promoting capacity, and the strain negatively interacted with 1 species and positively interacted with 11 species under NaHCO 3 stress. As a result, application of NYJ into soil affected microbial community in rhizospheric soil, increased the activities of urease, phosphatase, sucrase, and CAT, changed the symbiotic relationship between plant and microbes, modulated plant immune system, and elevated the relative abundances of genes implicated in plant growth-promoting traits, thus mitigating NaHCO 3 stress in cucumber.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The sequence of a 16S rRNA gene fragment in the strain NYJ and the raw data of the 16S rRNA amplicon sequencing and the shotgun metagenomic sequencing for rhizospheric microbes have been 19 submitted to the NCBI database under accession number MN197602, PRJNA615957, and PRJNA616167, respectively.

Competing interests
The authors declare that they have no competing interests.  Relative abundance of microbial phyla in rhizospheric soil based on the shotgun metagenomic sequencing (a) and 16S rRNA amplicon sequencing (b). Control, supplemented with water; NaHCO3, supplemented with NaHCO3; NYJ, supplemented with NYJ; NYJ+NaHCO3, supplemented with NYJ and NaHCO3.
28 Figure 3 Heatmap of the relative abundance of microbial genera in rhizospheric soil by using the shotgun metagenomic sequencing. Control, supplemented with water; NaHCO3, supplemented with NaHCO3; NYJ, supplemented with NYJ; NYJ+NaHCO3, supplemented with NYJ and NaHCO3. Heatmap of the relative abundance of genes, which are implicated in osmotic stress and annotated by mapping the metagenomic sequencing reads against the GO database.
Control, supplemented with water; NaHCO3, supplemented with NaHCO3; NYJ, supplemented with NYJ; NYJ+NaHCO3, supplemented with NYJ and NaHCO3. 32 Figure 7 Heatmap of the relative abundance of genes, which are implicated in the adaptation to sodium ion and annotated by mapping the metagenomic sequencing reads against the GO database. Control, supplemented with water; NaHCO3, supplemented with NaHCO3; NYJ, supplemented with NYJ; NYJ+NaHCO3, supplemented with NYJ and NaHCO3. Heatmap of the relative abundance of genes, which are implicated in plant growthpromoting traits and annotated by mapping the metagenomic sequencing reads against the GO database. Control, supplemented with water; NaHCO3, supplemented with NaHCO3; NYJ, supplemented with NYJ; NYJ+NaHCO3, supplemented with NYJ and NaHCO3.


extfieldsofstudy: Medicine Biology

s2fieldsofstudy:
Agricultural And Food Sciences
Biology

provenance: AgriculturalAndFoodSciences,2019-2020-Agricultural And Food Sciences-2020.jsonl.gz.gz

text:
A Partially Phase-Separated Genome Sequence Assembly of the Vitis Rootstock ‘Börner’ (Vitis riparia × Vitis cinerea) and Its Exploitation for Marker Development and Targeted Mapping

Grapevine breeding has become highly relevant due to upcoming challenges like climate change, a decrease in the number of available fungicides, increasing public concern about plant protection, and the demand for a sustainable production. Downy mildew caused by Plasmopara viticola is one of the most devastating diseases worldwide of cultivated Vitis vinifera. In modern breeding programs, therefore, genetic marker technologies and genomic data are used to develop new cultivars with defined and stacked resistance loci. Potential sources of resistance are wild species of American or Asian origin. The interspecific hybrid of Vitis riparia Gm 183 x Vitis cinerea Arnold, available as the rootstock cultivar ‘Börner,’ carries several relevant resistance loci. We applied next-generation sequencing to enable the reliable identification of simple sequence repeats (SSR), and we also generated a draft genome sequence assembly of ‘Börner’ to access genome-wide sequence variations in a comprehensive and highly reliable way. These data were used to cover the ‘Börner’ genome with genetic marker positions. A subset of these marker positions was used for targeted mapping of the P. viticola resistance locus, Rpv14, to validate the marker position list. Based on the reference genome sequence PN40024, the position of this resistance locus can be narrowed down to less than 0.5 Mbp on chromosome 5.

INTRODUCTION
Downy mildew of grapevine is caused by the oomycete Plasmopara viticola. It is a serious threat for viticulture especially in humid and warm environments. Intensive chemical plant protection, frequently through the use of copper, is necessary to protect grapevines under disease-promoting conditions. Resistance breeding allows for the selection of new cultivars that require a reduced number of chemical treatments. Genetic analyses of resistance loci have promoted such breeding initiatives that are today oriented towards stacking (combining) of several resistance loci (Töpfer and Eibach, 2016). In the recent past, therefore, a number of resistance loci for downy mildew have been identified, and these have been localized within grapevine genomes. Around two dozen loci designated Rpv (Resistance Plasmopara viticola) are listed in the Vitis International Variety Catalogue [VIVC, www.vivc.de;and (Hausmann et al., 2019)]. Some of them are of elite genetic background, such as Rpv1 (Merdinoglu et al., 2003), Rpv3 (Welter et al., 2007;Bellin et al., 2009), Rpv10 (Schwander et al., 2012), and Rpv12 (Venuti et al., 2013), and these can be used in breeding programs for selection of new cultivars. Other loci are less well characterized and are essentially known only from wild species (Marguerit et al., 2009;Blasi et al., 2011;Ochssner et al., 2016). To allow their application, these loci need to be introgressed into a favorable Vitis vinifera genetic background. Such an introgression can be accelerated considerably by using marker-assisted selection or marker-assisted back crossing (Herzog et al., 2013). To start such an approach, a locus should be well characterized, and sufficient closely linked markers should be available. Within plant breeding programs, simple sequence repeats (SSRs) and single nucleotide variant (SNVs) markers are widely applied. The discovery and development of this marker types by using NGS data has been successfully demonstrated for various plant species (Taheri et al., 2018).
One of the genetically useful genotypes as a source for resistance is the rootstock cultivar 'Börner'. It is an interspecific hybrid of V. riparia Gm 183 x V. cinerea Arnold; 'Börner' is thus derived from individuals of two Vitis species native to North America that carry several resistances. Among these are resistances to phylloxera [Rdv1 (Zhang et al., 2009)], black rot [Rgb1 and Rgb2 (Rex et al., 2014)], and to downy mildew [Rpv14 (Ochssner et al., 2016)]. To support the molecular analysis of these and other loci, we generated a first draft genome sequence assembly of the highly heterozygous 'Börner' genome. We produced comprehensive lists of SSR and SNV marker positions based on NGS data from 'Börner,' and we validated these marker positions randomly as well as by using Rpv14 on chromosome 5 (chr.5) as an example target.
The first grapevine reference genome sequence was generated based on Sanger sequencing technology using the 'Pinot noir' inbreeding line PN40024, which is estimated to be homozygous for 93% of the genome, thus enabling a successful assembly [(Jaillon et al., 2007), current version 12x.v2]. However, PN40024 suffers from inbreeding depression, and selfing grapevine for homozygosity is a difficult and a time-consuming process. Up until 2018, other available grapevine genome sequences were derived from heterozygous genotypes (Velasco et al., 2007;Adam-Blondon et al., 2011;Di Genova et al., 2014), and this posed a challenging sequencing and assembly problem that is not yet fully solved. In 2007, it was impossible to achieve a high-quality genome assembly for 'Pinot noir' due to its heterozygosity since merged allelic sequences made it impossible to separate both haplophases. As NGS techniques developed progressively, the number of sequenced plant genomes has rapidly increased, breaching the 100 mark some time ago (Michael and Jackson, 2013;VanBuren et al., 2015). Only a few genome sequences are based purely on the Sanger sequencing technique (e.g., Arabidopsis and grapevine). A larger number was generated based on mixed data types, including Sanger,454,and Illumina (e.g.,potato,sugar beet,and Brassica rapa) or only Illumina-derived sequences (e.g., eggplant and Norway spruce). The quality of the latter, generated by early NGS technologies, does not reach the quality of the older BAC-based Sanger assemblies despite higher read data coverage (Alkan et al., 2011;Burton et al., 2013). Although long reads became the best option for generating high-quality genome sequences (Li et al., 2018), the high throughput and low costs of Illumina sequencing has made it economical to sequence plenty of genomes at least to a draft state or to re-sequence individual genotypes. A major problem and current challenge in plant genome sequence assembly is resolving the repeat structures and proper continuous separation of haplophases within heterozygous or polyploid plants (VanBuren et al., 2015). Repeat sequences are very common in plant genomes, including repeats that are more than 10 kbp in length. Therefore, resolving plant genome complexity with short NGS reads, mainly 100-300 nt in length, is a big challenge. Such assemblies result in collapsed sequences and breakdown at the borders of repetitive regions, and this causes shortening of contigs. Some improvements were reached by using a variation of paired-end libraries (VanBuren et al., 2015) combined with mate-pair libraries. Long-insert libraries were of particular value to bridge such gaps up to about 12 kbp.
Here, we utilized sequencing data mainly generated by 454 and Illumina technologies to generate a draft assembly of the 'Börner' genome sequence with the main goal of allowing easy access to markers. The approach was complemented by sequencing BAC ends from a library created from 'Börner' DNA. We hypothesized that the interspecific nature of 'Börner' and the genetic distance of its two wild ancestor species, V. riparia and V. cinerea, could be helpful to dissect and assign the haplophases of 'Börner' in its parental genomes even with Illumina reads. Having at least to a certain extent phased contigs at hand, several downstream applications become feasible. The alignment of 'Börner' contigs to the reference PN40024 sequence revealed many polymorphic sites, including those at SSR positions, as being sources of potential high value for, to give an example, the purpose of grapevine resistance breeding. The detection of single SNVs or groups of nearby SNVs within a multiple alignment of two 'Börner' contigs and the PN40024 reference sequence are both considered to be useful in a bulked segregant analysis (BSA). The identification of markers linked to loci or genes causing disease resistance by applying BSA has been shown in plants (Michelmore et al., 1991). Pooling equal amounts of DNA is a practical way to reduce the cost of large-scale association studies to identify susceptibility loci for diseases (Sham et al., 2002;Zou et al., 2016). The pooling allows us to measure allele frequencies in groups of individuals by using far fewer genotyping assays than required when genotyping all individuals. Here, we demonstrate the power of heterozygous positions within one species to physical narrow down a candidate locus for downy mildew resistance (Rpv14).

BAC Library
To construct a BAC library of the grapevine cv. 'Börner' young frozen leaf material was sent to Bio S&T Inc. (Montreal, Canada) as a service provider for the generation of BAC libraries. The DNA was inserted into the vector pIndigoBAC-5 (Epicentre, Madison, USA; ENA/GenBank accession number EU140754) cut with HindIII. Competent cells of Escherichia coli strain DH10B were transformed with the ligation mix and plated on LB medium supplemented with 12.5 µg/mL chloramphenicol (CM), 40 µg/mL X-GAL and 0.2 mM IPTG. White colonies were picked and transferred in 384-well microtiter plates filled with LB freezing medium (Zimmer and Verrinder Gibbins, 1997) containing CM. Altogether, 159 microtiter plates with 384 wells were generated, and these were duplicated and stored at -80°C. The insert size was estimated by pulsed-field gel electrophoresis to be on average 93 kbp. The library coverage is almost 12 haploid genome equivalents based on a genome size of about 500 Mbp (Lodhi and Reisch, 1995).

Processing of BAC End Sequences (BES)
Trace files were processed by the base calling program phred (www.phrap.org, version 0.071220.c). Vector sequences were trimmed using the program blastall (Altschul et al., 1990) in version 2.2.24 with blast parameters set to the ones suggested by NCBI (www.ncbi.nlm.nih.gov/tools/vecscreen/about/). Lowquality sequences were masked by applying a sliding window approach (size: 10 bp, minimal average score > 13), and the longest unmasked subsequence was taken as a high-quality sequence if the length was ≥ 80 bp. Slippage reads were removed as described before (Dohm et al., 2012). The dataset was filtered for contamination with E. coli sequences (blastall matches with an e-value ≤ 1e-50). Afterwards, the sequences were singularized per BAC end. A total of 69,444 high-quality 'Börner' BES were submitted to ENA/GenBank and can be accessed under accession numbers KG622866 -KG692309.

NGS Sequencing Library Preparation
'Börner' DNA, extracted from young leaf material with a commercial kit (Qiagen) according to the instructions of the manufacturer, was sequenced using 454 GS-FLX Titanium (Roche), IonTorrent PGM (life technologies), Illumina GAIIx, MiSeq, and HiSeq1500 instrumentation. The preparation of Illumina paired-end (PE) and single-end (SE) sequencing libraries was performed according to the Illumina TruSeq DNA Sample Preparation v2 Guide. DNA was fragmented by nebulization. After end repair and A-tailing, individual indexed adaptors were ligated to the DNA fragments, thus allowing for multiplexed sequencing. The adaptor-ligated fragments were size selected on a 2% low melt agarose gel to a size of 300-600 bp. After enrichment PCR of fragments that carry adaptors on both ends, the final libraries were quantified with PicoGreen. The average fragment size of each library was determined on a BioAnalyzer High-Sensitivity DNA chip. Libraries were sequenced either 100 nt on an SE run using one lane of a GAIIx run (ENA/SRA accession numbers ERR2729619 and ERR2729620) or sequenced 2 x 100 nt PE on GAIIx (ERR2729739, ERR2729740, ERR2729741, ERR2729743, and ERR2729744) or HiSeq1500 (ERR2729742) sequencers. A mate-pair library was constructed according to the Illumina Mate Pair Sample Preparation v2 Guide. DNA was fragmented to a size of 3-4 kbp by using Hydroshear. Fragment ends were repaired and biotinylated. Afterwards, fragments of 3.5 kbp size were selected by gel electrophoresis. The resulting molecules were circularized. Removal of linear fragments was achieved by DNA exonuclease treatment. The circular molecules were randomly fragmented by nebulization. The biotinylated fragments were purified using streptavidin-coated magnetic beads. Several steps were carried out as described for the PE library construction. Finally, the library was sequenced in pairedend modus using an Illumina MiSeq generating 2 x 36 nt PE reads (ERR2729745).
GS-FLX sequencing libraries (454) were prepared from 10 micrograms of DNA according to the GS-FLX-Titanium Rapid Library Preparation Method Manual. DNA was fragmented by nebulization. Purified DNA fragments were subjected to end repair and subsequently to adaptor ligation. Large DNA fragments were selected via AMPure Beads. After emPCR and enrichment of DNA carrying beads, sequencing was performed on several runs on a Roche/FLX sequencing instrument using Titanium XLR70 sequencing kits (ERR2728011-ERR2718029). Mate-pair libraries (ERR2728288-ERR2728292) with varying insert length (multi-span PE reads, 3 kbp, and 8 kbp) were prepared from 15-20 micrograms of genomic DNA according to GS-FLX-Titanium manual method and using the GS FLX Titanium Paired End Adaptors kit and sequenced on a Roche/ FLX sequencer.
The library construction for IonTorrent PGM (life technologies) sequencing was performed according to the Ion Xpress™ Plus gDNA and Amplicon Library Preparation guide for a 200 bp library. DNA was fragmented by nebulization. The ends were repaired and ligated to adapters. Size selection of DNA molecules was performed via E-Gel ® SizeSelect™ Agarose Gel to an average insert size of 350 bp. After enrichment of the template carrying Ion Sphere particles on an Ion OneTouch, sequencing was performed using a 318 chip on an Ion Personal Genome Machine™ system (ERR2729812).
In addition to 'Börner' we used the V. vinifera cultivar 'Pinot noir' for DNA extraction, library construction, and Illumina HiSeq1500 sequencing (2 x 100 nt PE) according to the workflow described above. The read data were submitted to the ENA study PRJEB36326 with the read accessions ERR3834961 and ERR3834962).

NGS Read Pre-Processing
Illumina reads were quality trimmed and adapter clipped using Trimmomatic (Bolger et al., 2014) with all available adapter sequences from Illumina. Reads from the 454-and PGM/Ion Torrent-sequencing platform were trimmed with the CLC Genomics Workbench using the trim points of the SFF-files. Reads shorter than 36 nt were discarded. After pre-processing, 75.62 Gbp raw data (=90.9% of the raw bases) were available for de novo assembly.

Assembly, Post-Processing, and Validation
A de novo hybrid assembly of the 'Börner' genome was generated with the CLC Genomics Workbench (v. 8.0) using the certain parameters: length fraction 0.5, similarity fraction 0.8, word size 20, and bubble size 50. For the re-mapping of the reads to the contigs that the CLC assembler carries out, the filtering thresholds were optimized to a length fraction of 0.8 (default: 0.5) and a similarity fraction of 0.95 (default: 0.9). From 75.62 Gbp of read data, representing nearly 80-fold coverage for each of the two 'Börner' haplophases, 338,484 contigs with a total length of 652.17 Mbp were calculated. A total of 91.3% of the reads could be mapped back to the assembly, indicating a high probability that the sequence data was integrated correctly into the contigs. This raw assembly was designated BoeWGS0.8.1 (Supplementary Table 2).
The assembly metadata include a value for the number of reads that cover a given position and averages of these values for each contig. The value, which is essentially the total length of aligned reads divided by the assembly size, is referred to below as "alignment depth" (AD) and was used to deduce if a given contig represents a single haplophase or contains merged sequence information from both haplophases.
To test the completeness of the assembly we mapped the 69,444 'Börner' BES on BoeWGS0.8.1. Mapping was done with the program BLAT (Kent, 2002). The majority (43,008 BES = 61.9%) mapped with 99% identity or better and for a length fraction of 0.5. About 4.5% BES matched the end of contigs with between 20 to 50% of their length, but their mapping was not unique to one contig within the assembly. A small proportion of these (733 BES) produced multiple matches, suggesting that these contain mainly repetitive sequences. To avoid connecting contigs from different phases, we refrained from scaffolding since we had no information about the phase from which the BACs/ BES were derived.
The raw assembly BoeWGS0.8.1 was checked for contaminations by searching against a custom database with sequences of human, yeast, mold, several fungal and viral pathogens, and insects, such as aphids and mites, with blast+ (e-value <=1e-50) (Camacho et al., 2009). A BLAST (Altschul et al., 1990) run querying GenBank's non-redundant protein database (nr) completed the filtering process by using the same cut off criteria. Contigs matching to contaminants were also characterized by very low AD (1-20x). The removal of PhiX as well as sequences resembling cloning vectors caused a reduction of 213,834 bp and a total removal of three contigs. Plastid sequences from V. vinifera were detected in 67 contigs and removed as well. In addition, contigs smaller than 500 bp were discarded. Although nearly 30% of all initial contigs were rejected by this filter, the total assembly size was only reduced by about 5%. A total of 238,193 contigs remained after the different filters and were included in the final assembly, which was designated BoeWGS1.0 (Biosample PRJEB28084/http:// www.ebi.ac.uk/ena/data/view/PRJEB28084). The results were evaluated for standard assembly statistics parameters (Supplementary Table 2).
The BoeWGS1.0 contigs were aligned to the reference genome sequence PN40024 12x.v2 (https://urgi.versailles.inra. fr/download/vitis/12Xv2_grapevine_genome_assembly.fa.gz) with blast+ with the dust filter turned off and using an e-value cutoff of 1e-50 and a culling limit of 1. As a result of this analysis, each mapping contig from BoeWGS1.0 receives location information derived from PN40024. If exactly two BoeWGS1.0 contigs were assigned to one locus of the reference sequence, we took this as indication that the assembly separated the two 'Börner' haplophases. We have referred to these contigs, which may contain allelic sequence information, as "contig pairs".
BWA-MEM (Li, 2013) was applied to map all PE Illumina reads to the final assembly sequence using the -M option to mark short splits as secondary. The resulting mapping was subjected to inspection by REAPR (Hunt et al., 2013) to assess the assembly quality based on coverage, distance of PE reads, and orientation of PE reads.
BUSCO v3 (Simão et al., 2015) was deployed to analyze the recovery of the 'embryophyta odb9' benchmarking sequences in the assembly. This analysis was performed in the 'genome' mode using an e-value cutoff of 1e-3 and considering up to three candidate regions per BUSCO hit.

Calculation of Heterozygosity
Jellyfish v2.2.4 (Marcais and Kingsford, 2011) was applied on the reads of ERR2729742 to calculate histograms for k-mer lengths 19, 21, 23, and 25. The results were individually subjected to GenomeScope (Vurture et al., 2017) using default parameters for calculation of the heterozygosity. The average of all four resulting predictions was calculated to avoid k-mer length dependency of the final result. For comparison, we included Illumina read data from 'Pinot noir' (V. vinifera). As for the analysis of 'Börner,' about 40 Gbp Illumina PE data were used and the analyses were performed with the same parameters.

Genome-Wide Variant Detection
A strict mapping of 245,325,650 'Börner' PE reads (2x100 nt) on the PN40024 12x.v2 reference sequence with an average coverage of 50x formed the basis for the variant detection. The mapping was done with a CLC Genomics Workbench toolkit, the similarity fraction set to 0.95 and the length fraction to 0.9. Non-specific matches were discarded. Variants [single nucleotide variants (SNVs), insertions and deletions (InDels), as well as multi nucleotide variations (MNVs)] were detected with CLC's Probabilistic Variant Caller using default parameters. Only biallelic SNVs with a read coverage between 10 and 90 were kept for downstream analyses, whereas tri-allelic and poly-allelic variants, InDels, and MNVs were ignored.

Simple Sequence Repeat Detection
SSR detection was done with MISA [MIcroSAtellite identification tool; http://pgrc.ipk-gatersleben.de/misa/misa. html, (Thiel et al., 2003)] using default parameters. Possible SSR marker positions were evaluated by comparison of the candidate SSRs plus adjacent +/-200 bases from BoeWGS1.0 with the SSRs of the PN40024 12x.v2 reference sequence. The comparison was performed with blastn with dust filter turned off and 1e-50 as e-value cutoff. SSR candidates where a BoeWGS1.0 contig pair matched one coordinate in the reference sequence and differences did not exceed 40 bp (identity >=90%) were classified as useful for SSR marker development. The SSR candidates were named "SSR [chr#]_[contig-ID]" with a trailing number if more than one SSR candidate was detected for one contig pair.

SSR Marker Validation
For the validation of heterozygous SSRs within the BoeWGS1.0 assembly and variants in relation to the reference sequence, 53 SSR candidate positions were randomly selected for SSR marker validation. The SSRs from different genomic regions were selected with a minimum of two nucleotides and of six contiguous repeat units in the reference sequence. Primer3 (Untergasser et al., 2012) was used to design the amplimers (two primers directing towards each other based on the source sequence) fitting to the BoeWGS1.0 contig sequences as well as to the reference sequence PN40024 (Supplementary Table 1). The expected amplicon size of the markers was set from 150 to 400 bp and the primer size from 18 to 27 nt. DNA from 'Börner' and V3125 (see below) was used to check the polymorphism of the SSR marker sequences by PCR and gel electrophoresis on a 3% agarose gel. A few PCR products with expected smaller sequence variations were applied to Sanger sequencing and analyzed at the bp level by multiple sequence alignments. Verified heterozygous markers either in 'Börner' and/or in relation to V3125 were analyzed by fragment analysis using DNA of the crossing parents (V3125 and 'Börner') as well as the parents of 'Börner' as described before (Ochssner et al., 2016)

Locus-Guided Variant Calling and Marker Generation
The BoeWGS1.0 contig mapping on the PN40024 12x.v2 reference sequence was evaluated for a genomic region on chr.5. The BLAST results from the post-processing analysis described above were used. Regions where a multiple alignment composed of exactly two BoeWGS1.0 contigs (i.e., contig pairs) and parts of chr.5 from the PN40024 12x.v2 reference sequence was built were considered as genome regions representing both haplophases if the contig pairs selected showed the expected AD of about 75x (75 ± 25).
Each contig pair was aligned with edialign (from the EMBOSS suite, version 6.2) using default parameters. Alignments without N-stretches and with lengths of at least 500 bp were subjected to SNV detection. The list of contig pairs used for SNV detection on chr.5 is available as Supplementary Table 3. The downstream analysis focused on grouped SNVs between BoeWGS1.0 contig pairs that are bi-allelic after comparison with the PN40024 reference sequence.

Amplimer Design for the Rpv14 Locus
BoeWGS1.0 contig pairs on chr.5 were used for SNV-based marker development addressing the Rpv14 locus. Amplimers were designed using Primer3 (Untergasser et al., 2012) preferably for contig pair alignments with a rate between 0.4 and 1 SNVs per 100 bp (Supplementary Table 1 and Supplementary Table  3). We targeted an expected amplicon size between 400 and 530 bp. Deduced amplimers were tested to be unique in BoeWGS1.0 and PN40024 12x.v2.

Plant Material, Bulk Set Up, DNA Extraction, and PCR
A biparental F1 mapping population comprising 202 individuals of the grapevine breeding line V3125 ('Schiava Grossa' x 'Riesling') and 'Börner' (V. riparia Geisenheim 183 x V. cinerea Arnold), cultivated in the field at the JKI Institute for Grapevine Breeding Geilweilerhof, Siebeldingen, Germany, has been phenotyped three times for downy mildew leaf resistance as described previously (Ochssner et al., 2016). For SNV genotyping by amplicon sequencing, genomic DNA was extracted from the parents V3125, and 'Börner' and selected F1 genotypes of the mapping population. Building bulks of DNA from all the individual samples and typing one SNV marker at a time can save valuable template DNA and has been successfully utilized in microsatellite markers (Barcellos et al., 1997) and SNPs (Shaw et al., 1998). Two bulks with 10 individuals each were generated from the mapping population. The first (R) and second (S) bulks encompass resistant and susceptible genotypes, respectively.
Genomic DNA was extracted from young leaf material with a commercial kits (Qiagen) according to the instructions of the manufacturer. DNA samples were quantified using the Qubit™ dsDNA BR Assay Kit with the Qubit ® 2.0 Fluorometer (Invitrogen, Life Technologies, Darmstadt, Germany).
Genomic DNA from the parents of the mapping population V3125 x 'Börner', as well as from the two bulks, was used as template in 25 µl PCR reactions with the two marker-specific primers and 1 ng DNA per individual under standard conditions. The amplicons obtained were purified and sequenced from both directions using the ABI Prism BigDye Terminator chemistry on an ABI Prism 3730 sequencer (Applied Biosystems).

Determination of SNV/Allele Ratios
For the estimation of allele ratios, only bi-allelic BoeWGS1.0 SNVs that differentiate the two 'Börner' phases were considered, i.e., SNVs between BoeWGS1.0 contig pairs. Since only bi-allelic SNVs were considered, these positions are often homozygous in the V3125 genotype for one of the two 'Börner' alleles. Exceptions to these criteria (V3125 displays identical SNV heterozygosity to 'Börner') were excluded from the analysis. As a result, the allele frequency contributed from the 'Börner' parent has to be 50% and is 100% for the V3125 parent when counting the variant that is present in V3125. Following these assumptions, the expected frequency for each single allele in a bulk of the unselected F1 individuals of this mapping population is 75%. We have estimated the SNV frequencies within the parental genotypes in order to confirm the determination and to gain hints for problematic amplimers that amplify the haplophases with a bias for one of the two.
The estimation of the SNV frequencies among the pooled DNA samples (Pool R = resistant; Pool S = susceptible) and comparison with the parental lines was carried out by using the tool QVSanalyzer (Carr et al., 2009). The Sanger sequence trace files for determination of the relative proportions of two sequence variants were analyzed per batch. The generated output files contain details of the examined sequence variant ratios for individual samples as well as summary statistics. The area below the peak at the position of the targeted SNV was calculated and set into context with the surrounding peaks for each sample and the corresponding trace file. For graphical presentation of the results, the ratios were converted into percent values for SNV or allele frequency. SNVs considered in the analysis are listed in Supplementary Table 3).

A Draft Assembly of the 'Börner' Genome Sequence
The data for the 'Börner' genome sequence were obtained by whole genome shotgun (WGS) sequencing. The assembly, designated BoeWGS1.0, has a total size of 618.3 Mbp, and the N50 sequence size was 4,255 bp (Supplementary Table 2). We evaluated the average AD for the contigs in the assembly and classified the contigs with respect to the expected depth of aligned reads for paired contigs (separated haplophases) and contigs with merged sequences (both haplophases combined into one contig, see Supplementary Figure 1). A significant proportion of the contigs showed an AD between 50 and 100 (75 ± 25), which can be considered as the expected value for a haplophase-specific contig. The range has been deduced from the total amount of reads (75 Gbp), 2 x 500 Mbp for the expected fully diploid sequence length, and an interval of ±33%.
REAPR flagged a total of 14,007 regions (5.9%) as erroneous. Fragment coverage distribution errors within a contig accounted for 4,250 cases (1.8%), while the same error type over a gap added 5,674 cases (2.4%). The remaining cases were contributed by low coverage within a contig or across a gap with 3,723 (1.6%) and 360 (0.2%) cases, respectively. In summary, the assessment by REAPR indicated a short-read assembly from heterozygous material of acceptable to good quality. The heterozygosity of 'Börner' was calculated to be about 3.1% independent of the kmer length used in the analysis. This is a quite high value in comparison to the V. vinifera variety 'Pinot noir' for which we calculated a heterozygosity of 1.5% regardless of k-mer length (Supplementary Table 4).
BUSCO, i.e., the Benchmarking of Universal Single-Copy Orthologs, revealed the presence of 55.9% of all 1,440 benchmarking genes from the reference set for embryophyta, 42% of the total 1,440 as single copy and 13.3% as duplicated. In addition, fragments of 19.1% benchmarking genes were detected. Only 25% of the benchmarking genes were not detected in the assembly. It should be noted that the duplicated BUSCO genes can be explained by detection of two allelic versions in BoeWGS1.0 contig pairs. Taken together, the BoeWGS1.0 assembly represents a typical short-read assembly of a heterozygous genotype. Since we wanted to address marker development, we optimized the assembly for phase separation and not for continuity.

BoeWGS1.0 Contigs Representing Different Haplophases
The BoeWGS1.0 contigs were mapped against the V. vinifera reference genome sequence PN40024 12x.v2. About 210,444 (88.3%) could be mapped successfully with at least 30% of their length (https://doi.org/10.4119/unibi/2938185). The ratio of the total length of the mapped BoeWGS1.0 contigs to the length of the reference chromosomes is 1.14 on average (1.02 to 1.31; Supplementary Table 5), indicating uniform assembly quality and an overall homogeneous synteny between 'Börner' and PN40024. The median of the average contig ADs for all 19 chromosome mappings was 86.3 (Supplementary Table 5). Conspicuously, it is only the genetically unassigned chromosome "Ukn," or "unknown," where the mapped BoeWGS1.0 contigs have an AD of 248 on average. The remaining, unmapped contigs contain sequences that are too diverse from the PN40024 reference sequence to be mapped.
More than one quarter (29.2%; 142 Mbp) of the PN40024 reference sequence was covered by BoeWGS1.0 contig pairs, again suggesting that the two 'Börner' haplophases were partially separated in the assembly. Most of these contigs show on average an AD in the expected range. Another quarter of the reference (28.7%; 139 Mbp) was covered by one BoeWGS1.0 contig, and these contigs often (but not always) display higher AD values. In these cases, either the two 'Börner' haplophases were merged into one contig during assembly, or only one of the two haplophases displayed sufficient similarity to become mapped.
Due to repetitive sequences, a small fraction of the PN40024 reference sequence (3.7%; 18 Mbp) was matched by more than two contigs. Nearly 38.5% (187 Mbp) of the reference sequence (including 3% N-stretches) was not covered by any stringently matching BoeWGS1.0 contig (Figure 1).

Homozygous and Heterozygous Variant Frequencies
Genomic variants including SNVs contribute not only to intraspecific diversity but they are also candidates for valuable molecular markers. By mapping 'Börner' NGS reads against the reference sequence and subsequent variant calling, almost 5 million highly reliable bi-allelic SNVs were detected (https:// doi.org/10.4119/unibi/2938178 and https://doi.org/10.4119/ unibi/2938180).
Half of the 4,996,490 SNVs are heterozygous SNVs (2,536,406), which means one of the BoeWGS1.0 contigs shows the same nucleotide as the reference at a given position. The other half (2,460,084) are homozygous SNVs, describing variants only existing between PN40024 and BoeWGS1.0. The frequency is on average 1 per 77 bp (Supplementary Table 6), ranging from 1/70 to 1/82 bp for chr.3, chr.8, and chr.16. In addition to SNVs, MNVs and small InDels were called, dropping the overall variant frequency slightly to 1 variant per 68 bp.

SSR Detection for Marker Development
Simple sequence repeat (SSR) markers are often used to study molecular diversity or heterozygosity as well as for genetic mapping in grapevine. By comparing SSR positions in the grapevine reference genome sequence with BoeWGS1.0, a total of 10,820 putative SSR marker positions ("candidates") with different unit sizes (two to six) and repeat numbers (up to 27) were deduced ( Table 1). Regions of the reference matched by a BoeWGS1.0 contig pair were exploited. Out of the 10,820 positions, 12% (1,313) were monomorphic and 38% (4,110) were tri-allelic in an alignment of the three sequences. The remaining 50% were either bi-allelic or showed more than one motif in the SSR region. The more than 4,000 tri-allelic SSR candidates are the most valuable ones with regard to genetic mapping in a broad range of accessions (Supplementary Table  7). For various reasons, not all of them are new; for example, candidate SSR chr1_1203 corresponds to the established marker GF01-03, SSR chr1_1083 to marker GF01-21, and SSR chr1_1248 to marker GF01-22 . However, the identification of already existing markers convincingly indicates that our candidate list is valid.
A set of 45 tri-allelic and eight bi-allelic SSR candidates with a unit size of two or three randomly selected from all  chromosomes was used to validate the power and reliability of SSR candidate prediction. The bi-allelic candidates were included to estimate the chance of false negative predictions. In this subset of eight candidates, we tested six candidates; the reference matched only one contig of BoeWGS1.0, and two candidates were not found on a contig pair but shared a unit size and number. Out of the total set of 53 SSR sequences, 38 could be verified by either PCR and subsequent agarose gelelectrophoresis, Sanger sequencing, or fragment analysis (Supplementary Table 8). The 15 negative cases were caused by unspecific primers that amplified various fragments instead of allelic amplicons. We did not optimize the amplimers empirically in a second round of amplimer design. As expected, 'Börner' amplicons with predicted small sequence variations between the alleles were analyzed more successful by Sanger sequencing than those with larger sequence variations. In six out of eight cases, the predicted bi-allelic SSR sequences were verified as monomorphic within 'Börner'. Analysis of the sequenced nucleotides before and after the SSR region confirmed the assumption that these genome regions were less heterozygous than others (or monomorphic), and they were therefore only represented by a single contig in the BoeWGS1.0 assembly. For the SSR chr1_43023 candidate, the new marker GF01-59 was established, although for this candidate only one BoeWGS1.0 contig matched to the reference. Both Sanger sequencing and fragment analysis confirmed the existence of two 'Börner' alleles for GF01-59.
Sanger sequencing of amplicons from SSR candidate loci with larger differences (e.g., in unit size, unit number, etc.) was less successful. In addition, we observed in these cases higher heterozygosity between the alleles even outside the SSR itself. However, validation by agarose gel-electrophoresis worked out smoothly for these SSR candidates.
Since gel-electrophoresis did not allow us to verify exact allele sizes, these were determined by fragment analysis in 14 cases (Supplementary Table 8). Fragment analysis was performed with DNA of the exact parental lines of 'Börner' and V3125 in addition to 'Börner' itself. For 12 out of 14 SSR candidates different fragment sizes for the two 'Börner' alleles were detected that fit well to the predictions for both 'Börner' alleles. The validated SSR candidates-for which we generated complete documentation-received marker designations. Seven of these SSR markers turned out to be highly or fully informative because they discriminate all four alleles of the biparental cross V3125 x 'Börner'. One marker (GF01-60, derived from candidate SSR chr1_100137_2) was scored as monomorphic in this analysis, likely due to the small sequence difference of 2 bp in fragment size. In fact, exactly the 2 bp difference in fragment length were properly detected by Sanger sequencing, these but are surely below the reliable resolution for a marker assay based on fragment analysis. The last of the 14 SSR markers, GF06-19 (SSR chr6_800_1), did not show any scoreable fragment for 'Börner,' but amplified a single scoreable amplicon in V3125. Taken together, every method seems to have its limitations, but the different validation approaches demonstrated the high reliability and usability of our prediction approach for new SSR markers. Further investigation of the highly informative markers in the grapevine varieties 'Schiava Grossa', 'Riesling', 'Villard Blanc', 'Calardis Musque' and 'Pinot noir', as well as in the parents of 'Börner' showed a transferability of about 70 to 80%. Most of the assayed SSR positions were heterozygous within and between the different cultivars.

Targeted Mapping of the Downy Mildew Resistance Rpv14 by Bulk Segregant Analysis
To assist the introgression of the downy mildew locus Rpv14, which has been mapped in 'Börner' to the lower arm of chr.5 (Ochssner et al., 2016), into V. vinifera, closely linked markers are required. For SNV marker-based association studies for Rpv14, 25 amplimers were designed targeting the lower arm of chr.5 of 'Börner', and 17 of those turned out to be functional. As controls, two additional amplimers were designed, one at the top of chr.5 opposed to Rpv14 and one on chr.14. All primers were designed based on BoeWGS1.0 contigs and SNV data (see Methods). The 76 amplicons obtained for 19 markers and 4 template DNAs (susceptible V3125, resistant donor 'Börner', and the F1 bulks S and R) were sequenced, and the trace files were subjected to peak area determination and evaluation.
All predicted SNVs within or between the BoeWGS1.0 contigs could be verified. In one case, we observed a tri-allelic SNV and removed the respective data points from further analyses. For SNVs from physically distinct contigs, the peak ratios or SNV frequencies varied as initially expected for a bulked segregant analysis. By ordering BoeWGS1.0 contig pairs along the chr.5 of the reference genome sequence and allocating the corresponding SNVs or allele frequencies to these contigs, the changes in the frequencies of the pools are obviously correlating with the candidate locus ( Figure 2). For the control regions on chr.14 and at the top of chr.5, we hypothesized no selective pressure and therefore expected an allele frequency for the bulks of 75%. The observed frequencies on chr.14 fulfilled the expectation. The frequencies in the north of chr.5 in the R bulk (55-60%) differ to some extent from the expectation where, as in the S bulk, the SNV frequency is as expected. This may be due to the low number of genotypes within the bulks.
Even if only one out of the ten genotypes has no recombination between the target locus and the terminal part of the chromosome, the estimated frequency value changes theoretically by almost 10% of the expected value. Within the region of the candidate locus between approximately 19 and 21 Mbp, the resistant bulk has a nearly uniform estimated SNV frequency of 50%. This indicates that all resistant individuals in the bulk carry the same allele of the resistant 'Börner' parent (inherited from V. cinerea) (Ochssner et al., 2016). Since there are only a limited number of recombination events possible, the large interval is not surprising. At the same time, the SNV frequency in the susceptible bulk converges (continuously) to about 95%. Our results indicate that the variants of 'Börner', located between 19.67 Mbp and 20.6 Mbp in the corresponding reference sequence PN40024 are highly and specifically linked to the resistant trait. The center of the candidate locus with the highest linkage, showing SNV frequencies over 85% for the S bulk is an interval of 330 kbp between 20.31 and 20.6 Mbp. This is the region flanked by the BoeWGS1.0 contigs c35489 and c244127 on one side, and c12059 and c12060 on the other (NCBI accessions CCJE01069526.1, CCJE01173868.1, CCJE01026433.1 and CCJE01045859.1).

Genome Sequencing
For the development of molecular markers and to gain knowledge about genome regions as well as loci relevant for grapevine breeding, we analyzed the genome of the resistant interspecific hybrid 'Börner' (V. riparia x V. cinerea) by WGS sequencing. The draft assembly BoeWGS1.0 covers a total of 618 Mbp, which corresponds to 62% of the expected 1 Gbp long diploid 'Börner' genome. Assessment by REAPR revealed some critical regions that could be caused by the incompletely resolved haplophases. Mapping of reads to at least in parts very similar sequences, such as separated alleles, still poses a challenge. However, most of the contigs pass the read mapping-based assembly evaluation.
A major challenge during the de novo assembly of this dataset was to optimize the separation of sequences into contigs derived from both haplophases for subsequent generation of highly informative SSR and SNV marker assays. Different evaluations for the success of this separation result in somehow related-but in detail different-numbers. BUSCO finds about 14% duplicated benchmarking genes. Some genes in the reference set are,  however, not detected at all (BUSCO completeness is about 75%). This is probably an underestimation for the contig level since the relatively small contigs complicate gene calling. The total length of the assembly is about 120 Mbp longer than the haploid Vitis genome sequence, indicating that 240 Mbp of the 618 Mbp represent sequences from separated haplophases. The analyses of contig mappings against the reference genome sequence PN40024 indicates that 142 Mbp of the reference are covered by two BoeWGS1.0 contigs. Taken together, we assumed that the BoeWGS1.0 assembly contains phase-separated sequence information for one quarter of the 'Börner' genome. Note that this information is not "phased throughout" since, even if a contig pair is detected, we lack information about which of the two belongs to the V. cinerea and which belongs to the V. riparia haplophase. Anyway, since the regions for which separation worked are distributed throughout the genome, the data are a very good source for marker development.
Although the separation of sequences from both haplophases is great for molecular breeding applications and several genomewide investigations, the short-read-derived BoeWGS1.0 assembly is still a draft genome sequence. The short contig length significantly limits downstream analyses like gene prediction or approaches to detect genomic rearrangements. The main reason for the fragmentation of the BoeWGS1.0 assembly are regions where the ancestor genomes of V. cinerea and V. riparia are quite similar. In fact, the frequency of SNVs and MNVs detected after read mapping between the two haplophases of 'Börner' is in the same range as between 'Börner' and PN40024. However, it should be considered that read mapping requires quite similar sequences to be specific.
Currently, substantially longer read lengths (mean >20 kbp) than those that were used here are provided by third-generation single-molecule sequencing technologies like those offered by Pacific Biosciences (PacBio) or Oxford Nanopore Technologies (ONT). The long reads can be used to solve the problem of bridging repeat sequences or transposable elements. The first de novo assembly of domesticated grapevine using PacBio sequencing was from V. vinifera cv. 'Cabernet Sauvignon' (Chin et al., 2016) followed recently by the diploid assembly of the cultivar 'Chardonnay' (Roach et al., 2018;Zhou et al., 2019). The first and only wild Vitis genome sequence so far comes from the rootstock Vitis riparia 'Gloire de Montpellier' and reached considerable continuity by incorporating PacBio and 10X Chromium (Girollet et al., 2019). However, assembling a highly heterozygous genotype, accurate phasing, and aligning of Vitis haplotypes is still challenging. With regard to the data on hand and the complexity of this assembly, the application of dedicated scaffolding tools like scaffoldScaffolder (Bodily et al., 2015) or SSPACE (Boetzer et al., 2011) comes at a high risk of generating erroneous connections between contigs. Therefore, we calculated the BoeWGS1.0 assembly with very stringent parameters using the CLC Genomics Workbench toolkit, which contains an implementation of a de Bruijn Graph assembler, and refrained from scaffolding.
In addition to the room for improvement for the assembly parameters, there is a significant fraction of the BoeWGS1.0 genome sequence that is currently not able to align, and it is thus not assignable to the reference genome sequence. This genome fraction likely holds sequence information underlying important traits. Chromosomal anchoring of these contigs is a task that could be realized by the use of the provided SSR or SNV marker sequences of this study in an existing F1 mapping population (V3125 x 'Börner') . Finally, the generation and incorporation of long read data for improvement of the 'Börner' genome sequence is an option for future studies. We expect that a long-read assembly with good continuity would also allow to study 'Börner' sequences that display high divergence from the reference.

Deduction of SSR Sequences for Marker Development
The detected polymorphic SSR candidate loci are a very valuable resource for genome-wide marker development and genotyping. SSR markers are still the most abundant molecular marker type used for mapping in Vitis species and marker-assisted selection in grapevine breeding programs. Because of their very polymorphic nature, SSRs often allow us to clearly distinguish between more than two alleles. This is highly necessary when using F1 mapping populations derived from a cross of highly heterozygous parents, when following an allele in a phylogenetic tree, or when identifying accessions (Cipriani et al., 2011). All these conditions fit to grapevine.
In this study, we were able to reliably predict thousands of candidates for SSR markers. Furthermore, we validated high prediction reliability, which is essentially limited only by the design of specific primers, and demonstrated their applicability in Vitis genotyping approaches. The prediction schema for selection of promising SSR candidates that relied on aligning contigs representing separated haplophases instead of large amounts of sequence reads turned out to be very successful. We were able to reliably detect, predict, and select for longer sequence variations caused by the variable number of repeated units. This is favorable for marker detection by fragment analysis as the current predominant method used for SSRbased genotyping.

Heterozygosity of 'Börner'
For V. riparia (cv. 'Gloire de Montpellier'), the mean distance between heterozygous SNVs is 217 bp (Girollet et al., 2019). For the heterozygous 'Pinot noir,' a frequency of 1 nucleotide variant per 100 bp and 1 InDel per 0.45 kbp was described (Velasco et al., 2007). The overall variant frequency that we obtained between the 'Börner' genome sequence and the reference PN40024 was 1 variant per 68 bp. This number is likely to be a significant underestimation because of the stringent read mapping parameters and the narrow confines for the coverage filter. A calculation of heterozygosity using GenomeScope indicated a heterozygosity of 3.1% for 'Börner' and 1.5% for 'Pinot noir' as an example for a prominent V. vinifera cultivar. This result gives additional indication for the highly heterozygous nature of the 'Börner' genome compared to V. vinifera varieties.

Application of Phased Contig Data for Targeted Mapping of the Downy Mildew Resistance Rpv14 of Vitis
Single SNVs or groups of nearby SNVs within the multiple alignments of two 'Börner' contig pairs and the PN40024 reference sequence were shown to be useful for the molecular identification and localization of a quantitative trait locus. The highly reliable SNVs detected in our study can function as molecular markers in a bulked segregant analysis as shown for Rpv14 on chr.5 in our analysis. Rpv14 was selected as an example for proof of concept because the alternative Rdv1, which causes phylloxera resistance, has already been restricted to an area of about 350 kbp-a value deduced from the physical distance of the published flanking genetic markers (Zhang et al., 2009;Hausmann et al., 2014) relative to the PN40024 reference sequence. The other option, namely black rot resistance mediated by Rgb1/Rgb2, was not selected because there are at least two segregating loci (Rex et al., 2014). We were able to physically reduce the size of the Rpv14 target region to less than 500 kbp. In relation to the results from Ochssner et al. (Ochssner et al., 2016), the center of the target region is narrowed down from the north and the south with the center being still at 20.06 Mbp. One limitation for further reduction of the target region was the number of F1 individuals in our R and S pools or BSA bulks. A recent review on BSA (Zou et al., 2016) mentioned an optimal bulk size for monogenic traits of 10 to 20% of the segregating population for each pool or "tail" of the trait value distribution. However, a BSA study in rice (Wambugu et al., 2018) was successful with 10 or fewer individuals per pool. In addition, equal size of the contrasting pools is considered to be important. In our study, the pool size of 10 resulted from the number of F1 genotypes with reliable and consistent phenotypic scoring results within the population. Obviously, this limits the resolution since only a relatively small number of recombination events were evaluated.
Further research should focus on the identification of new F1 individuals with additional recombinations in the target region. This will probably require enlargement of the number of F1 individuals in the mapping population and additional phenotyping for downy mildew resistance. Subsequently, the work presented here allows for quick access to additional markers, which can be used to determine the recombination points that then can be correlated with resistance genes from candidate gene predictions.

AUTHOR CONTRIBUTIONS
DH and LH conceived and planned the experiments. DH, PV, and LH designed and performed the experiments. TR, DH, and BP calculated the data. RT, BW, and DH conceived the original idea and supervised the project. DH, RT, and BW wrote the manuscript with input from all authors. All authors interpreted and discussed results.

FUNDING
The project was supported by the German Ministry of Education and Science (BMBF) grant FKZ 0315460A and 0315460B (Acronym: GrapeReSeq) to BW and RT, respectively. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.