Genome-wide impacts of alien chromatin introgression on wheat gene transcriptions

Agronomic characteristics and tolerance to biotic and abiotic stresses in hexaploid wheat can be drastically improved through wheat-alien introgression. However, the transcriptional level interactions of introduced alien genes in the wheat genetic background is rarely investigated. In this study, we report the genome-wide impacts of introgressed chromosomes derived from Ae. longissima on gene transcriptions of the wheat landrace Chinese Spring. RNA-seq analyses demonstrated 5.37% and 4.30% of the genes were significantly differentially expressed (DEGs) in CS-Ae. longissima disomic 3Sl#2(3B) substitution line TA3575 and disomic 6Sl#3 addition line TA7548, respectively when compared to CS. In addition, 561 DEGs, including 413 up-regulated and 148 down-regulated or not transcribed genes, were simultaneously impacted by introgressed chromosomes 3Sl#2 and 6Sl#3, which accounts for 41.25% of the DEGs in TA3575 and 38.79% in TA7548. Seventeen DEGs, annotated as R genes, were shared by both introgression lines carrying chromosomes 3Sl#2 and 6Sl#3, which confer resistance to powdery mildew. This study will benefit the understanding of the wheat gene responses as result of alien gene(s) or chromosome intogression and the plant defense response initiated by powdery mildew resistance genes in chromosomes 3Sl#2 and 6Sl#3.


Results
Karyotype analysis of TA3575. Karyotype analysis of TA3575 was identified by integration analyses of PCR amplification using ten 3B-specific and 12 3S l -specific markers (Supplementary Table S1) and in situ hybridization of root-tip cell mitotic metaphase chromosomes. The PCR results showed that no 3B-specific bands were amplified using 3B-specific primer pairs in both TA3575 and the N3B-T3D control, where a pair of chromosome 3B was replaced by two pairs of 3D. Conversely the 3S l #2-specific amplifications resulted only from TA3575, indicating the absence of 3B chromosomes and the presence of alien chromosome 3S l #2 in TA3575. The in situ hybridization results confirmed that in TA3575, 20 pairs of wheat chromosomes plus a pair of 3S l #2 chromosomes from Ae. longissima replaced the pair of 3B chromosomes (Fig. 1, Supplementary Fig. S1). Thus, TA3575 was redesignated as a CS-Ae. longissima disomic substitution (DS) 3S l #2(3B) line, rather than the disomic addition (DA) 3S l #2 line previously identified by the Wheat Genetic Resource Center (WGRC), Kansas State University, Manhattan, KS.
RnA-seq quantity analysis, sequence assemble. A total of 77,962,967, 77,492,142 and 76,266,826 clean reads (>100 bp) were obtained from RNA sequencing for CS, TA3575 and TA7548, respectively. Sequencing quality scores Q30, which inferred a base call accuracy of 99.9%, were greater than 93.06% for all samples, which signifies that the RNA-seq quality was appropriate for subsequent sequence assembly and analysis.
Clean reads of all the samples were further assembled into a total of 158,953 unigenes with a median length of 1,247.95 bp using the short-read assembly software, Trinity. The contig N50 value was 1,725 bp in length. We mapped 91.29% clean reads back to the reference assembly.
Further analysis of DEG distribution on wheat chromosomes showed that most DEGs mapped to wheat chromosome 3B; 577 (31.43%) out of the 1,839 genes mapped on chromosome 3B in CS were differentially transcribed in TA3575, with 33.43% on the short and 30.18% on the long arm of chromosome 3B (Fig. 2). Analysis of the expression pattern of these 577 DEGs showed that the majority (461, 79.90%) were non-transcribed and 100 genes (17.33%) were down-regulated. Only 16 (2.77%) DEGs were up-regulated, indicating that at least 34.57% (461/1,839) of genes on the missing chromosome 3B were not genetically compensated by the introgressed 3S l #2 chromosome from Ae. longissima, regardless of gene expression level differences. The DEGs from others chromosome-arms averaged 4.01% (1,413/35,219), ranging from 1.47% (7/475) for the short arm of 5 A (5AS) to 7.23% (86/1,190) for the long arm (5AL) and 7.55% (52/689) for the short arm of chromosome 3D (Fig. 2). Major transcription pattern changes were different for different wheat chromosomes. Besides chromosome 3B, 80% of the DEGs for chromosome 4BS were non-transcribed or down-regulated genes (Fig. 2), whereas, for chromosome 6DS, at least 75% of the DEGs were up-regulated in the presence of chromosome 3S l #2 in TA3575 (Fig. 2).
Validation of differentially transcribed genes in TA3575. We selected nine genes mapped to chromosome 3B and three genes from other chromosomes to validate gene transcription changes in TA3575 using 18S rRNA as a control (Supplementary Table S1). PCR amplification using genomic DNA (gDNA) from control CS, N3B-T3D, and TA3575 as templates showed that all nine 3B genes were present in CS and absent in N3B-T3D, as expected. However, in TA3575, six of the nine 3B-specific genes, including CL22200Contig1, CL27257Contig1, CL24323Contig1, CL29051Contig1, CL53819Contig1 and CL51762Contig1, were amplified and the remaining three (CL40247Contig1, CL35206Contig1 and CL23073Contig1) were absent (Fig. 3, Supplementary Fig. S3).
Based on semi-quantitative RT-PCR, which was performed using nine 3B-specific primer pairs, four (CL24323Contig1, CL40247Contig1, CL35206Contig1, and CL23073Contig1) were not expressed, two (CL22200Contig1 and CL27257Contig1) were down-regulated and transcribed, and three, including CL29051Contig1, CL53819Contig1, and CL51762Contig1, were up-regulated and expressed in TA3575 www.nature.com/scientificreports www.nature.com/scientificreports/ as compared to CS (Fig. 3, Supplementary Fig. S3). On the other hand, genes not located on chromosome 3B, CL53958Contig1 and CL78654Contig1 were down-regulated and expressed, and CL23519Contig1 was up-regulated in TA3575, compared with CS. The results of both PCR amplifications using gDNA and semi-quantitative RT-PCR using RNA as templates were consistent with the RNA-seq data analyses and suggested that RNA-seq analyses in this study were valid and reproducible.

Distribution of DEGs in TA7548.
A total of 1,592 (4.3%) DEGs, including 396 NEGs (24.87%), 526 DRGs (33.04%) and 670 URGs (42.09%), were identified based on a pairwise comparison of expressed genes in CS-Ae. longissima DA 6S l #3 line TA7548 and CS on comparison of 37,058 mapped genes (Supplementary Tables S2, S3). These 1,592 DEGs were also asymmetrically distributed on different chromosomes in TA7548. The least affected genes were those located on the short arm of 5B; only nine out of 506 genes (1.78%) showed a change in transcription in TA7548 compared to that in CS. On the contrary, chromosome 6A of wheat contained the most DEGs (190/1,385, 7.50%), which included 54 (54/190, 28 Table S2). On an average, about 90% of the DEGs were down-regulated or non-transcribed on both the short and the long arms of chromosome 6A, followed by chromosome 4BS, where more than 80% of DEGs were negatively regulated. In contrast, close to 70% of DEGs on the short arm of chromosome 4D were positively impacted by chromosome 6S l #3 in TA7548 (Fig. 4, Supplementary Table S2).
Unlike in TA3575, where the highest ratio of non-transcribed to transcribed genes, were those genes that were located to the missing chromosome 3B. The percentage of DRGs mapped to chromosome 6A in TA 7548 was significantly higher than that for any other chromosomes. The DEGs included 13.76% on the short arm and 13.65% on the long arm of chromosome 6A (Fig. 4). In addition, down-regulated genes (DRGs) were a major DEG type on both arms of chromosome 6A (Fig. 4), indicating that genes on both arms were negatively affected to the some extent by introgressed chromosome 6S l #3 in TA7548.
Transcription-affected genes shared by both introgressed 3S l #2 and 6S l #3 in CS. Based on a pairwise comparison of DEGs between TA3575 and TA7548, 561 DEGs were simultaneously impacted by introgression of both chromosomes 3S l #2 and 6S l #3, which accounted for 41.25% and 38.79% of the total, excluding those on chromosome 3B in TA3575 (1,360) and 6A in TA7548 (1,446  www.nature.com/scientificreports www.nature.com/scientificreports/ The 561 DEGs shared by chromosomes 3S l #2 and 6S l #3 from Ae. longissima were annotated by blastx alignment against the GO database. The results indicated that most of the enriched, shared DEGs were categorized biologically as negative regulators of programmed cell death, followed by those involved in protein phosphorylation, defense response to fungus and response to chitin. Of those categorized as cellular components, cytoplasmic membrane-bounded vesicle and an integral component of the plasma membrane were the top terms. The molecular function category included protein serine/threonine kinase activity and transferase activity (   www.nature.com/scientificreports www.nature.com/scientificreports/ Impact of introgressed chromosomes 3S l #2 and 6S l #3 on putative disease resistance genes of cS. Both chromosomes 3S l #2 and 6S l #3 from Ae. longissima carry resistance gene(s) to powdery mildew of wheat 41 . The impact of 3S l #2 in TA3575 and 6S l #3 in TA7548 on the transcription of CS-derived putative plant disease resistance (R) genes was analyzed in this study. A sequence homology search of DEGs in the PRGdb discovered 34 (1.71%) R genes in seven classes from 1,991 DRGs in TA3575, with the maximum number of R genes in the RLP class (16,47.06%), followed by those in classes NL (8,23.53%), CNL and T (3, 8.82% each), and N, TNL and other (1, 2.94% each) ( Table 2, Supplementary Table S6). Of these 34 putative R genes, 26 (76.47%) were up-regulated and eight (23.53%) down-regulated by chromosome 3S l #2. There were also two up-regulated R genes (RLP class) (CL77976Contig1 and CL29051Contig1) mapped on chromosome 3B (in CS) in TA3575 where chromosome 3B was missing, indicating that these two R genes are derived from 3S l #2 rather than CS wheat chromosome 3B ( Table 2, Supplementary Table S6).
A pairwise comparison of effected R genes between TA3575 and TA7548 further revealed that 17 R genes were impacted together by the introgression of both chromosomes 3S l #2 and 6S l #3. These 17 shared R genes were grouped into six classes, including one R gene for each class CNL, TNL and others, four in NL, seven in RLP and three in class T. Twelve (70.59%) R genes were transcribed with up-regulation, whereas the remaining five (29.41%), which all belong to class RLP, were down-regulated in both TA3575 and TA7548 (Supplementary  Table S6). Transcription patterns of these 17 shared R genes were similar in both TA3575 and TA7548.
Identification of a CS-Ae. longissima DS 6S l #3(6A) line from TA7548. As shown in TA3575, if a pair of wheat chromosomes is missing, non-transcribed genes account for the majority of DEGs located on the substituted chromosomes. In TA3575, where a pair of 3B chromosomes is replaced by a pair of 3S l #2 chromosomes, 461 out of 577 DEGs located to 3B were NEGs, accounting for as high as 79.90%. However, in TA7548 an unusually high ratio of DRGs (119/190, 62.63%) rather than NEGs (54/190, 28.42%) of chromosome 6A-specific DEGs led to the suspicion that chromosomes 6A might be missing in some TA7548 individuals (progenies) used for RNA-seq in this study.
Molecular markers analysis of seven, randomly selected individual plants from TA7548 using five chromosome 6A-specific SSR and five chromosome 6S l -specific markers showed that all these chromosome 6A-specific SSR markers were present in five of the seven plants (5/7, 71.43%); two individuals (2/7, 28.57%) presented 6S l -specific markers but lacked 6A-specific markers (Fig. 6a,b, Supplementary Fig. S4). Subsequent genomic in situ hybridization (GISH) analysis of these individual plants indicated that the two individuals missing 6A-specific markers had 42 chromosomes, 40 wheat chromosomes plus a 6S l chromosome pair. Further fluorescent in situ hybridization (FISH) analysis verified that a pair of 6A chromosomes was replaced by a 6S l #3 chromosome pair in these two plants, which resulted in a disomic 6S l #3(6A) substitutions (Fig. 6c,d). The remaining five plants were confirmed to be disomic 6S l #3 addition lines, which had 42 wheat chromosomes plus a pair of 6S l chromosomes (Fig. 6e,f). Therefore, the TA7548 population used in this study was composed of about 71.43% DA 6S l #3 and 28.57% DS 6S l #3(6A) plants. This novel CS-Ae. longissima DS 6S l #3(6A) line, spontaneously formed during the long process of reproduction of the original DA 6S l #3 plants, was named TA7548-6A. The substitution of chromosome 6A by chromosome 6S l #3 in some individuals of TA7548 led to the unusually high ratio of DRGs of chromosome 6A in TA7548 used for RNA-seq in this study.

Discussion
Linkage drag, caused by alien chromosome segments introduced into wheat genetic backgrounds, has led to unfavorable agronomic and end-use quality traits and, thus, limited the utilization of alien genes in wheat improvement programs. The linkage drag caused by deleterious genes associated with targeted genes had been widely demonstrated 21,50-52 . However, linkage drag that may be attributed to the inability of alien genes to replace missing wheat genes, or the interactions between the introduced alien genes and recipient wheat genes, has, until now, been relatively under-researched 38,53 . In this study, we verified the genome-wide impact on wheat gene expression by the presence of introduced alien chromosomes. Transcription of 1,413 (4.01%) of 35,219 genes mapped to non-3B chromosomes in CS were significantly changed by the substitution of a pair of chromosome 3B by chromosome 3S l #2 derived from Ae. longissima in TA3575. Whereas 1,402 (3.93%) of 35,673 non-6A genes in CS were affected by introgression of chromosome 6S l #3 in TA7548. The ratios of affected genes in this study were higher than the 3% (960/35,301) reported by Rey et al. (2018), investigating the effect of introgression of barley telosomic 7HL in CS 38 . Further analyses of DEGs on different wheat chromosomes demonstrated that the impact of introgressed chromosomes 3S l #2 and 6S l #3 was obviously different for wheat chromosomes. The minimum ratio of DEGs was associated with chromosome 5AS in TA3575 (1.47%) and chromosome 5BS in TA7548 (1.78%); whereas the highest ratios were on chromosomes 3DS in TA3575 (7.55%) and 6DL in TA7548 (6.35%). Those negatively impacted genes, including non-transcribed and down-regulated, accounted for 48.90% (691/1,413) in TA3575, and 53.42% (749/1,402) in TA7548. Whether these impacts may lead to linkage drag of agronomic traits, or the effects can be decreased in different recipient wheat backgrounds, are currently under investigation. The present study will benefit the understanding of gene interactions between recipient wheat and the alien donor species. www.nature.com/scientificreports www.nature.com/scientificreports/ Of those DEGs caused by introgressed chromosomes 3S l #2 and 6S l #3 from Ae. longissima, 17 genes in six R gene classes were annotated as putative R genes from 561 DEGs shared by introgressed 3S l #2 in TA3575 and 6S l #3 in TA7548. Transcription analysis of these 17 R genes, based on transcriptome analyses, showed that 12 shared R genes (70.59%) were significantly up-regulated and transcribed in both TA3575 and TA7548, whereas only five genes (29.41%) in class RLP were down-regulated. Because both chromosomes 3S l #2 and 6S l #3 confer resistance to powdery mildew of wheat 41 , whether the transcription changes of these shared R genes were involved in the pathogen-defense processes deserved to be further investigated. Molecular marker could be developed from these R gene for mapping of the two powdery mildew resistance genes introgressed from Ae. longissima.
RNA-sequencing techniques have provided a very useful means for examining chromosome structural changes 38,54,55 . In this study, the discovery of an abnormally high percentage of DRGs on both the short and the long arm of wheat chromosome 6A, based on RNA-seq analyses, implied that chromosome variation of chromosome 6A might exist in the CS-Ae. longissima DA 6S l #3 line TA7548. We further verified that 28.57% of the plants in the TA7548 population were, in fact, spontaneously formed DS 6S l #3(6A) by integration analyses of RNA-seq data, molecular markers and in situ hybridization. This novel DS 6S l #3(6A) line (TA7548-6A) will be useful for further transfer of 6S l #3 resistance genes to powdery mildew into wheat by inducing wheat-Ae. longissima 6S l -6A translocations. In a similar way, we selected a spontaneous Robertsonian translocation T4S l #2·4BL carrying the powdery mildew resistance gene Pm66 56 . In summary, integrated use of genome sequencing, molecular markers, and classic cytogenetic techniques can speed up the introgression of alien targeted genes, thus promoting the utilization of wild relatives of wheat in breeding programs.

conclusion
In the current study, transcriptome analyses were performed using two CS-Ae. longissima introgression lines carrying powdery mildew resistance gene(s) and their genetic background line wheat CS. The results showed that introgression of chromosomes 3S l #2 and 6S l #3 derived from Ae. longissima had genome-wide impact on gene transcription in wheat. A total of 5.37% and 4.30% genes (or 4.01% of non-3B-specific genes and 3.93% of non-6A-specific genes) were differentially transcribed due to the introduction of chromosomes 3S l #2 in TA3575 and  www.nature.com/scientificreports www.nature.com/scientificreports/ 6S l #3 in TA7548, respectively. Furthermore, 17 putative R genes of wheat were significantly impacted together in both TA3575 and TA7548, which carry chromosomes 3S l #2 and 6S l #3. The majority of these shared R genes (70.59%) were significantly up-regulated by introgression of chromosomes 3S l #2 and 6S l #3. Whether these shared putative R genes were involved in the defensive processes initiated by powdery mildew resistance genes located on chromosomes 3S l #2 and 6S l #3 needs to be further investigated.

Materials and methods
plant materials. Wheat landrace Chinese Spring (CS); nullisomic 3B-tetrasomic 3D (N3B-T3D, TA3272), where a pair of chromosome 3B is replaced by two pairs of 3D; N6A-T6B (TA3152) where 6A is substituted by two pairs of 6B; a CS-Ae. longissima disomic 3S l #2(3B) substitution line TA3575, where a pair of 3S l #2 chromosomes replaces the 3B chromosome pair in CS; and a CS-Ae. longissima disomic 6S l #3 addition line TA7548, where a pair of 6S l #3 chromosomes are added to the chromosome complement of CS; were used in this study. Both chromosomes 3S l #2 and 6S l #3 confer resistance to powdery mildew 41 . The number following the chromosome designation is used to distinguish between the same Ae. longissima chromosome derived from different Ae. longissima accessions 57 . All plant materials were kindly provided by the Wheat Genetic Resources Center (WGRC) at Kansas State University, USA, propagated and stored at the experimental station of Henan Agricultural University. plant growth and tissue collection. Seeds of the two wheat-Ae. longissima introgression lines (TA3575 and TA7548) and CS were sterilized with 5% sodium hypochlorite solution for 10 min at room temperature and then planted in 8-cm diameter pots filled with an autoclaved vermiculite potting medium. The seedlings were grown in an illuminated incubator at 18-20 °C, 18 h light and 6 h dark, and 75% relative humidity. Pooled, equal-length segments from the first leaf of 10 individuals for each line were then immediately frozen in liquid nitrogen for subsequent RNA isolation. Two independent, biological replicates were collected for subsequent cDNA library construction and RNA sequencing.  www.nature.com/scientificreports www.nature.com/scientificreports/ were sequenced on the Illumina sequencing platform (HiSeqTM 2500 or Illumina HiSeq X Ten) by OE Biotech (Shanghai), China, and 125 bp/150 bp paired-end reads (raw reads) were generated.
RnA-seq data analysis. Paired-end reads were pretreated using the NGS QC Toolkit software to remove sequences containing adapters or poly-N above 5%, and low-quality reads to produce valid data (clean data) for downstream sequence assembly 58 . Valid ratio, Q30 and GC content of each sample were calculated after pretreatment.
De novo assembly of clean data was applied using Trinity (version: trinityrnaseq_r20131110) 59 , followed by the removal of abundant sequences using TGICL 60 . Annotation and function of these genes were assigned with a cut off e-value <1e −5 based on blastx alignment against protein sequences in public databases NR (NCBI non-reduction protein sequences, ftp://ftp.ncbi.nih.gov/blast/db), GO (Gene Ontology, http://www.geneontology.org) and PRGdb (Plant Resistance Gene Database) 61 .
For DEGs analysis, the read counts of genes were first normalized as FPKM. Then, the expression significance of genes of CS vs. TA3575 and CS vs. TA7548 were calculated using DESeq (http://bioconductor.org/ packages/release/bioc/html/DESeq.html) with a threshold FDR ≤ 0.05 & |log 2 (fold change)| ≥ 1. Genes with an FPKM ≥ 0.5 were considered expressed, whereas genes with FPKM = 0 were not transcribed. Genes of CS were pair-wise compared with those of TA3575 and TA7548 to identify genes present in CS and absent in TA3575 and TA7548. All the genes were assigned to chromosome arms based on blastn alignment against wheat reference genomic sequences (CS RefSeq v1.0) 62 with a cutoff of expect < =1e-10&qcov > =75% at UGRI BLAST (https:// urgi.versailles.inra.fr/blast/blast.php). In total, 44 molecular markers were used in the study (Supplementary Table S1). Of these, 10 3B-specific SSR and 12 3S l #2-specific EST markers were used for the karyotype analysis of TA3575 63,64 . Another nine PCR primer pairs of 3B-specific and three non-3B chromosomes were designed based on transcription sequences of CS and used for validation of gene expression changes in TA3575. The remaining five 6A-specific SSR and five 6S l #3-specific markers were used for molecular marker analyses of TA7548 63,65-67 . PCR reaction mixture preparation and amplification by "F50SSR" (for SSR markers) or "Touch-down 63" (for the remaining markers) followed Liu et al. 14 . The PCR products were digested with restriction enzymes followed by Liu et al. 14 . The PCR products or restricted fragments were resolved in 2.5% agarose gels and visualized by ethidium bromide staining under UV light. chromosome preparation and in situ hybridization. Root tip collection, nitrous oxide treatment, fixation, and slide preparations were according to Liu et al. 14 . GISH using fluorescein-12-dUTP labeled genomic DNA of Ae. longissima as the probe, was followed the procedure of Liu et al. 14 . The ratio of genomic Ae. longissima DNA and CS blocking DNA was 1:120. FISH using FAM (6-carboxyfuorescein) and TAMRA (6-carboxytetramethylrhodamine) modified oligonucleotide probes followed Huang et al. 68 . After hybridization and slide washing, 25-30 µl of Vectashield mounting medium containing 1 µg/ml DAPI (Vector Laboratories Inc, Burlingame, CA, USA) was added to each slide and then covered with a 24 × 30 cm glass coverslip. Observation of fluorescent images was on a Zeiss Axio Scope A1 fluorescence microscope (Germany) and captured with an AxioCam MRc5 CCD camera. Images were further processed with Adobe Photoshop CS3 (Version 10.0.1) (Adobe Systems Inc., San Jose, CA, USA).