Evidence of a possible double hybrid origin for two Malagasy species of Piper L. (Piperaceae)

Two new species of genus Piper L. from Madagascar: Piper malgassicum Papini, Palchetti, M. Gori & Rota Nodari and Piper tsarasotrae Papini, Palchetti, M. Gori & Rota Nodari, were analyzed to investigate their phylogenetic position and evolutionary history. Both plastidial and nuclear markers were used for sequencing. The plastidial markers (ndhF and trnL intron) showed a close relationship between the two species with respect to the other species of Piper. Both species appeared phylogenetically related to the African P. guineense and the Malagasian/Mascarenhas endemic P. borbonense. The nuclear marker (G3pdh) amplification produced two separate sets of sequences: “long” sequences and “short” sequences, characterized by some long deletions. Analyzing together the nuclear sequences, we observed that the “long” sequence of P. tsarasotrae had a stricter relationship to the African accessions of P. guineense, while the accession of P. malgassicum was more strictly related to P. borbonense. On the contrary both “short” sequences of P. malgassicum and P. tsaratsotrae resulted phylogenetically related to Asian accessions and more distantly related to the formerly cited species. This unexpected result was tentatively explained with a more ancient hybridization event between an ancestor of P. malgassicum and P. tsarasotrae (and possibly P. borbonense) and an Asian species of Piper. The Asian contribution would have produced the ancestors carrying the “short” sequences. A more recent hybridization event would have led to the separation of P. malgassicum from P. tsarasotrae with an African pollen-derived genome contribution from P. guineense or, more probably, an ancestor thereof, to an ancestor of P. tsarasotrae. The chromosome numbers of P. tsarasotrae (2n = about 38) and P. malgassicum (2n = about 46), were more similar to the Asian species than to the American species. Unfortunately, no chromosome number of the African species P. guineense is currently available, to compare the chromosomal numbers.


Introduction
Genus Piper L. (Piperaceae) is one of the largest genera of Angiosperms, with more than 2000 species [1] and were considered belonging to a basal group of angiosperms, the so called "paleoherbs" [2].
Piper is a pantropical genus developing highly variable growth forms [3], with the highest biodiversity in the American continent with a number of species ranging from 500 [4] [5], to 1100 [6], later increased to more than 1800 [7], many of them with a small distribution area [1].
The separation of species is often tricky, due to the small dimension of the floral parts and hence the number of synonyms may be high [8], while other species tend to get naturalized [9]. While only two species are known as native of the African continent, P. guineense and P. capense, more species are known of Madagascar, even if some of them are known only for a single herbarium sample. The currently recognized species in Madagascar are P. heimii C. DC, P. pachyphyllum Baker and possibly P. borbonense (Miq.) C. DC., described for the Bourbon island, nowadays La Reunion [10], belonging to Mascarenhas Islands. However, its presence in Madagascar was affirmed by De Candolle [11] [12]. The fact that P. borbonense is cultivated makes more complex to understand its real distribution area [13].
Piper malgassicum Papini, Palchetti, M. Gori, Rota Nodari and Piper tsarasotrae Papini, Palchetti, M. Gori, Rota Nodari, were recently described as new Malagasy species [14] and are of economic interest, since their dried fruits are often mixed with P. borbonense to produce the typical Malagasy spice called in local language "voatsiperifery" pepper.

Results
Amplification of two plastid fragments named ndhF and trnL intron was carried on and the amplicons correctly sequenced producing reads of 1860 bp and 920 bp, respectively. Cloning of the amplicon of the nuclear gene G3pdh of P. malgassicum and P. tsarasotrae allowed to isolate two haplotypes, which were named "long" (1060bp for P. malgassicum and 1045bp for P. tsarasotrae) and "short" 965 bp (for both species) after their size. We used a total of 71 sequences, considering separately the short and long sequences of P. malgassicum and P. tsarasotrae for G3pdh and the plastid sequences matrix. The total alignment of the G3pdh region was 1127 nucleotides long including gaps. The final parts of the sequences were very variable and hence the alignment was ambiguous. For this reason, we excluded the characters from position 957 to 1127. The rest of the alignment was used for indels (gap) coding (with the software gapcoder), resulting in further 99 characters that were inserted after the nucleotide sequences. The plastid genes ndhF and trnL were inserted one after the other in the sequence, producing an aligned matrix of 2016 characters. The coding of indels resulted in further 115 characters. RAxML applied on the nuclear G3pdh matrix (indels coding excluded) produced a maximum likelihood tree with bootstrap support obtained with 1000 replicates (Figure 1). Figure 1. Maximum likelihood tree produced by RAXML with nuclear sequences. The supports above or below the branches are, respectively, the bootstrap resampling support with maximum likelihood criterion produced by RAXML, and the bayesian support calculated including the information derived from indels. In case the bayesian support is lower than 50, it is not indicated on the figure.
The support on branches corresponds to maximum likelihood bootstrap support (left) and Bayesian support with gaps (on the right). The same method was using for the plastid matrix ( Figure  2).  Comparing the two maximum likelihood trees, the one based on nuclear DNA data (G3pdh sequences) and that obtained with plastid markers, we could observe that in the first case P. malgassicum, clustered together and as sister group of P. borbonense (Figure 3), another species from the La Reunion Island, which lies relatively close to Madagascar. This relationship is corroborated by 100% maximum likelihood bootstrap (MLS) and bayesian (BS) support. The other Malagasy species, P. tsarasotrae, typical of arid forest, was more strictly related to the entries of the African species P. guineense, with 100% MLS and 100% BS. All these species formed a well characterized clade with 89% MLS and 100% BS and their closest species appeared to be Asian species P. caninum, (Figure 3). The BS without considering gaps coding gave the same support in this clade. The "short" sequences of G3pdh of both P. tsarasotrae and P. malgassicum clustered together within a group of Asian species, mainly originating from Malaysia and Australia with 98% MLS and 100% BS (Figure 1).
The (phylogenetic) story told by the data obtained from chloroplast genome sequences was quite different: the Malagasian species P. tsarasotrae and P. malgassicum clustered together with the phytogeopraphically close P. borbonense with 90% MLS and 100% BS, while the 5 accessions of the African P. guineense were in a more external condition with respect to the former group and separated in two groups, one from Cameroon (NW Africa) and one from Uganda/Kenya (Central-East Africa). All these species together formed a monophyletic group with 64% MLS and 96% BS (95% bayesian support in the analysis without gaps). Also in this case P. caninum, together with P. rothianum, was the outgroup to the African + Malagasy species (Figure 2) with 80% MLS and 98% BS (99% without gaps). Adding indels data to the matrix did not appear to increase the support value of nodes in the plastidial genes tree.
The counted chromosome numbers varied from 2n=46±2 in P. malgassicum ( Figure 3A) to 2n=36± 2 in P. tsarasotrae ( Figure 3B). The uncertainty in the counts, that should be taken only as preliminary result, derived from the small size of the chromosomes (many of them less than 1 μm of length), the low amount of metaphases in the root tips of the plants cultivated in Florence and the apparently small dimension of the mitotic spindle, leading to partial overlapping of many of the small chromosomes.

Discussion
The fact that the phylogenetic history based on the chloroplast markers told a different tale with respect to the tree produced with nuclear markers may be explained with a possible ancient hybridization/introgression event with pollen coming from an ancestor of the African P. guineense and reaching the ancestor of P. tsarasotrae, that would hence share some part of the nuclear genome with the African species. The only species of Piperaceae analyzed under the point of view of the type of plastid inheritance was a species of Peperomia, resulted only with maternal plastidial inheritance [14]. The presence of the short G3pdh nuclear sequences may be related to a still more ancient hybridization event involving the ancestor of the Malagasy species and some ancestor of Asian origin. Also in this case probably, with Asian pollen entering in contact with the ancestor´s stigma of the Malagasian species. As a matter of fact the closest relatives to the African species sensu lato (including the Malagasy and the Reunion species) are Asian, with the closest species (among those here sampled) apparently from Malaysia (Figure 3 and Figure 4). Apparently interspecific hybrids can be obtained in genus Piper also experimentally [21], while the hybrid origin of several Andean species was already proposed by Quijano-Abril et al. [1].
The presence of paralogs of G3pdh in angiosperms may represent a problem in several phylogenetic analysis [15] [16] [17]. However, here most of the indels were found in the introns of the gene and hence we are not able to assess the functionality of the short sequences.
The preliminary results about the chromosome numbers scored about 2n=46+-2 in P. malgassicum and 2n=36+-2 in P. tsarasotrae. The uncertainty in the counts was due to the small dimension of the chromosomes that were observed in most of the species of the genus, together with stickiness [18], the low amount of metaphases in the root tips of the plants cultivated in vitro and the apparently small dimension of the mitotic fuse, leading to partial overlapping of many of the small chromosomes. The mitotic spindle can reach dimensions up to 60 μm [19] [20], while in P. malgassicum and P. tsarasotrae it was about 15-20 μm (see Figure 5).
The chromosome numbers in genus Piper are very variable, ranging from 2n=26 to 2n=104, with some species apparently able to possess several possible chromosome numbers [18]. Most new world species show a karyotype of 2n=26, while in Asia tetraploids 2n=52 would prevail [18]. No data was available for African and Malagasy species up to the here presented results. However, the clear difference in karyotype between P. tsarasotrae and P. malgassicum, two species otherwise strictly related, may confirm a possible hybridization/introgression event with a species with a different chromosome number with respect to the ancestor of the Malagasian species.
A discordance between plastid and nuclear inheritance inferred through DNA sequencing has been often related to reticulate evolution and species of hybrid origin [22] [23] [24], even if the karyological data that may be a decisive evidence, has been rarely used in relationship to the DNA sequence evidence as, for instance in Selvi et al. [25]. The relationship between the Malagasian species and P. guineense with the asian species as members of Piper s. s. was already proposed by Jaramillo and Callejas [26] and Jaramillo et al. [27] as a result of a dispersal event, and our results do not disagree with this position. Apparently, in Africa and Madagascar the conditions leading to the wide diversification observed in South-american Piper [28] are lacking or less capable of influencing the speciation process.

Materials and Methods
A first round of sample collection within the internal area of Madagascar was conducted in 2016 and the samples have been submitted to analyses. The results have been reported in Palchetti et al. 2018 [13] but, in order to get a deeper knowledge about the genetic asset of the two species and to confirm the obtained results e second round of sample collection has been carried on in 2019 ( Figure  4). 4 new plants were collected in two different areas of the Ambositra region in Madagascar, the first 2 plants belonging to the P. malgassicum type were collected into the tropical rainy forest of Vohiday and the second 2 plants, belonging to the P. tsarasotrae type in the semi-dry area of Tsaratsotra village these plants were compared with the samples of P. tsarasotrae and P. malgassicum which have been used for the previous research [13]. Samples were conserved either in ethanol 96% either as herbarium sample by the ET (Tropical Herbarium of Florence, CSET, https://www.bio.unifi.it). Some seeds were also germinated in Florence for karyotyping. The DNA used for this work was extracted from tissue conserved in ethanol 96% [29] [30].
. Figure 4. Geographical localization of the sampling area for P. tsaratsotrae (yellow) and P. malgassicum (white). Area of the sampling campaign of 2018 (a and b) and 2019 (A and B) for P. tsaratsotrae and P. malgassicum respectively. DNA was extracted from 40 mg of the ethanol preserved leaves after drying under vacuum. The starting material was inserted in 2 ml tube, together with tungsten carbide beads, frozen in liquid nitrogen and finely ground in a tissue homogenizer (Tissue Lyser ®, Qiagen). DNA was extracted using Invisorb Spin Plant Mini kit (Stratec molecular®) according to the manifacture's guidelines.
Two new primer pairs were designed using the chloroplast genome sequence of P. kadsura (GenBank: KT223569.1) as template to cover the entire NADH dehydrogenase F (ndhF) plastid gene: ndhF-F3_forward 5'-AGGTTCTTATCGAGCCGCTT-3' and ndhF-F3_reverse 5'-GTAAGAAGAAATGCGCCCCC-3' and ndhF-F10_forward 5'-CTTCGCCGTATGTGGGCTTT-3' and ndhF-F10_reverse 5'-TCGACCAAAAGCAAGCAAGAG-3'. The amplicons have been directly and bi-directionally sequenced by using the corresponding primers for each amplified sequence. Since direct sequencing of G3pdh showed fragments of extra peaked sequencing data, we proceeded with cloning with InsTAclone PCR Cloning Kit (Thermo Scientific®) of the G3pdh amplification products. Several colonies for each cloned sample were amplified using T7 and SP6 primers whose sites are located at the boundaries of the cloning vector. PCR products were purified using the QIAquick PCR Purification Kit (Quiagen) and sent to the University of Florence internal sequencing service CIBIACI (www.cibiaci.unifi.it). Manual correction and assembly of the sequences was performed using the software Multalin [33] and MEGA7 [34]. Unexpectedly, two DNA sequences were obtained, after removing the cloning vector fragments, showing a different size: 965bp and 1058bp which were named "short" and "long" sequences respectively ( Figure 5). Figure 5. Alignment the long and short G3pdh sequences isolated from P .malgassicum and P. tsarasotrae using BioEdit software [35]. Shaded fragments represent the primers used for amplification.
At a first sight only the long sequences of G3pdh have been considered as the right ones [13] because, as observed by Smith et al. 2008 [9], no paralogs have been detected in a great deal of other Piper species and therefore other results have been discarded as PCR artifacts. In the present work a second thorough revision of the sequencing output has been carried on and showed the presence of overlapping peaks in all the samples and an additional round of analysis of the colonies confirmed the presence of the short sequences. To rule out any doubt, two additional specimen for each species has been collected and submitted to amplification and cloning in order to confirm the presence of the short sequence. As all the samples showed the same pattern, we decide to use also this "short" sequence to study the phylogeny of these two piper species, by comparing with all the accession present online.

Phylogenetic analysis
The DNA sequences were aligned with CLUSTALX 2.0 [36] and checked by eye for manual adjustment. The plastidial and the nuclear sequences were aligned separately to produce matrices that were later combined with the software combinex2_0.py ( [38]. The phylogenetic analysis was executed on both cpDNA (ndhF and trnL) and nuclear sequences (G3pdh). Maximum parsimony analysis was performed with PAUP* 4.0b1 [39] [40]. The genbank sequences of P. humistratum Görts & K. U. were used as outgroups both in the nuclear and the plastid genes matrix, following the previous phylogenetic analysis by Smith et al. (2008) [9]. This sequences used as outgroup resulted belonging to the sister clade with respect to the clade containing the African species and the other related clades. References of the other species used in the analysis are summarized (with genbank codes) in Table 1 in Smith et al. [9]. All characters had equal weight and unordered state transitions. Gaps were coded with the "simple indel coding" model [38], with the software Gapcoder [42] and added to the final matrix after the DNA sequences as in Papini et al. [43].
The evolutionary model implemented in Mrbayes for treating gaps was the same as that proposed by Lewis et al. [44] for treating morphological data, the Mk model.
We used MrMODELTEST 2.0 [45] to choose the best evolutionary model of DNA sequences on the basis of the Akaike information criterion [46]. The best model was used as settings with MrBayes 3.2.7 [47] for Bayesian Inference. A maximum likelihood (ML) phylogenetic analysis was carried out with RaxML [48]and the resulting trees were edited with Figtree [49]. We mapped the support on the tree branches with the results of the Bayesian phylogenetic analysis after removing the first trees with low likelihood values as "burn-in", as in Papini et al. [50] [51]. The remaining trees were used to produce a 50% majority-rule consensus tree in which the percentage indicated on branches was used as a measure of the Bayesian posterior probability.

Karyological analysis
Chromosomes images were obtained from somatic mitoses recorded from root tips of only one plant living in a pot. The procedure was the same as in Mosti et al. [52] and Mousavi et al. [53], with a pretreatment in 8-hydroxychinoline and fixation in Carnoy. Then the material was hydrolyzed in HCl and then stained with Lacto-propionic-orceine.
We observed metaphase plates of meristematic cells, with the technique of fresh squashes of root tips. Chromosome counts were made during direct observations with the microscope, and later recounted on enlarged digital images. Images were recorded with a microscope Leica DM RB Fluo.

Conclusions
The surprising discrepancy between the nuclear and the plastid phylogeny could be explained with an ancestral introgression event due probably to pollen contribution from an ancestor of the African mainland P. guineense towards the ancestor of P. tsarasotrae. The presence of possible paralogs Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 October 2020 of the nuclear gene G3pdh, clustering together with more distantly related Asian species lead to the hypothesis that a second more ancient hybridation/introgression event would have occurred between south Asian species and the ancestor of the Malagasian species. The chromosome numbers observed in the Malagasian species would confirm different evolutionary history.
Further studies about the karyotypes of the Malagasy species, the african P. guineense and P. borbonense will be necessary together with the investigation of the possible presence of short sequences in P. borbonense.