Complete Chloroplast Genomes and Comparative Analysis of Sequences Evolution among Seven Aristolochia (Aristolochiaceae) Medicinal Species

Aristolochiaceae, comprising about 600 species, is a unique plant family containing aristolochic acids (AAs). In this study, we sequenced seven species of Aristolochia, and retrieved eleven chloroplast (cp) genomes published for comparative genomics analysis and phylogenetic constructions. The results show that the cp genomes had a typical quadripartite structure with conserved genome arrangement and moderate divergence. The cp genomes range from 159,308 bp to 160,520 bp in length and have a similar GC content of 38.5%–38.9%. A total number of 113 genes were identified, including 79 protein-coding genes, 30 tRNAs and four rRNAs. Although genomic structure and size were highly conserved, the IR-SC boundary regions were variable between these seven cp genomes. The trnH-GUG genes, are one of major differences between the plastomes of the two subgenera Siphisia and Aristolochia. We analyzed the features of nucleotide substitutions, distribution of repeat sequences and simple sequences repeats (SSRs), positive selections in the cp genomes, and identified 16 hotspot regions for genomes divergence that could be utilized as potential markers for phylogeny reconstruction. Phylogenetic relationships of the family Aristolochiaceae inferred from the 18 cp genome sequences were consistent and robust, using maximum parsimony (MP), maximum likelihood (ML), and Bayesian analysis (BI) methods.


Introduction
Aristolochia sensu lato, comprising about 500 species, is the largest genus of Aristolochiaceae, with a broad distribution range from tropical to subtropical, extending to temperate regions [1,2]. Several species of Aristolochia, such as Aristolochia moupinensis, Aristolochia tagala, and Aristolochia mollissima, have been reported as traditional Chinese medicines [3,4]. Aristolochiaceae is a unique plant family containing aristolochic acids (AAs), and their derivatives are widely implicated in liver cancers [5,6]. However, current studies have demonstrated that AAs are of nephrotoxicity, carcinogenicity, and mutagenicity [7,8]. The sale and use of AA-containing herbal preparations have been restricted in many countries [9].
The monophyly of Aristolochiaceae was well supported in most analysis, and was divided into two subfamilies, Asaroideae and Aristolochioideae [10,11]. The studies recognized two genera Saruma and Asarum in Asaroideae [10][11][12][13]. The genus Aristolochia of subfamily Aristolochioideae Int. J. Mol. Sci. 2019, 20, x FOR PEER REVIEW 6 of 23 C Figure 1. Gene maps of the complete cp genome of seven species of Aristolochia. Gene map of cp genome of (A) Aristolochia manshuriensis; (B) Aristolochia kaempferi, Aristolochia macrophylla, Aristolochia mollissima and Aristolochia kunmingensis; (C) Aristolochia tagala and Aristolochia tubiflora. Genes on the inside of the circle are transcribed clockwise, while those outside are transcribed counter clockwise. The darker gray in the inner circle corresponds to GC content, whereas the lighter gray corresponds to AT content.

IR Contraction and Expansion
The IR regions are expanded in five species of subgenus Siphisia compare with other two species (A. tagala and A. tubiflora) of subgenus Aristolochia, indicated by different duplication genes in the IR regions, where eight or seven tRNA genes were duplicated, respectively ( Figure 1, Table 2). The size of the IR region of subgenus Siphisia varies from 25,664 bp to 25,700 bp, and is 25,242 bp and 25,431 bp in the two plastomes of subgenus Aristolochia (Table 1).

IR Contraction and Expansion
The IR regions are expanded in five species of subgenus Siphisia compare with other two species (A. tagala and A. tubiflora) of subgenus Aristolochia, indicated by different duplication genes in the IR regions, where eight or seven tRNA genes were duplicated, respectively ( Figure 1, Table 2). The size of the IR region of subgenus Siphisia varies from 25,664 bp to 25,700 bp, and is 25,242 bp and 25,431 bp in the two plastomes of subgenus Aristolochia ( Table 1).
Fluctuation of IR-SC borders, together with the adjacent genes, were examined among seven Aristolochia species and six plastomes retrieved from GenBank (including Aristolochia contorta: NC_036152.1, Aristolochia debilis: NC_036153.1, Asarum canadense: MG544845-MG544851, Saruma henryi: MG520100, Piper auritum: NC_034697.1, and Drimys granadensis: NC_008456.1) (Figure 3). The LSC-IRb border, was located within the genic spacer of rps19-trnH for A. kaempferi, A. macrophylla, and A. mollissima (Type I), within the rps19 gene for A. kunmingensis and A. moupinensis (Type II), while in the rps19-rpl2 spacer for A. tagala and A. tubiflora (Type III). There were two types of SSC-IRa border among 13 detected species. In the three plastomes (A. moupinensis, A. tubiflora, and A. tagala), which ycf1 gene was fully located in the SSC region, and 25-43 bp apart from the SSC-IRa border. The SSC-IRa border was situated in the coding region ycf1 gene in the other 10 sequenced species, which spanned into the IRa region. Among the 10 detected species, the pseudogene ycf1 in the IRb region with the same length as far as the IRa expanded into ycf1 gene, and the length ranged from 153 bp to 2271 bp. The ndhF gene was entirely located in the SSC region in 10 species of Aristolochiaceae, but varied in distance (11-80 bp) from the IRb-SSC border. The LSC-IRa border in the species of subgenus Aristolochia was situated in the trnH gene with 10 bp into the IRa region (Type III), while the border was located in the trnH-psbA spacer in subgenus Siphisia species (Type I and II) ( Figure 3).

Codon Usage
All the protein-coding genes were composed of 26,398 codons in the cp genomes of the seven species of Aristolochia. The codon usages of protein-coding genes in the cp genomes are summarized in Figure 4 and Table S2. Among these codons, the most common amino acid in the protein-coding genes is leucine, which appears 2775 times in A. kaempferi and A. mollissima. The relative synonymous codon usage (RSCU) value analysis showed that almost all amino acids have more than one synonymous codon, except methionine and tryptophan. Nearly all of the proteincoding genes of Aristolochia species had the standard ATG start codon (RSCU = 1). About half of codons have RSCU > 1, and most of those (29/31, 93.5%) end with base A or T. About half of the codons have RSCU < 1, and most of those (28/31, 90.3%) end with base C or G.

Codon Usage
All the protein-coding genes were composed of 26,398 codons in the cp genomes of the seven species of Aristolochia. The codon usages of protein-coding genes in the cp genomes are summarized in Figure 4 and Table S2. Among these codons, the most common amino acid in the protein-coding genes is leucine, which appears 2775 times in A. kaempferi and A. mollissima. The relative synonymous codon usage (RSCU) value analysis showed that almost all amino acids have more than one synonymous codon, except methionine and tryptophan. Nearly all of the protein-coding genes of Aristolochia species had the standard ATG start codon (RSCU = 1). About half of codons have RSCU > 1, and most of those (29/31, 93.5%) end with base A or T. About half of the codons have RSCU < 1, and most of those (28/31, 90.3%) end with base C or G.

Positive Selection Analysis
We compared the ratio of non-synonymous (dN) and synonymous (dS) substitution for 79 protein-coding genes among seven species, including A. kunmingensis, A. kaempferi, A. tagala, A. debilis, As. canadense, S. henryi, and P. auritum within Piperales. The statistical neutrality test showed that five genes in the seven cp genomes are under significant positive selection, and these genes are involved in the synthesis of ribosomal small and large subunit protein (rps12, rps18, and rpl20) or unknown function (ycf1 and ycf2) ( Table 4). Likelihood ratio tests (M1a vs. M2a, M7 vs. M8) supported the presence of positively selected codon sites (p < 0.05) (Table S3). According to the M2a and M8 models, the rpl20 harbored three or four sites under positive selection. The gene ycf1 harbored one or three sites under positive selection based on two models, respectively. In addition, we identified rps12 gene with one positive selection site.

Repeat Structure and Simple Sequence Repeats Analyses
Repeats in ten cp genomes were analyzed using REPuter, including seven species of Aristolochia, S. henryi, P. auritum, and D. granadensis ( Figure 5, Table S4). The results showed that A. macrophylla had the greatest number of repetitive elements in cp genome, comprised of 25 forward, 26 palindromic, 21 reverse, and eight complement repeats. The size of the most repeats were 30-39 bp, and the repeats with the length > 49 bp only occurred in cp genomes of S. henryi and P. auritum. The

Positive Selection Analysis
We compared the ratio of non-synonymous (dN) and synonymous (dS) substitution for 79 protein-coding genes among seven species, including A. kunmingensis, A. kaempferi, A. tagala, A. debilis, As. canadense, S. henryi, and P. auritum within Piperales. The statistical neutrality test showed that five genes in the seven cp genomes are under significant positive selection, and these genes are involved in the synthesis of ribosomal small and large subunit protein (rps12, rps18, and rpl20) or unknown function (ycf1 and ycf2) ( Table 4). Likelihood ratio tests (M1a vs. M2a, M7 vs. M8) supported the presence of positively selected codon sites (p < 0.05) (Table S3). According to the M2a and M8 models, the rpl20 harbored three or four sites under positive selection. The gene ycf1 harbored one or three sites under positive selection based on two models, respectively. In addition, we identified rps12 gene with one positive selection site.

Repeat Structure and Simple Sequence Repeats Analyses
Repeats in ten cp genomes were analyzed using REPuter, including seven species of Aristolochia, S. henryi, P. auritum, and D. granadensis ( Figure 5, Table S4). The results showed that A. macrophylla had the greatest number of repetitive elements in cp genome, comprised of 25 forward, 26 palindromic, 21 reverse, and eight complement repeats. The size of the most repeats were 30-39 bp, and the repeats with the length > 49 bp only occurred in cp genomes of S. henryi and P. auritum. The longest repeats, with a length of 1591 bp, was detected in S. henryi. The total numbers of SSRs were also identified in the cp genomes of the ten species ( Figure 6 and Table S5). Mononucleotide repeats were the largest in a number of these SSRs, with 88% and 85% found in A. tubiflora and A. tagala, respectively. A/T repeats were the most common of mononucleotides, while AT/TA repeats are the majority of dinucleotide repeat sequences (96.3%-100%). The trinucleotide in the five species of subgenus Siphisia were only comprised of AAT/ATT repeats (100%), while A. tubiflora and A. tagala of subgenus Aristolochia also comprised AAC/GTT and AAG/CTT repeats. longest repeats, with a length of 1591 bp, was detected in S. henryi. The total numbers of SSRs were also identified in the cp genomes of the ten species ( Figure 6 and Table S5). Mononucleotide repeats were the largest in a number of these SSRs, with 88% and 85% found in A. tubiflora and A. tagala, respectively. A/T repeats were the most common of mononucleotides, while AT/TA repeats are the majority of dinucleotide repeat sequences (96.3%-100%). The trinucleotide in the five species of subgenus Siphisia were only comprised of AAT/ATT repeats (100%), while A. tubiflora and A. tagala of subgenus Aristolochia also comprised AAC/GTT and AAG/CTT repeats. Figure 5. Repeat sequences in ten cp genomes. REPuter was used to identify repeat sequences with length ≥ 30 bp and sequence identity ≥ 90% in the cp genomes. F, P, R, and C indicate the repeat types F (forward), P (palindrome), R (reverse), and C (complement), respectively. Repeats with different lengths are indicated in different colors. Figure 6. Frequency of simple sequence repeats (SSRs) in the ten cp genomes.

Comparative Genomic Divergence and Hotspots Regions
The SC and IR regions of cp genomes of the seven species (including A. moupinensis, A. kunmingensis, A. tagala, A. contorta, S. henryi, As. canadense, and P. auritum) were compared using the mVISTA program to detect hyper-variable regions (Figure 7). The alignment revealed high sequence conservatism across the cp genomes of A. moupinensis and A. kunmingensis of subgenus Siphisia. The comparison among seven cp genomes showed that the IR region was more conserved than the SC regions. The most divergent regions located in the intergenic spacers, and the most divergent coding regions were ndhF and ycf1. longest repeats, with a length of 1591 bp, was detected in S. henryi. The total numbers of SSRs were also identified in the cp genomes of the ten species ( Figure 6 and Table S5). Mononucleotide repeats were the largest in a number of these SSRs, with 88% and 85% found in A. tubiflora and A. tagala, respectively. A/T repeats were the most common of mononucleotides, while AT/TA repeats are the majority of dinucleotide repeat sequences (96.3%-100%). The trinucleotide in the five species of subgenus Siphisia were only comprised of AAT/ATT repeats (100%), while A. tubiflora and A. tagala of subgenus Aristolochia also comprised AAC/GTT and AAG/CTT repeats. Figure 5. Repeat sequences in ten cp genomes. REPuter was used to identify repeat sequences with length ≥ 30 bp and sequence identity ≥ 90% in the cp genomes. F, P, R, and C indicate the repeat types F (forward), P (palindrome), R (reverse), and C (complement), respectively. Repeats with different lengths are indicated in different colors. Figure 6. Frequency of simple sequence repeats (SSRs) in the ten cp genomes.

Comparative Genomic Divergence and Hotspots Regions
The SC and IR regions of cp genomes of the seven species (including A. moupinensis, A. kunmingensis, A. tagala, A. contorta, S. henryi, As. canadense, and P. auritum) were compared using the mVISTA program to detect hyper-variable regions (Figure 7). The alignment revealed high sequence conservatism across the cp genomes of A. moupinensis and A. kunmingensis of subgenus Siphisia. The comparison among seven cp genomes showed that the IR region was more conserved than the SC regions. The most divergent regions located in the intergenic spacers, and the most divergent coding regions were ndhF and ycf1.

Comparative Genomic Divergence and Hotspots Regions
The SC and IR regions of cp genomes of the seven species (including A. moupinensis, A. kunmingensis, A. tagala, A. contorta, S. henryi, As. canadense, and P. auritum) were compared using the mVISTA program to detect hyper-variable regions (Figure 7). The alignment revealed high sequence conservatism across the cp genomes of A. moupinensis and A. kunmingensis of subgenus Siphisia. The comparison among seven cp genomes showed that the IR region was more conserved than the SC regions. The most divergent regions located in the intergenic spacers, and the most divergent coding regions were ndhF and ycf1. Figure 7. Sequence identity plot compared seven cp genomes with A. moupinensis as a reference by using mVISTA. Grey arrows and thick black lines above the alignment indicate genes with their orientation and the position of the IRs, respectively. A cut-off of 70% identity was used for the plots, and the Y-scale represents the percent identity from 50% to 100%.
Comparative analysis among our seven sequenced species within Aristolochia was conducted of the entire cp genomes, LSC, SSC, IR, and CDS regions, respectively ( Table 5). The nucleotide diversity (Pi) value was also calculated to evaluate the sequence divergence among these cp genomes, and their values varied from 0 to 0.07746 (Figure 8). The analysis revealed that the SSC region, compared with other regions, exhibited the highest levels of divergence (Pi = 0.03114). These values of the LSC region, varied from 0.00175 to 0.07746, with the mean value of 0.02182. The IR region exhibited the lowest Pi values varying from 0 to 0.01056, with the mean of 0.00411, indicating that IR region was the most conserved one. Furthermore, we identified 16 hotspot regions (Pi > 0.04, the mean value = 0.05413) with the full length of 20,296 bp, including rps16-trnQ-psbK, psbI-trnS-trnG, atpH-atpI, psbM-trnD, rps4-trnT-trnL, trnF-ndhJ, ndhC-trnV, accD-psaI, petA-psbJ, rps18-rpl20, trnN-ndhF, rpl32-trnL-ccsA, and four regions of ycf1 coding gene ( Table 6). Ten of these (rps16-trnQ-psbK, psbI-trnS-trnG, atpH-atpI, psbM-trnD, rps4-trnT-trnL, trnF-ndhJ, ndhC-trnV, accD-psaI, petA-psbJ, and rps18-rpl20) are located in the LSC, and six (trnN-ndhF, rpl32-trnL-ccsA and ycf1) in the SSC region, which could be utilized as potential markers for the phylogeny reconstruction and species identification of this subgenus in further studies.  . Sequence identity plot compared seven cp genomes with A. moupinensis as a reference by using mVISTA. Grey arrows and thick black lines above the alignment indicate genes with their orientation and the position of the IRs, respectively. A cut-off of 70% identity was used for the plots, and the Y-scale represents the percent identity from 50% to 100%.

Phylogenetic Analyses
The phylogenetic relationships of Aristolochiaceae were constructed based on six datasets (entire cp genome sequences except a copy of IR, LSC, SSC, IR, and CDS regions and combining 16 hotspots) of 18 samples, using three methods of ML, MP, and BI, respectively ( Figure 9). The robust topologies were consistent for most clades of cp genomes, LSC, SSC, CDS, and hotspots datasets, with the high bootstrap values for most of the branches ( Figure 9A). From these six different datasets, the phylogenetic analysis showed that the genera Asarum and Saruma represented by seven species formed a clade with posterior probabilities (PP) = 1 based on BI, bootstrap values (%) (BS) =100 based on ML and BS =100 based on MP methods. However, the tree constructed using sequences of the IR region failed to resolve the phylogeny position of Asarum epigynum and As. canadense ( Figure 9B

IR Contraction and Expansion
Taken another two reported species (A. debilis and A. contorta) of subgenus Aristolochia into account, although genomic structure and size were highly conserved, the IR-SC boundary regions were variable between these nine cp genomes of Aristolochia (Figure 3). In general, contraction and expansion at the borders of IR regions are common evolutionary events and may cause IR size variation of plastomes [29,[34][35][36]. The length of the IR regions of five Siphisia species, varying in the range of 25,664-25,700 bp, was longer than those of the four species of subgenus Aristolochia, which varied from 25,175 bp to 25,459 bp (Table 1) [28]. We identified three types of the IR-SC junctions from the nine Aristolochia species, according to the organization of genes ( Figure 3). Within five detected species of subgenus Siphisia, its patterns were Type I and II, while the Type III only occurred in the four species of subgenus Aristolochia. Type I was found in A. mollissima, A. macrophylla, and A. kaempferi, and was characterized by trnH gene in IR region and LSC-IRb border located in rps19-trnH spacer. Type II was only found in A. moupinensis and A. kunmingensis and refers to LSC-IRb border within the rps19 gene. The trnH gene is intact and located upstream of rpl2 in IRb region for type I and II. Type III pattern was found in the four species of subgenus Aristolochia, characterized by LSC-IRb and SSC-IRa border in the rps19-rpl2 spacer and trnH gene, respectively. The trnH gene spanned the junction between IR-LSC regions in the four species of subgenus Aristolochia.
The shift of IR-LSC borders, caused by contraction and expansion of the gene trnH, is one of major differences between the plastomes of the subgenera Siphisia and Aristolochia. The whole gene duplication of trnH was detected in most monocots (e.g., Acorus, Phalaenopsi and Dioscorea), D. granadensis (Winteraceae) of magnoliids, and basal eudicots (Ranunculus japonica and Ranunculus macranthus) [34,[37][38][39][40][41]. Wang et al. (2008) conducted RT-PCR assays and deduced that the duplicated trnH genes in most of non-monocots and monocots were regulated by different expression levels of promoters, and had distinct fates [37]. Within the family Aristolochiaceae, the trnH gene was located in the LSC region of S. henryi, 128 bp away from the border of LSC-IR, and was also a single copy in the six cp genomes of Asarum, but not sure the positions of the gene [29]. Furthermore, the study proposed that the low-complexity trnH region and ultimately inversion of a portion of the LSC were due to an AAT repeat. For inversion of a large portion of the LSC region, there were genes rearranged in SC-IR borders of sequenced species of Asarum, the IR boundaries of cp genomes of Asarum were highly

IR Contraction and Expansion
Taken another two reported species (A. debilis and A. contorta) of subgenus Aristolochia into account, although genomic structure and size were highly conserved, the IR-SC boundary regions were variable between these nine cp genomes of Aristolochia (Figure 3). In general, contraction and expansion at the borders of IR regions are common evolutionary events and may cause IR size variation of plastomes [29,[34][35][36]. The length of the IR regions of five Siphisia species, varying in the range of 25,664-25,700 bp, was longer than those of the four species of subgenus Aristolochia, which varied from 25,175 bp to 25,459 bp (Table 1) [28]. We identified three types of the IR-SC junctions from the nine Aristolochia species, according to the organization of genes ( Figure 3). Within five detected species of subgenus Siphisia, its patterns were Type I and II, while the Type III only occurred in the four species of subgenus Aristolochia. Type I was found in A. mollissima, A. macrophylla, and A. kaempferi, and was characterized by trnH gene in IR region and LSC-IRb border located in rps19-trnH spacer. Type II was only found in A. moupinensis and A. kunmingensis and refers to LSC-IRb border within the rps19 gene. The trnH gene is intact and located upstream of rpl2 in IRb region for type I and II. Type III pattern was found in the four species of subgenus Aristolochia, characterized by LSC-IRb and SSC-IRa border in the rps19-rpl2 spacer and trnH gene, respectively. The trnH gene spanned the junction between IR-LSC regions in the four species of subgenus Aristolochia.
The shift of IR-LSC borders, caused by contraction and expansion of the gene trnH, is one of major differences between the plastomes of the subgenera Siphisia and Aristolochia. The whole gene duplication of trnH was detected in most monocots (e.g., Acorus, Phalaenopsi and Dioscorea), D. granadensis (Winteraceae) of magnoliids, and basal eudicots (Ranunculus japonica and Ranunculus macranthus) [34,[37][38][39][40][41]. Wang et al. (2008) conducted RT-PCR assays and deduced that the duplicated trnH genes in most of non-monocots and monocots were regulated by different expression levels of promoters, and had distinct fates [37]. Within the family Aristolochiaceae, the trnH gene was located in the LSC region of S. henryi, 128 bp away from the border of LSC-IR, and was also a single copy in the six cp genomes of Asarum, but not sure the positions of the gene [29]. Furthermore, the study proposed that the low-complexity trnH region and ultimately inversion of a portion of the LSC were due to an AAT repeat. For inversion of a large portion of the LSC region, there were genes rearranged in SC-IR borders of sequenced species of Asarum, the IR boundaries of cp genomes of Asarum were highly variable and experienced positional shifts at borders. Such as there was an entirety of the SSC of As. canadense and As. sieboldii has been incorporated into the IR, and the boundary of the LSC-IR was found within rpl2 or rpl14 gene [29]. Within the species of S. henryi, rps19 pseudogene existed in the IRa region, with the length of 183 bp. The trnH-rps19 gene cluster had been used to distinguish monocots from other angiosperm for the organization of gene flanking the IR-SC junction [37,39]. The events of contraction or expansion of the IR regions also can be used to distinguish the species within Aristolochiaceae.

Inferring the Phylogeny and Species Identification of Aristolochia
Chloroplast genomes provide abundant resources significant for evolutionary, taxonomic, and phylogenetic studies [42][43][44]. The whole cp genomes and protein-coding genes have been successfully used to resolve phylogenetic relationships at multiple taxonomic levels during the past decade [45,46]. Repeats can lead to changes in genomic structure, and can be investigated to population genetics of allied taxa [47][48][49][50]. Repeats in ten cp genomes revealed that the repeats had a great number, comprised of 38-80 repeats ( Figure 5 and Table S4), 66 and 138 repeats were respectively detected in A. debilis and A. contorta [28]. Given the variability of these repeats between lineages, they can be informative regions for developing genomic markers for phylogenetic analysis. SSRs, known as microsatellites, are tandemly repeated DNA sequences that consist of one-six nucleotide repeat units and are ubiquitous throughout the genomes [51]. A total number of 95-142 SSRs were identified in the seven cp genomes detected ( Figure 6 and Table S5). According to the analysis of high variable regions, the hotspot regions within seven cp genomes also provide sufficient information sites to reveal phylogeny structure among species of family Aristolochiaceae, especially for the spacer ycf1 and rpl20, with high nucleotide diversity and under positive selection ( Table 4). The ycf1 gene could be served as the barcode of land plants, and was also recognized as the most variable regions in plastid genome [50,52]. The gene rpl20 is an important part of protein synthesis, and is involved in translation [53]. This study will also provide a reference for phylogenomic studies of closely related lineages among Aristolochia and other genera.
Furthermore, we can design effective markers for clarifying the phylogenetic relationships of Aristolochia and elucidating the evolutionary history of species complex of Aristolochia at the population level, based on the analysis of SSR and SNP sites. Understanding genetic variation within and between populations plays an important role in improving genetic diversity and is essential for future adaptive changes, reproduction patterns, and its conservation [20,54,55]. The cpDNA and B-class gene PISTILLATA (PI) have been used to investigate taxonomy at the species complex, such as Aristolochia kaempferi group, and these studies revealed that its DNA barcoding and taxonomy are difficult to assess for multiple hybridization and introgression events in the group [56,57]. More genes under selection and neutral markers should be used to clarify those multiple diversification events. It will better to apply the full genome information or hyper-variation regions to elucidate the species diversity of Aristolochia.

Plant Material, DNA Extraction, and Sequencing
We selected seven species according to their potential medicinal uses, including A. kaempferi, A. kunmingensis, A. macrophylla, A. mollissima, and A. moupinensis from subgenus Siphisia and A. tagala and A. tubiflora of subgenus Aristolochia (Table 7). Genomic DNA was isolated from silica-gel dried leaf tissue or herbarium specimens using Plant Genomic DNA Kit (TIANGEN, Beijing, China). DNA integrity was examined by electrophoresis in 1% (w/v) agarose gel and their concentration was measured using a NanoDrop spectrophotometer 2000 (Thermo Scientific; Waltham, MA, USA). The DNA was used to construct PE libraries with insert sizes of 150 bp and sequenced according to the manufacturer's manual for the Illumina Hiseq X.
The cp genomes of the seven species was annotated using the online program Dual Organellar GenoMe Annotator (DOGMA) (University of Texas at Austin, Austin, TX, USA) [60], Annotation of Organellar Genomes (GeSeq) [61] and Chloroplast Genome Annotation, Visualization, Analysis, and GenBank Submission (CPGAVAS) (Institute of Medicinal Plant Development, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China) [62]. The tRNA genes were confirmed using tRNAscan-SE software (v2.0, University of California, Santa Cruz, CA, USA) [63]. Plastome annotations were manually corrected with the software Artemis [64]. The gene map was drawn using the Organelle Genome DRAW (OGDRAW) [65,66] with default settings and checked manually. The complete cp genome sequences of the seven species were deposited in GenBank, accession numbers are MK503927-MG503933 (Table S6).

Genome Structure Analyses
The distribution of codon usage was investigated using the software CodonW (University of Texas, Houston, TX, USA) with the RSCU value [67]. GC content was analyzed using Molecular Evolutionary Genetics Analysis (MEGA v6.0, Tokyo Metropolitan University, Tokyo, Japan) [68]. REPuter program (https://bibiserv.cebitec.uni-bielefeld.de/reputer) (University of Bielefeld, Bielefeld, Germany) [69] was used to identify the size and location of repeat sequences, including forward, palindromic, reverse, and complement repeats in the seven cp genomes. For all repeat types, the minimal size was set as 30 bp and the two repeat copies had at least 90% similarity. Perl script MISA (https://webblast.ipk-gatersleben.de/misa/) [70] was used to detect microsatellites (mono-, di-, tri-, tetra-, penta-, hexanucleotide repeats) with the following thresholds (unit size, min repeats): ten repeat units for mononucleotide SSRs, five repeat units for dinucleotide SSRs, four repeat units for trinucleotide SSRs, and three repeat units each for tetra-, penta-, and hexanucleotide SSRs.

Positive Selection Analysis
To identify the genes under selection, we scanned the cp genomes of seven species within Piperales using codeml of the package PAML4 [71,72]. The software was used for calculating the non-synonymous (dN) and synonymous (dS) substitution rates, along with their ratios (ω = dN/dS). The analyses of selective pressures were conducted along the ML tree in Newick format (S7), which based on the whole CDS region was used to determine the phylogenetic relationships of these seven species. Each single-copy CDS sequences was aligned according to their amino acid sequence. We used the site-specific model with five site models (M0, M1a & M2a, M7 & M8) were employed to identify the signatures of adaptation across cp genomes. This model allowed the ω ratio to vary among sites, with a fixed ω ratio in all the branches. Comparing the site-specific model, M1a (nearly neutral) vs. M2a (positive selection) and M7 (β) vs. M8 (β & ω) were calculated in order to detect positive selection [73]. Likelihood ratio test (LRT) of the comparison (M1a vs. M2a and M7 vs. M8) was used respectively to evaluate of the selection strength and the p value of Chi square (χ 2 ) smaller than 0.05 is thought as significant. The Bayes Empirical Bayes (BEB) inference [74] was implemented in site models M2a and M8 to estimate the posterior probabilities and positive selection pressures of the selected genes.

Genome Comparison and Nucleotide Variation Analysis
The whole-genome (minus a copy of IR region) alignment for the cp genomes of the seven species including our A. moupinensis, A. kunmingensis, A. tubiflora and four reported species (A. contorta, As. canadense, S. henryi and P. cenocladum) of Piperales, was performed and plotted by the mVISTA program (http://genome.lbl.gov/vista/mvista/submit.shtml) in Shuffle-LAGAN model [75,76], and with A. moupinensis as the reference. The seven cp genomes of Aristolochia were first aligned using MAFFT v7 [77] and then manually adjusted using BioEdit v7.0.9 [78]. Variable sites and nucleotide variability across complete cp genomes, LSC, IR, SSC, and CDS regions of seven species were calculated using DnaSP v5 [79]. Furthermore, for the seven cp genomes minus a copy IR region, a sliding window analysis was conducted to evaluate the nucleotide variability using DnaSP software. The step size was set to 200 base pairs, and the window length was set to 600 base pairs.

Phylogenetic Analyses
To estimate phylogenetic relationships within the Aristolochiaceae, plastomes of 18 taxa were compared, including nine samples from Aristolochia, six and one cp genomes from Asarum and Saruma, respectively (Table S5). A total of 11 cp genomes were downloaded from the NCBI database. In the phylogenetic analyses, P. auritum and P. cenocladum of Piper were used as outgroup. Phylogenetic trees were constructed by MP, ML and BI methods using the cp genomes, LSC, SSC, IR, CDS and hotspots regions. The sequences of the involved regions were aligned using MAFFT v7. MP analysis was performed with PAUP*4.0b10 [80], using a heuristic search performed 1000 replications and tree bisection-reconnection (TBR) branch swapping. BI was conducted using the program MrBayes v3.2 [81] with the GTR+I+G model at the CIPRES Science Gateway website (http://www.phylo.org/) [82]. The Markov Chain Monte Carlo (MCMC) analysis was run for 2,000,000 generations, sampling every 1000 generations. The posterior probabilities (PP) of the phylogeny and its branches were determined from the combined set of trees, discarding the first 25% trees of each run as burn-in, as determined by Tracer v1.7 [83]. Maximum likelihood phylogenies were constructed by a fast and effective stochastic algorithm using IQ-TREE v1.6.2 [84] with the Best-fit model by ModelFinder [85] according to Bayesian information criterion (BIC) and the robustness of the topology was estimated using 2000 bootstrap replicates. Figtree v1.4 (http://tree.bio.ed.ac.uk/software/figtree/) [86] was used to visualize and annotate trees.

Conclusions
The complete cp genomes of A. kaempferi, A. kunmingensis, A. macrophylla, A. mollissima, and A. moupinensis of the subgenus Siphisia, and A. tagala and A. tubiflora of the subgenus Aristolochia were reported in this study. The cp genomes length and gene content of the genus Aristolochia were comparatively conserved. Although genomic structure and size were highly conserved, the IR-SC boundary regions were variable between these nine cp genomes of Aristolochia. The whole duplicated trnH gene within five species of Siphisia is one of major differences between the plastomes of the subgenera Siphisia and Aristolochia. We also identified SSR sites, five positive selection sites and 16 variable regions, which provide a reference for developing tools to further study Aristolochia species. Furthermore, the phylogenetic constructions with six datasets of 18 cp genomes illustrated robust and consistent relationships with high supports.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/xxx/s1. Author Contributions: X.L. performed the experiments, analyzed the data, and wrote the manuscript; Y.Z. assembled sequences and revised the manuscript; X.Z. and S.L. collected, identified plant materials and gave suggestions to the manuscript; J.M. revised the manuscript. All authors have read and approved the final manuscript.
Funding: This work was supported by the National Natural Science Foundation of China (No. 31370225).

Acknowledgments:
The authors give special thanks to Shuwan Li, Yuan Wang, Zhanghua Wang, and Zhixiang Hua for collecting plant material. We acknowledged someone for their assistance with fieldwork, for data analysis, for giving comments on the manuscript paper. Our sincere thanks are also to the anonymous reviewers for their comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.