Qiqiang Guo·Huie Li·Zengqiang Qian·Jie Lu·Weilie Zheng

Abstract Five Larix species (L.griffi thii,L.speciose,L.himalaica,L.kongboensis,and L.potaninii var.australis),have survived on the Qinghai-Tibet Plateau (QTP) under specific climate conditions for decades.The lack of genomic information seriously hinders research on the evolution,conservation and ecology of these Larix resources.In this study,complete chloroplast (cp) genomes of the 5 species were assembled and compared based on next generation sequencing technology combined with polymerase chain reaction validation.The results show that the 5 cp genomes are relatively conservative in size,gene content and arrangement,and border variation.Phylogenetic analysis showed that the species are closely related as well as to seven other species of the same genus.In addition,the 5 cp genomes contained few simple sequence repeats and relatively low nucleotide variability;thus,12 candidate polymorphic cp DNA markers will be helpful for further research on relevant population genetics.These results will provide valuable genetic information for the conservation,evolution and ecology of these species and their relatives.

Keywords Larix species·Plastid genome ·Microsatellite·Phylogenetic analysis·Qinghai-tibet plateau

Introduction

The genusLarix(family Pinaceae),commonly known as larch,represents a group of sun-loving,cold-tolerant deciduous conifers.Due to such characteristics as strong timber,f ine structure and resistance to decay,Larixspecies are widely used for papermaking,buildings and bridge construction,and furniture-making.The species of this genus are generally 10—20 m tall and straight with symmetrical crowns,well-developed root systems.Therefore,they are also an excellent landscaping species.Local residents in the Qinghai-Tibet Plateau (QTP) have traditionally used the darkred corewood for making barrels for ghee,a type of cooking oil.

The QTP in southwestern China has high altitudes and low temperatures.As one of the global biodiversity hotspots,the Plateau’s uplift has resulted in many isolated habitats,a prerequisite for the development of numerous species through speciation.Larixspecies are widely distributed in North America and Eurasia in the temperate and cold temperate high-altitude mountains;the QTP is the southernmost boundary of their distribution.Five species,Larix griffi thiiHook.f.,L.speciose,L.himalaicaW.C.Cheng &L.K.Fu,L.kongboensisR.R.Mill,andL.potaniniivar.australisA.Henry ex Handel-Mazzetti,have survived in the QTP for decades.However,morphology-based taxonomy remains controversial,especially as to the interrelationship ofL.griff ithii,L.kongboensis,andL.potaniniivar.australis(Liu et al.1978;Fu 1983;Fu et al.1999).

Modern taxonomy and related studies based on genomic information have significant advantages such as simplicity,rapidity,accuracy,and independence from tissue types,developmental stages and environmental conditions.The complete chloroplast (cp) genomes ofLarixpossess unique and strict paternal genetic characteristics through pollen(Bondar et al.2019).To date,cp genomes of some species and their relatives in the family have been deciphered and have been applied to studies of their evolution,conservation and ecology (Lin et al.2010;Wu et al.2011;Hao et al.2016;Han et al.2017;Ishizuka et al.2017;Qiu et al.2018;Bondar et al.2019).However,little is known about these fiveLarixspecies native to the QTP,and this has seriously hindered these studies.To address this gap,cp genomes of these five species were assembled and characterized to provide valuable genomic information for their conservation,evolution and ecology.

Materials and methods

Plant materials and sequencing

Fresh needle leaves of the 5 species were collected from Tibet between August and October 2017.The specific locations are:L.griffi thii—Yadong County,Shigatse City(27°32′29.21″N,88°59′40.94″ E;3080 m a.s.l.);L.speciose—Bomi County,Nyingchi City (29°50′20.82" N,95°50′31.74″ E;2650 m a.s.l.);L.himalaica—Jilong County,Shigatse City (28°33′8" N,85°19′6″ E,3410 m a.s.l.);L.kongboensis—Milin County,Nyingchi City (29°1′9.1" N,93°53′13.7″ E;3210 m a.s.l.);andL.potaniniivar.australis—Mankang County,Qamdo City (29°19′35.09" N,98°40′27.94″ E;3958 m a.s.l.).The voucher specimens are held in the herbarium of the Institute of Plateau Ecology of Tibet Agriculture &Animal Husbandry University,Nyingchi.Total DNA was extracted from leaves using a modified CTAB method (Healey et al.2014).The libraries(insert size=500) were constructed and sequenced on highthroughput Illumina HiSeq X Ten platform at Novogene Inc.(Beijing,China).

Genome assembly and annotation

The resultant paired-end 150-bp raw reads were trimmed,f iltered with an average quality less than Q5 or N content greater than five using Trimmomatic v0.36 (Bolger et al.2014),and imported and assembled into contigs using the de novo assembly tool of the CLC Genomics Workbench v10 with default settings (QIAGEN;Aarhus,Denmark).Reads coverage was checked using Bowtie 2 (Langmead and Salzberg 2012),and contigs with a coverage of less than 20 were filtered out.The top 10 contigs in length were blast against NCBI Nucleotide database to screen cp genomic contigs with a homology over 97% with deposited cp genome of other species,among which the longest one showed its best hit of cp genome ofLarix potaniniivar.chinensisunder the accession number KX880508.Thus,this genome was used as an initial reference for assembling complete cp genomes of the fiveLarixspecies using MITObim v1.8 with default settings (Hahn et al.2013).The initial assembly products were corrected using the contigs from CLC Genomics Workbench v10,and the low-quality regions and boundaries were verif ied by PCR and Sanger sequencing.The genome annotation was completed with DOGMA (Wyman et al.2004)and GENEIOUS R10 (Kearse et al.2012),combined with manual adjustments.Construction of each genome map was completed using the web tool OGDRAW v1.3.1 (Greiner et al.2019).

Comparative and phylogenetic analysis of the cp genomes

Simple sequence repeats (SSR) were analyzed using Sci-RoKo v3.4 (Kof ler et al.2007) under the perfect MISAmode with minimum repeat times of mono-,di-,tri-,tetra-,penta-and hexanucleotide set to 10,7,5,4,4,and 4,respectively.Genomes alignment was done using the LAGAN module of the web tool mVISTA (http://genom e.lbl.gov/vista/mvist a/submi t.shtml) with the annotation of cp genome ofL.griffi thiias reference.DNA polymorphism analysis was performed using DnaSP v5 with default settings(Librado and Rozas 2009),and the results were displayed with 500 bp for window length and step size.Expansion and contraction of the Inverted Repeats (IRs) were manually analyzed.Cp genomes of species ofLarix,including the five species in this study,as well as species ofCedrus,Pinus,andPiceawere selected for phylogenetic analysis within the Pinaceae family.A phylogenetic tree of these cp

genomes was constructed using MEGA v6 (Tamura et al.2013) based on the Maximum-Likelihood (ML) analysis of the concatenated coding sequence of protein-coding genes with a bootstrap value of 1000.

Results

Genome sequencing assembly and characterization

The complete cp genomes of the 5Larixspecies were successfully retrieved from Illumina reads sequenced from the whole-genome DNA.The quantities of isolated cp genomic reads were 199,897 forL.griffi thii,84,436 forL.speciosa,293,314 forL.himalaica,52,468 forL.kongboensis,and 73,173 forL.potaniniivar.australis.The 5 cp genomes range in size from 121,334 bp (L.kongboensis) to 121,971 bp (L.himalaica) (Fig.1).Reads coverage of all cp genomes are uniform and the mean coverage is 80.1 × forL.griffi thii,104.7 × forL.speciosa,360.9 × forL.himalaica,150.1 × forL.kongboensis,and 90.3 × forL.potaniniivar.australis.They each encode the same panel of 72 protein-coding genes,34 tRNAs,and 4 rRNAs,and are also conservative in gene type and arrangement.The annotated sequences have been deposited into Genbank of the National Center for Biotechnology Information in the United States under accession numbers MN822882-MN822886.The corresponding raw data has also been deposited into the SRA database of Genbank under PRJNA604660.

Fig.1 Physical maps of the chloroplast genomes of the 5 Larix species

Simple sequence repeats analysis

SSR scanning revealed low occurrence frequencies of SSRs in the 5 cp genomes,and all are mononucleotide (T)and dinucleotide (AT) repeats (Table 1).L.griffi thii,L.himalaica,andL.speciosacontained three SSRs in similar positions,whereasL.kongboensisandL.potaniniivar.australiscontained four SSRs.Three of the detected SSRs are polymorphic across the five cp genomes.

Table 1 Description of SSRs detected in the chloroplast genomes of the 5 Larix species

Alignment of the genomes

Alignment analysis revealed high similarities among the cp genomes of the 5Larixspecies.Obvious divergences were detected in the protein-coding sequence ofycf1and in the intergenic regions ofccsA-rpl32,rpl32-rps12,andpsbI-atpE.Almost all other divergences fell within the intergenic regions (Fig.2).DNA polymorphism analysis showed 241 sequence polymorphisms detected among the five cp genomes,and nucleotide variability (Pi) ranged from 0 to 0.148 (Table S1),indicating that DNA polymorphism among the 5 cp genomes is not robust.Nucleotide variabilities of 15 were above 0.003,with six located within theycf1proteincoding region,one within thepsbNcoding sequence and the remaining eight within the intergenic regions (Fig.3).Given that the abundance of repetitive fragments in theycf1protein-coding sequence often prevents the design of suitable universal primers,the polymorphic SSRs and nucleotides located within non-ycf1protein-coding regions were manually screened as candidate cp DNA markers for future studies (Table 2).

Table 2 Candidate polymorphic DNA markers among chloroplast genomes of five Larix species from the QTP

Fig.2 Alignment of chloroplast genomes of the 5 Larix species.The Y-scale represents the percent identity between 50 and 100%.Arrows represent gene sequences and transcriptional directions,blue represents exons and red represents conserved non-coding sequences (CNS)

Fig.3 Sliding window of nucleotide variability among the chloroplast genomes of 5 Larix species from the QTP.X-axis represents midpoint position of a window,Y-scale the nucleotide variability of each window.Peaks with a π value above 0.003 were annotated.Position is based on the alignment file

Inverted repeats (IRs) contraction and expansion

The contraction and expansion of IRs often leads to changes in the borders between an IR and a single copy (SC).However,in this study,the borders between IRs and SCs showed clear conservativeness,with the only difference at the IRBSSC border (Fig.4).Only the lengths fromtrnF(located within SSC) to the SSC-IRB borders were divergent among the 5 cp genomes,ranging from 1792 (L.griffithii) to 1806 bp (L.speciosa).There was no other change in the borders,implying that the IRs among the five cp genomes have experienced no drastic contraction or expansion (Fig.4).

Fig.4 Contraction and expansion of IRs in the chloroplast genomes of the 5 Larix species;the red box shows the divergent regions among the chloroplast genomes

Phylogenetic analysis

To examine the phylogenetic relationship among the cp genomes of these 5Larixspecies and several related species,a phylogenetic tree was constructed based on coding sequences from cp genomes of 16 species of Pinaceae with a bootstrap value of 1000.The cp genomes ofLarixwere clustered into 2 major clades,with one clade comprising onlyL.gmeliniivar.japonicaand the other clade the remaining 12 taxa (Fig.5).Furthermore,this clade could be further divided into 2 minor clades with the 5Larixspecies in this study clustered into the same minor clade,indicating a relatively closer interrelationship.It should be noted that one variety ofL.potaninii(L.potaniniivar.australis) is relatively distant from two other varieties (L.potaniniivar.chinensisandL.potaniniivar.macrocarpa) as reported previously (Han et al.2017;Qiu et al.2018) (Fig.5).

Fig.5 Phylogenetic trees of 16 taxa within the family Pinaceae based on the protein-coding sequences of the cp genomes;the 5 Larix species in this study are highlighted with a light blue background

Discussion

The continuous isolation in the mountainous QTP results in distinct microclimates for plant growth,providing the conditions for species differentiation.The cp genomes of the 5Larixspecies in the QTP are relatively conservative,especially in genome length,gene content and arrangement.The sizes of the 5 cp genomes (121,334-121,971 bp) are similar to those of otherLarixtaxa (−120 kb) (Lin et al.2010;Wu et al.2011;Han et al.2017;Kim et al.2018;Qiu et al.2018;Zimmermann et al.2019),indicating that all cp genomes ofLarixspecies are relatively conserved.However,there are still obvious sequence variations among the cp genomes of the 5Larixspecies in this study.The top 15 variations are almost in the intergenic region except for the ones inycf1protein-coding regions,andycf1variation was obvious among the 5 cp genomes.This gene,so far with unknown function,is often large and a variation hotspot across many plant cp genomes.For example,approximately 8% of all variable positions are located within theycf1region among 19 cp genomes ofLarixspecies from the southern lowlands of the Taymyr Peninsula and from other Pinaceae plants (Parks et al.2009;Lin et al.2010;Wu et al.2011;Kim et al.2018;Bondar et al.2019;Zimmermann et al.2019).Theycf1sequence is also highly variable in the cp genomes of other angiosperms (Gao et al.2018).One of the very important reasons for theycf1variation is that it contains many repetitive fragments,and the repetition times are highly variable,even among the cp genomes of close species (Ishizuka et al.2017).ycf1contains numerous repetitive fragments that may affect the genetic diversity of the species,and was recommended as plastid DNA markers for inferring intergeneric or interspecific phylogeny (Dong et al.2015).However,our point of view is in line with that of Zimmermann et al.(2019) who proposed that long repetitive fragments often block the design of the ideal primers or cause nested peaks in sanger sequencing results of amplicons.Thus,ycf1sequence variation is not of practical value and is not recommended as a marker for studies of evolution.This is also the main reason whyycf1was not considered as a candidate DNA marker in this study.

Polymorphic cpSSRs are some of the most common molecular markers which have been widely identified and used in population genetic assays of woody plants (Ni et al.2018;Wang et al.2018;Lee et al.2019;Liu et al.2019;Worth et al.2019;Yang et al.2020).However,cpSSRs were in low frequency and occurred only in the form of mononucleotide (T) and dinucleotide (AT) in this study.It consists with previous reports that most SSRs in cp genomes are A/T enriched (Gao et al.2018,2019;Zhou et al.2019).Given the relatively low SSR polymorphism and nucleotide variability of the five cp genomes,we suggest the combined use of these 2 types of polymorphisms,and a selected 12 candidate cp genome DNA markers would be helpful for future relevant population genetic studies.

Conclusions

Complete chloroplast genomes of 5Larixspecies from the Qing-Tibet Plateau were assembled and compared in this study.The 5 cp genomes are relatively conservative in genome size,gene content and arrangement,and border variation.Phylogenetic analysis revealed a close relationship among these 5 species.In addition,to improve the accuracy of further population genetic studies,12 candidate cp DNA markers were characterized based on polymorphic sites and SSRs among the 5 cp genomes.These results will provide genetic information for the conservation,evolution and ecology of the 5Larixspecies and their relatives.