Genomic Characteristics of Elite Maize Inbred Line 18-599 and Its Transcriptional Response to Drought and Low-Temperature Stresses

Elite inbred line 18-599 was developed via triple test cross from introduced hybrid P78599 and used as parents of dozens of maize hybrids adapting to the diverse ecological conditions of the maize ecological region in Southwest China. In this study, its genomic DNA was resequenced and aligned with the B73 genome sequence to identify single nucleotide polymorphism (SNP), and insertion (In) and deletion (Del) loci. These loci were aligned with those between B73 and 1020 inbred lines in the HapMap database to identify specific variation loci of 18-599. The results showed that there were 930,439 specific SNPs and 358,750 InDels between 18-599 and the 1020 lines. In total, 21,961 of them showed significant impacts on the functions of 12,297 genes, such as frameshift, change of splicing site, stop gain, change of start site, and stop loss. Phylogenetic analysis showed that 18-599 was closely related to inbred lines ZEAxujRAUDIAAPE and 2005-4, but far from some inbred lines directly isolated from P78599. This result indicated that 18-599 not only pyramided the elite genes of P78599, but also acquired genetic divergence during the repetitive backcrosses of triple test cross to confer its elite agronomic characteristics. Subsequently, the RNA of 18-599 was sequenced. The aligned 9713 and 37,528 of the 165,098 unigenes were screened and aligned with annotated transcripts of the B73 genome differentially expressed under drought and low-temperature stress, respectively, and their functions were involved in the responses to these stresses. The quantitative PCR results of fourteen random genes verified the RNA sequencing results. These findings suggest that the transcriptional responses of many resistance-related genes were an important mechanism for 18-599 to adapt to diverse ecological conditions.


Introduction
The conventional method of maize breeding is to isolate homozygous plants from germplasm resources and develop inbred lines with elite agronomic characteristics, wide adaptability, good quality, and high reproductive capacity. These inbred lines are evaluated for their combining ability by test cross and used as parents of hybrids. After evaluation for their high yield, adaptability, and quality by multiyear and multilocation experiments, superior hybrids are approved for commercial dissemination [1][2][3][4]. Therefore, the germplasm resources present the most fundamental source of elite agronomic characteristics of these superior hybrids.
Maize is originally an allogamous plant. After long-term domestication and selection under natural and cultivation conditions, abundant genetic divergences have accumulated, and a variety of germplasm resources with different characteristics have been inherited. However, the breeding and dissemination of hybrids between inbred lines have failed to protect germplasm resources while doubling maize yield. As a result, the germplasm

Specific Variation Loci of 18-599
Of the 4,125,020 SNPs and 588,952 InDels called by the alignment between ge sequences of inbred lines 18-599 and B73, 3,213,863 SNPs and 230,202 InDels matched with the allelic variation loci aligned between B73 and the 1020 inbred li the HapMap database (v3.21). The other 911,157 SNPs and 358,750 InDels were iden as non-allelic specific loci of 18-599. Of the matched SNPs and InDels, the genoty 19,282 loci were found to be different from those aligned between B73 and the 1020 i lines and identified as specific genotypes of 18-599. Of the total 930,439 specific SNP 358,750 specific InDels between 18-599 and the 1020 lines, 21,961 loci showed sign impacts on the functions of 12,297 genes, such as frameshifts, change of splicing site gains, change of start sites, and stop losses (Table 1).

Identical Variation Loci and Phylogenetic Tree
The alignment showed that 18-599 shared 30-40% identical variation loci with the 1020 inbred lines, and more than 50% identical variation loci with 31 of them 2). After removing the deletion and heterozygous genotypes, these 31 inbred lines used to construct a phylogenetic tree with 18-599 ( Figure 2). The result showed th 599 was closely related to inbred lines ZEAxujRAUDIAAPE and 2005-4, and CA

Specific Variation Loci of 18-599
Of the 4,125,020 SNPs and 588,952 InDels called by the alignment between genomic sequences of inbred lines 18-599 and B73, 3,213,863 SNPs and 230,202 InDels were matched with the allelic variation loci aligned between B73 and the 1020 inbred lines in the HapMap database (v3.21). The other 911,157 SNPs and 358,750 InDels were identified as non-allelic specific loci of 18-599. Of the matched SNPs and InDels, the genotypes of 19,282 loci were found to be different from those aligned between B73 and the 1020 inbred lines and identified as specific genotypes of 18-599. Of the total 930,439 specific SNPs and 358,750 specific InDels between 18-599 and the 1020 lines, 21,961 loci showed significant impacts on the functions of 12,297 genes, such as frameshifts, change of splicing sites, stop gains, change of start sites, and stop losses (Table 1).

Identical Variation Loci and Phylogenetic Tree
The alignment showed that 18-599 shared 30-40% identical variation loci with 743 of the 1020 inbred lines, and more than 50% identical variation loci with 31 of them (Table 2). After removing the deletion and heterozygous genotypes, these 31 inbred lines were used to construct a phylogenetic tree with 18-599 ( Figure 2). The result showed that 18-599 was closely related to inbred lines ZEAxujRAUDIAAPE and 2005-4, and CAU 178 isolated from hybrid P78599, as well as 78,599 itself, but situated far from inbred lines Cheng 18, Dan 599, R150, and other inbred lines isolated from hybrids P78599 and P78641.    ,602 clean reads were retained, respectively, accounting for more than 96% of the raw reads. About 70% of the clean reads were mapped to the B73 genomic sequence (Table S1). From the unmapped, a total of 165,098 unigenes with a total length of 74,602,848 bp and a contig N50 of 496 bp were assembled by TRINITY.

Transcriptional Response to Drought and Low-Temperature Stress
By using Stringtie and DESeq v2. 2, 12,914 differentially expressed genes (DEGs) were identified in response to drought stress (Table S2). Among them, the expressions of 6630 genes were up-regulated with an average multiple (log 2 ) of 2.33 times. The upregulated multiple (log 2 ) of gene Zm00001d015053 was as high as 24.35 times. The expressions of 6374 genes were down-regulated with an average multiple (log 2 ) of 2.44 times. The down-regulated multiple (log 2 ) of gene TRINITY_ DN38188_ c10_ g1 was as low as 21.18 times. A total of 46,698 differentially expressed unigenes were identified in response to low-temperature stress (Table S3). Among them, the expressions of 21,532 genes were up-regulated, with an average multiple (log 2 ) of 2.48 times. The upregulated multiple (log 2 ) of gene TRINITY_DN75068_c0_g1 was as high as 22.61 times. The expressions of 25,252 genes were down-regulated, with an average multiple (log 2 ) of 2.87 times. The down-regulated multiple (log 2 ) of gene Zm00001d038996 was as low as 143.89 times.

RT-qPCR Verification
The RT-qPCR results showed that the relative expression levels of genes TRINITY_ DN36046_c1_g3, TRINITY_DN13689_c0_g1 and TTRINITY_DN39631_c1_g1 were significantly down-regulated, and those of genes TRINITY_DN33619_c1_g1, TRINITY_DN72435_ c0_g1, TRINITY_DN20047_c0_g1 and TRINITY_DN16287_c0_g1 were significantly upregulated in response to drought stress ( Figure 3a). In response to low-temperature stress, the relative expression levels of genes TRINITY_DN43011_c0_g1, Zm00001d043044, and TRINITY_DN63145_c0_g1 were significantly down-regulated, and those of genes TRIN-ITY_DN73242_c0_g1, TRINITY_DN20047_c0_g1, Zm00001d053091, and TRINITY_DN17398_ c0_g1 were significantly up-regulated ( Figure 3b). These results verify the above transcriptomics analysis results.
ITY_DN72435_c0_g1, TRINITY_DN20047_c0_g1 and TRINITY_DN16287_c0_g1 were significantly up-regulated in response to drought stress ( Figure 3a). In response to low-temperature stress, the relative expression levels of genes TRINITY_DN43011_c0_g1, Zm00001d043044, and TRINITY_DN63145_c0_g1 were significantly down-regulated, and those of genes TRINITY_DN73242_c0_g1, TRINITY_DN20047_c0_g1, Zm00001d053091, and TRINITY_DN17398_c0_g1 were significantly up-regulated (Figure 3b). These results verify the above transcriptomics analysis results.
(a) (b) Figure 3. The relative expression levels of fourteen genes in response to drought and low-temperature stress. (a) Seven genes in response to drought stress; (b) Seven genes in response to low-temperature stress.

Discussion
The resequencing and genomics analysis showed that elite inbred line 18-599 had 21,961 specific variation loci in 12,297 genes while aligned against the 1020 inbred lines in the HapMap database ( Table 1). The phylogenetic analysis showed that 18-599 was closely related to ZEAxujRAUDIAAPE and 2005-4, as well as CAU 178 isolated from hybrid P78599, but situated far from 78599, Cheng18, Dan 599, R150, and other inbred lines

Discussion
The resequencing and genomics analysis showed that elite inbred line 18-599 had 21,961 specific variation loci in 12,297 genes while aligned against the 1020 inbred lines in the HapMap database (  Figure 2). This result indicated that 18-599 was isolated from the breeding population constructed with hybrids P78599 and P78641 by the triple test cross method and evaluated under diverse ecological conditions [17], which not only inherited some elite germplasm of P78599 and P78641 but also accumulated abundant genetic variation with 78599, Cheng 18, Dan 599, R150 and some other lines isolated from the same hybrids by conventional methods, involving many genes responding to drought and low-temperature stress (Figures 3 and 4). The elite genotypes of these genes were pyramided by the repetitive test crosses and backcrosses of the triple test cross and conferred high combining ability, high reproduction capacity, strong resistance, extensive adaptability, and elite agronomic characteristics of 18-599 [18,19].
The expressions of 37,528 genes were responsive to low-temperature stress. Of these, 1394 were annotated and involved in biological processes (845), molecular functions (368), and cell components (181). The function classification of 65 was responsive to abiotic stimulus, 22 responsive to temperature stimulus, and 13 responsive to cold stress ( Figure 5). Among them, the expression of gene TRINITY_DN75068_c0_g1 was upregulated 22.61 times. Its function has not been annotated. The expression of gene Zm00001d038996 was downregulated 143.89 times. It encodes glycosyltransferase 1. The relative expression levels of seven randomly sampled genes detected by RT-qPCR verified the results of RNA-Seq and transcriptomics analysis (Figure 3b). The three down-regulated genes TRIN-ITY_DN43011_c0_g1, Zm00001d043044, and TRINITY_DN63145_c0_g1 encode a chlorophyll A/B binding protein, an autophagy protein (Cost1), and MYB-like DNA-binding domain, respectively. The up-regulated genes TRINITY_DN20047_c0_g1, Zm00001d053091, and TRINITY_DN17398_c0_g1 encode an annexin, a flavin-binding kelch domain F box protein (FKF2) and a heat shock factor, respectively, except gene TRINITY_DN73242_c0_g1 which encodes an unknown protein. The orthologs of the seven differentially expressed genes with function annotations were well documented for their response to extreme temperatures and other abiotic stresses [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48].
Because drought and low-temperature stresses are major environmental factors that restrict maize growth and lead to yield decrease, the above results suggested that the transcriptional responses of many resistance-related genes were important mechanisms for 18-599 to adapt to the diverse ecological conditions in Southwest China [49,50].

Identification of Specific Variation Loci and Phylogenetic Analysis
A Perl script was used to identify specific SNPs and InDels of 18-599 by alignment of the SNPs and InDels between 18-599 and the B73 genomic sequence, and those between the 1020 inbred lines in HapMap database v3.21 (https://www.sanger.ac.uk/resources/ downloads/human/hapmap3.html, accessed on 20 April 2020) and the B73 genomic sequence [56]. The functions of the genes involving these SNPs and InDels were annotated by ANNOVAR (https://doc-openbio.readthedocs.io/projects/annovar/en/latest/, accessed on 10 May 2020). The genes with variation loci of significant impact on their function, such as frameshift, change of splicing site, stop gain, change of start site, and stop loss, were filtered. Inbred lines with an identity of SNPs and InDels more than 50% with 18-599 were filtered from the HapMap database v3.21. After removing deletion and heterozygous genotypes, they were used to construct a phylogenetic tree with 18-599 by PHYLIP (https://evolution.genetics.washington.edu/phylip.html, accessed on 16 June 2020).

RNA Preparation and Sequencing
In total, three independent biological replicates of 18-599 seedlings were stressed under simulative drought conditions of 16% (Polyethylene glycol, PEG), low temperatures of 8 • C 14 h/4 • C 10 h, and 28 • C and 300 mmol/m 2 ·s (blank control), respectively. Three days later, each of the replicates was ground in liquid nitrogen under RNase-free conditions, and used for RNA extraction by using RNAiso Plus (TaKaRa, Dalian, China). Possible contamination of genomic DNA was removed by gDNA Eraser Kit (TaKaRa, Japan). After detection of concentration in Qubit 2.0 fluorometer (Invitrogen, USA), OD260/OD280 ratio in Nanodrop spectrophotometer (Thermo Fisher Scientific, USA), and integrity by 0.8% agarose gel electrophoresis, the qualified samples were sequenced on an Illumina Hiseq xten platform at Meiyin Gene (China). Evaluation of sequencing quality and filtration of clean reads were similar as described above. The clean reads were mapped against the maize reference genome sequence of B73 (V4), indexed by SAMtools (http://samtools.sourceforge. net/, accessed on 18 September 2020), and stored in BAM format. The unmapped reads were de novo assembled by TRINITY (https://www.plob.org/tag/trinity, accessed on 5 January 2021) [57]. The longest assembled result of the same TRINITY clustering was filtered by a Perl script (Supplement S1). The unigenes were indexed by HISAT2 v 2.1.0 (https://www.biostars.org/p/288726/, accessed on 12 March 2021), integrated with GTF annotation file of B73 genome (Zea_mays.AGPv4.38.gtf). The unmapped reads of each individual were mapped onto unigene sequences. The result of alignment was merged with the alignment of B73 for each sample.

RT-qPCR Verification
To verify the differential expression displayed by the above transcriptomics analysis, fifteen pairs of primers (Table S4) were designed by Primer-BLAST (https://www.ncbi. nlm.nih.gov/tools/primer-blast, accessed on 15 November 2021) and used for RT-qPCR amplification of fourteen random DEGs as well as internal control gene ZmGAPDH. The three replicates of the RNA samples prepared above were reversely transcribed into cDNA with TransStartR Tip Green qPCR SuperMix (Transgen Biotch, Beijing, China), and used as templates for two-step amplification of RT-qPCR by using SYBR ® Premix Ex Taq TM II (Tli RNaseH Plus) (TaKaRa, Dalian, China) in CFX96 TM Real-Time System (Bio-Rad, Hercules, CA, USA). The temperature procedure was as follows: 95 • C for 30 s; 39 cycles of 95 • C for 5 s, and 50-57 • C for 30 s. At the end of the last cycle, the temperature was increased to 95 • C by 0.5 • C/s, so that a melting curve could be calculated and used to differentiate between specific and non-specific amplicons. The 2 −∆∆CT method of the CFX Manger TM software version 2.0 (Bio-Rad, Hercules, CA, USA) was used to normalize the expression differentiation between the sampled and the internal control genes [62]. The average relative expression levels of four technical and three biological replicates were calculated by their comparison to that of the internal control gene. The statistical significance was assessed via Student's t-test with IBM-SPSS software (http: //www-01.ibm.com/software/analytics/spss/, accessed on 26 February 2022).

Conclusions
Inbred line 18-599 not only pyramided the elite genes of P78599 but also acquired genetic divergence during the repetitive backcrosses of triple test cross, conferring its elite agronomic characteristics. The transcriptional responses of many resistance-related genes were important mechanisms for 18-599 to adapt to the diverse ecological conditions in Southwest China. However, a large number of the specific variations were found during the responsive DEGs and more in others. It was difficult to identify the direct connections between the specific variations and the responsive DEGs. The elite agronomic and adaptive traits remain to be elucidated, although they were thoroughly evaluated during the production practice. Ultimately, the results of this study still provide a reference for future research and breeding.