George P. Tiley, Akanksha Pandey, Rebecca T. Kimball, Edward L. Braun, J. Gordon Burleigh. 2020: Whole genome phylogeny of Gallus: introgression and data-type effects. Avian Research, 11(1): 7. DOI: 10.1186/s40657-020-00194-w
Citation: George P. Tiley, Akanksha Pandey, Rebecca T. Kimball, Edward L. Braun, J. Gordon Burleigh. 2020: Whole genome phylogeny of Gallus: introgression and data-type effects. Avian Research, 11(1): 7. DOI: 10.1186/s40657-020-00194-w

Whole genome phylogeny of Gallus: introgression and data-type effects

More Information
  • Corresponding author:

    J. Gordon Burleigh, gburleigh@ufl.edu

  • Received Date: 05 Nov 2019
  • Accepted Date: 06 Mar 2020
  • Available Online: 24 Apr 2022
  • Publish Date: 16 Mar 2020
  • Background 

    Previous phylogenetic studies that include the four recognized species of Gallus have resulted in a number of distinct topologies, with little agreement. Several factors could lead to the failure to converge on a consistent topology, including introgression, incomplete lineage sorting, different data types, or insufficient data.

    Methods 

    We generated three novel whole genome assemblies for Gallus species, which we combined with data from the published genomes of Gallus gallus and Bambusicola thoracicus (a member of the sister genus to Gallus). To determine why previous studies have failed to converge on a single topology, we extracted large numbers of orthologous exons, introns, ultra-conserved elements, and conserved non-exonic elements from the genome assemblies. This provided more than 32 million base pairs of data that we used for concatenated maximum likelihood and multispecies coalescent analyses of Gallus.

    Results 

    All of our analyses, regardless of data type, yielded a single, well-supported topology. We found some evidence for ancient introgression involving specific Gallus lineages as well as modest data type effects that had an impact on support and branch length estimates in specific analyses. However, the estimated gene tree spectra for all data types had a relatively good fit to their expectation given the multispecies coalescent.

    Conclusions 

    Overall, our data suggest that conflicts among previous studies probably reflect the use of smaller datasets (both in terms of number of sites and of loci) in those analyses. Our results demonstrate the importance of sampling large numbers of loci, each of which has a sufficient number of sites to provide robust estimates of gene trees. Low-coverage whole genome sequencing, as we did here, represents a cost-effective means to generate the very large data sets that include multiple data types that enabled us to obtain a robust estimate of Gallus phylogeny.

  • The Horned Larks (Eremophila alpestris) are a group of birds of wide distribution and more than 40 subspecies had been described (Drovetski et al. 2014). These larks represent an interesting model of study regarding the speciation process, as there are important morphological, biogeographic, and genetic differences between clades and they seem to have rapid phenotypic changes in response to habitat modification (Mason and Unitt 2018; Ghorbani et al. 2019).

    In 2013, a taxonomic study on the Alaudidae (larks) family using molecular techniques indicated that the Horned Larks could be considered as a species complex (where speciation boundaries are unclear and potentially different species are grouping under the same scientific name) (Alström et al. 2013). Later in 2014, a molecular analysis with samples of different geographic locations of the Horned Larks reinforced the need to split these birds into multiple species and proposed its division in seven different species, including six Palearctic (elwesi, bilopha, penicillata, atlas, brandti, and flava lineages) and one Nearctic (alpestris lineage) clades (Drovetski et al. 2014). However, a research in 2019 proposed the division of the Horned Larks into only four different species or clades and seven subclades, composed of (1) the Common Horned Lark (E. alpestris), subdivided into the Nearctic and the North Palearctic subclades (brandti, and flava lineages); (2) the Mountain Horned Lark (E. penicillata), subdivides into the penicillata and the atlas subclades; (3) the Temminck's Lark (E. bilopha); and (4) the Himalayan Horned Lark (E. longirostris), subdivided into the longirostris and elwesi subclades (Ghorbani et al. 2019).

    Nevertheless, there is a lack of information on several subspecies of the Horned Lark, so the taxonomy of this group of birds is still uncertain. The great majority of the Horned Larks have located in the Holarctic region, and both the Nearctic and Palearctic populations present a wide distribution and a high number of subspecies (several with migratory patterns). Nonetheless, in the Neotropics, there is only one particular and sedentary subspecies (E. a. peregrina) that is restricted to the high plateau of the Colombian Eastern Andes (Beason 1995; Cadena 2002; Valencia and Armenteras 2004), and it is endangered (EN) due to the urbanization and displacement of the native grasses (Renjifo et al. 2016; Zuluaga-Bonilla and Macana-García 2016). This bird is among the subspecies that have not been sampled in the genetic studies of the Horned Larks. So, there is no approximation of its speciation process or taxonomic status.

    Also, bioacoustics divergence analyses have not been included in the speciation studies of the Horned Larks. Suchinformation has proved to be useful in the analysis of speciation process and taxonomic revision of cryptic passerines (Toews and Irwin 2008). Song recognition is an important feature for birds evolutionary processes, as acoustic divergences preclude the conspecific recognition, causing barriers in gene flow, prevention of hybridization and finally promoting the speciation (Slabbekoorn and Smith 2002; Derryberry et al. 2018). Thus, analysis of song divergence and playback response has been important to study the taxonomic relation of birds (Alström and Ranft 2003). This approach has identified significant differences in song parameters and behavioral response to playback experiments among some subspecies of passerine birds, revealing significant premating reproductive isolation between populations (Dingle et al. 2010; Pieplow and Francis 2011; de Oliveira Gordinho et al. 2016; Slender et al. 2018). Further, evaluation of stimuli response in 72 pairs of closely related but geographically distant (allopatric) Neotropical passerines concluded that playback experiments are a suitable option in the study of the speciation process in birds with taxonomic uncertainty (Freeman and Montgomery 2017).

    We aimed to evaluate the differences in songs and behavioral responses of the Neotropical Horned Lark compared to the Nearctic subclade to expand the knowledge on the evolutionary history of this diverse and complex group of birds. To achieve this, we evaluated potential differences in the song's traits. Then, we tested if these variations in the songs could be biologically relevant by evaluating their effect on the behavioral response using playback experiments.

    The study was done in La Copa reservoir, a montane ecosystem, located at an elevation of 2660 m in the central Andean mountains at the municipality of Toca in the Boyacá state (Fig. 1). We selected this area because it harbors the largest populations of the Neotropical or Colombian Horned Lark and has reports of nesting (Valencia and Armenteras 2004; Botía-Becerra and Echeverry-Galvis 2010). A previous study by our research group showed that at the beginning of the reproductive season, both males and females presented an intense vocal activity and there were some fights between adult males (Arias-Sosa et al. 2021). Later the larks were observed mainly in couples and both parents were involved in the care of the nest and offspring, indicating that this bird is at least seasonally monogamous (Arias-Sosa et al. 2021).

    Figure  1.  Study area with the recording and playback points

    In North America, the Horned Larks tend to define territories during the reproductive season that are defended by the males (Beason and Franks 1974; Cannings and Threlfall 1981; DeGraaf and Yamasaki 2001). Nevertheless, in our study area, there were fights between adult males but they didn't define territories during the mating season; because the area was small and the population density was high, the larks were hinderedto have independent well-defined territories (Arias-Sosa et al. 2021).

    As we aimed to make playback experiments, we implemented a preliminary test with some recordings of males and females obtained during field observations. In this test, we realized that both males and females responded to the playback similarly, emitting responding songs and getting closer to the speaker. We distinguished males and females based on the distinctive facial marks of the adult individuals. About 60% of responses were from females and 40% from males. The majority of playback experiments only take into account male response, based on the assumption that the female response (mate choice) would follow the same pattern, but this may not be true in all cases (Matessi et al. 2000; Dingle et al. 2010; Pegan et al. 2015; de Oliveira Gordinho et al. 2016). Considering the female response is particularly important in the tropical region whereit is more common that female birds present a high song activity and response to playback stimuli (Slater and Mann 2004; Illes and Yunes-Jimenez 2009), we thereforedecided to include both sexes in the experiments.

    Also, during the playback preliminary test, we realized that multiple individuals responded during the same experiment (as the Horned Larks in the zone presented a gregarious behavior). It was difficult to differentiate the individual response due to the small size of the birds and because they usually move fast, which makes challenging to distinguish which individual of the group was emitting the response singing. Thus, we decided to implement the playback with all the birds in a radius of 25 m per experiment rather than individual analysis.

    During the preliminary experiments, we also observed the time it takes for the birds to return to their normal behavior after the acoustic stimuli. We observed that after 1 min the birds returned to look for food on the ground and move away from the speaker. Thus, we established that 2 min of silence between experiments would be adequate to avoid interference between them.

    As we mentioned previously during the preliminary studies, we could recognize the sex of several birds using binoculars at a close to medium distance (2‒10 m) and by following a single individual. However, when we tried to measure the playback response at a far distance (30 m, required for the playback experiments), it was difficult to accurately recognize the sex of each bird. Because the identification of the sexual dimorphism in this species requires great attention to particular marks on the face and chest, we thuscould not perform a separate sex analysis.

    The songs of the Neotropical Horned Lark were obtained by field recording in La Copa Reservoir from November 2017 to May 2018. The recordings were obtained from natural vocalizations using a unidirectional microphone (Takstar SGC 568) and were stored in a digital device (Tascam DR 07 MK II MDR12). The same device was used for all recordings, and files were converted to WAV 16-bit, 44.1 kHz. To avoid pseudoreplication we took the recordings in line transect across the 6 points studied (Fig. 1). We took recordings of one singing lark, and then we moved on to take a record of another individual. Overall, we got 40 recordings from E. a. peregrina.

    The most recent work in the taxonomy of the Horned Larks proposed it split into four different species or clades, three of them distributed in the western Palearctic and one distributed in the North Palearctic and the Nearctic (Ghorbani et al. 2019). This last species of the Horned Lark was named E. alpestris (sensu stricto) and is subdivided into two subclades: the Nearctic and the North Palearctic subclade. We chose to compare the songs of E. a. peregrina with the Nearctic subclade because it is the geographically closest one and the one with more recordings available. Also because the previous genetics studies indicated a single colonization process of the New-world and a common evolutionary history of the New-world larks (although they only considered some North American populations) (Ghorbani et al. 2019).

    We seached for recordings of the Nearctic E. alpestris in the Macaulay Library (Cornell Lab of Ornithology 2018; https://www.macaulaylibrary.org/) and Xeno-canto (Xeno-canto Foundation 2018; https://www.xeno-canto.org/). We had into account only recordings longer than 30 s. These files were obtained in MP3 format and converted to WAV 16-bit, 44.1 kHz. We retrieved 88 recordings.

    All recordings were digitized and analyzed using AvisoftSAS Lab Lite software (Sound Analysis and Synthesis Laboratory Version 5.2.12.01 December 2017) (Turčoková et al. 2010; Hamao et al. 2013). From an initial analysis of the recordings, we found two types of vocalizations of the Horned Lark, a short call composed of two or three notes with a length of 0.5 s, and another with a mean of 10 notes and a length of about 2 s. This last one was considered the principal song and was the only one used in the next analysis. We selected this song because it presented more notes and traits that allow better analysis. We observed that the short call is used to communicate between close individuals, especially between the parents with their offspring.

    To eliminate external interferences like wind, songs of other birds, and human noise we applied the infinite impulse response (IIR) and finite impulse response (FIR) filters to each recording. These are common and effective mathematical algorithms used in digital acoustic management and in birds' bioacoustics studies to reduce noise and avoid artifacts (Leonardo 2004; Mockford and Marshall 2009; Rimell et al. 2015). We selected the adequate cut-off frequencies for each song based on the sonogram as it varied between recordings (Leonardo 2004; Mockford and Marshall 2009).

    From the original 40 Colombian and 88 Nearctic recordings retrieved, we excluded 20 and 60 recordings respectively because they only contained short calls, had incomplete songs, presented too much background noise, or had a bad sound quality. Thus, we continued the analysis with only 20 Colombian and 28 Nearctic recordings (Additional file 1: Table S1).

    For each recording, we extrapolated the songs and excluded those that presented significant noise based on the clarity of the spectrogram. From this analysis, we got 71 songs (obtained from the 20 recordings) of the Colombian lark and 78 songs (obtained from the 28 recordings) of the Nearctic larks. We analyzed these songs using the software AvisoftSAS (with a – 5 to – 10 dB threshold) to get the number of notes, the number of notes per second, song length, the time between notes (average), the highest frequency, and the lowest frequency. We selected these song features because they have been described as parameters of divergence between passerines (Dingle et al. 2008, 2010).

    All statistical analyses were performed using the free software R (3.6 version) and the package R Commander (2.5–3 version), with a p < 0.05 as significant (R Core Team 2020). We obtained several songs from the same recording (and potentially the same individual) and we observed interesting variation between them. Thus, we made two separate analyses: (1) using all the obtained songs to get an analysis that included the individual variation and (2) using one randomly selected song per recording to get an analysis not affected by individual pseudoreplication.

    As the variables did not present normal distribution (evaluated by Shapiro–Wilk and Kolmogorov–Smirnov tests), the differences in acoustic features were analyzed by the nonparametric Mann–Whitney U test (Turčoková et al. 2010). Later we implemented a principal component analysis (PCA), using Past software (Hammer et al. 2001) to combine the acoustic variables into a single one named "song divergence". It corresponded to the first axis of the PCA (PC1) that presented an eigenvalue > 1 and was evaluated by the ANOVA test as it was normally distributed. The PCA was done because there was a probable correlation between the variables based on Bartlett's Test of Sphericity (BTS) and Kaiser-Meyer Olkin test (KMO) in both analyses: using all songs (BTS p < 001, KMO = 0.48) and one song per recording (BTS p < 001, KMO = 0.46).

    Finally, we plot the PC1 vs PC2 using a scatterplot to see the grouping of the Nearctic and Neotropical songs. As the Nearctic songs come from different locations in North America, we show the locations in the scatterplot classified in: Central-West (CW), North-East (NE), North-West (NW), South-West (SW), and South USA-Mexico (SM). This classification was based on the subdivision of the Nearctic Horned Larks proposed by Ghorbani et al. (2019) (except South USA-Mexico that was not included in this previous study).

    The playback experiments evaluated the degree of recognition and response strength by individuals of E. a. peregrina exposed to the conspecific and Neotropical songs stimuli. The experiments were carried out in the northern part of La Copa Reservoir, where our research group registered the highest densities of E. a. peregrina and several nests (Arias-Sosa et al. 2021).

    In each playback experiment, we tested the response to three stimuli: (1) songs of the conspecific populations (Neotropical), (2) songs of the allopatric populations (Nearctic), and (3) response to a different species (Rufous-collared Sparrow Zonotrichia capensis) as a control group. Z. capensis is a common bird of the Colombian Andes that cohabitates naturally in sympatry with the lark and present similar size and behavior (feeding on seeds and insects in grasslands). The songs from the Rufous-collared Sparrow were obtained from Xeno-canto recordings from Boyacá, Colombia. In the conspecific (Neotropical) stimuli, we used songs from other localities (the farthest away in the southern region of the reservoir) to avoid the familiarity of the larks with the local dialect.

    We implemented a total of 36 playback experiments. Each experiment consisted of the broadcasting of a recording, composed of 2 min mix of songs of each one of the three stimuli (Neotropical, Nearctic, and control) separated by 2 min of silence between each other. We prepared each recording with a different mix of songs to avoid potential bias (Pegan et al. 2015). To prepare each recording we had a set of 38 songs of each of the three stimuli. We randomly selected 10 songs per stimuli, added a random pause between songs, and randomly repeated them until complete the 2 min stimuli. All recordings were normalized in the range of 30–80 dB.

    We broadcast the recordings using an iPod Touch 4G 8 GB player through amplification with an EXTRA BASS™ XB32 portable speaker at 40–80 dB at 1 m, placed on the ground (because the Horned Larks stay most of the time foraging on the floor). We started the experiments only if the close birds were looking for food and not singing by themselves (indicators of a non-stimulated state), to reduce the risk of getting an intrinsic response (not related to the playback stimuli) (de Oliveira Gordinho et al. 2016; Slender et al. 2018).

    To reduce the effect of pseudoreplication the experiments were performed in different locations and with one-hour intervals to avoid interference between each other. To avoid potential bias the conspecific (local population) songs were always reproduced after the allopatric and control songs (Sagario and Cueto 2014). This because the conspecific song could overstimulate the response of the birds in the subsequent stimulus (Sagario and Cueto 2014). In the case of Nearctic and control treatment, we alternated its order. The response to the conspecific songs is unlikely to be affected by the reproduction order because as we observe in the preliminary observational study 2 min of silence is enough for the birds to return to their normal (unstimulated) behavior. During all the playback experiments we confirm that the target birds returned to their feeding behavior and get away from the speaker within the 2 min of silence.

    Once we started the experiment, we registered the behavioral responses at approximately 30 m from the speaker using binoculars. We followed the response during the 2 min of each stimulus and the subsequent 1 min of silence. As the Horned Larks in our study area presented a gregarious behavior that makes it difficult to take individual data, we retrieved the information of all responding larks per experiment. As we reproduce the 3 archives in a consecutive way the same number of individuals were exposed to each stimulus. In some of the experiments, we could observe that both sexes sang in response and approached the speaker; however, we could not differentiate the sex of most of the individuals and we were not able of making a sex-specific analysis. We evaluated the following response variables:

    1. The approach response: the number of birds that approached the speaker. This variable was subdivided into: approach within a 20-m radius and an approach within a 10-m radius. A higher number of larks approaching indicates a stronger response.

    2. The response vocalization: the number of songs emitted by the larks during the experiments on a radius of 25 m. A higher number of response songs indicates a stronger response.

    3. The latency to approach: the time (in seconds) that it took for the first birds to approached the speaker within a 5-m radius. When no bird approached, we reported the measure as 180 s (the 3 min of observation). A lower time is associated with a stronger response.

    These variables were selected because they have been described as consistent indicators of behavioral response to playback in passerine birds (Matessi et al. 2001; Jankowski et al. 2010; de Oliveira Gordinho et al. 2016; Freeman et al. 2016; Slender et al. 2018).

    The research was done under a research project made in collaboration with the "Corporación Autónoma Regional de Boyacá" the governmental entity that regulates environmental activities in the state. None individual was captured, manipulated, or suffered any damage.

    The differences in the playback variables were evaluated by the Kruskal–Wallis test with a Bonferroni post-hoc test. Because there is a probable correlation between the variables (BTS p < 001, KMO = 0.7), we used principal component analysis (PCA), using Past software to combinate the playback variables into a single one named "response intensity". This variable corresponded to the first axis of the PCA (PC1) with an eigenvalue > 1, and was evaluated by the ANOVA test with a Bonferroni post-hoc test as it was normally distributed. Also, we used a PC1 vs PC2 scatterplot to observe the difference in the response to the three stimuli.

    The songs of the Neotropical larks presented some general and notorious differences from those of the Nearctic populations when we observed the spectrograms (Additional file 2: Figure S1). The Neotropical songs were generally paused (with a longer pause between notes), shorter, and with a lower number of notes. While the Neotropical songs were longer, with more notes, and with a lower pause between notes.

    Both analyses (using all the retrieved songs and one song per recording) showed a significant difference in the number of notes, song length, the time between notes, highest frequency, and lowest frequency (Table 1; Fig. 2). On the other hand, the number of notes per second was significantly different in the analysis using all the retrieved songs (p < 0.001) but evidenced a notorious overlapping between the two groups (Fig. 2) and it was not significantly different when we take only one song per recording (p = 0.36). Thus, there seem to be differences in the six song traits. But the clearest divergent ones were the number of notes, the time between notes, and the frequency limits. Regarding these variables, the songs from the Neotropical Horned Lark presented a lower number of notes, a higher time between notes (longer pauses), and higher frequencies than the Nearctic songs (Table 1).

    Table  1.  Differences in Song features
    a. Analysis using all retrieved songs (n = 149)
    Variable Statistic result Neotropical (n = 71) Nearctic (n = 78)
    PCA F = 261.4; p < 0.001 NA NA
    Number of notes H (2) = 13.916; p < 0.001 7.75 ± 2.07 10.64 ± 2.59
    Number of notes per second H (2) = 0.8086; p = 0.3685 6.83 ± 2.8 7.30 ± 1.99
    Song length H (2) = 7.1654; p < 0.05 1.23 ± 0.35 1.54 ± 0.48
    Time between notes H (2) = 12.197; p < 0.001 0.62 ± 0.22 1.03 ± 0.53
    Highest frequency H (2) = 16.098; p < 0.001 7310 ± 441.17 6725 ± 448.56
    Lowest frequency H (2) = 28.251; p < 0.001 4230 ± 479.14 3250 ± 288.68
    b. Analysis using one song per record (n = 48)
    Variable Statistic result Neotropical (n = 20) Nearctic (n = 28)
    PCA F = 38.86; p < 0.001 NA NA
    Number of notes H (2) = 49.433; p < 0.001 7.93 ± 2.34 11.36 ± 2.59
    Number of notes per second H (2) = 11.868; p < 0.001 6.35 ± 2.18 7.49 ± 1.66
    Song length H (2) = 18.999;. p < 0.001 1.31 ± 0.36 1.57 ± 0.41
    Time between notes H (2) = 52.162; p < 0.001 0.53 ± 0.23 1.00 ± 0.42
    Highest frequency H (2) = 65.335; p < 0.001 7443.66 ± 365.57 6744.60 ± 466.45
    Lowest frequency H (2) = 94.058; p < 0.001 4245.07 ± 332.86 3256.76 ± 326.49
    The data are shown as mean±standard deviation. (a) Difference in the analysis using all retrieved songs. (b) Differences in the analysis using one song per record
     | Show Table
    DownLoad: CSV
    Figure  2.  Differences in the acoustic traits among the Nearctic and Neotropical Horned Larks. We showed a boxplot of each variable evaluated in: a the analysis using all retrieved songs, and b the analysis using one song per recording.* = p < 0.05; ** = p < 0.005; *** = p < 0.001

    For the PCA analysis, we found a significant difference in the PC1 component (named song divergence) between the Nearctic and Neotropical lark songs in both analyses (Fig. 3). Likewise, the scatterplot showed a clear differentiation between the Nearctic and Neotropical songs (Fig. 3). While, there was high overlap in the songs of the different geographic subdivisions of the Nearctic clade. This indicates that the Neotropical lark presents more divergent song traits than the observed within the different populations in North America (that have a recent genetic divergence). The loads and eigenvalues of the PCA are shown in Additional file 1: Tables S2–4.

    Figure  3.  Boxplot of the song divergence (PC1) between the Nearctic and Neotropical songs and PC1 vs PC2 scatterplot evaluated in: a the analysis using all retrieved songs, and b the analysis using one song per recording. The ellipses of the scatterplot were calculated automatically in the R software using the option of concentration ellipses (.5.8°). For the Nearctic songs, we show its geographic subclassification: Central-West (CW), North-East (NE), North-West (NW), South-West (SW), and South USA-Mexico (SM). *** = p < 0.001

    We still observed an important degree of approach response to the Nearctic stimuli but it was significantly lower (less lark approached) than the response to the conspecific songs within 25 m (p < 0.05) and 10 m radius (p < 0.001) (Fig. 4). Likewise, the larks exposed to the conspecific (Neotropical) stimuli emitted more reply songs and approached faster (p < 0.001) (Table 2). Compared to the experiments with the North American stimuli, there were about three times more larks approached to the 10 m radius in response to the conspecific song stimuli (3.8 ± 1.6 larks vs 1.3 ± 1.1 larks), they approached two times faster to a 5 m radius to the speaker (63.2 ± 27.7 s vs 121.3 ± 33.4 s) and emitted three times more reply songs (9.5 ± 2.9 songs vs 2.8 ± 1.4 songs). On the other hand, the response to the control experiment was minimal as expected (Fig. 4). This indicates that although the Colombian Horned Lark seems to respond to the Nearctic songs, the strength of this response is notoriously lower.

    Figure  4.  Differences in the playback response to the Nearctic, Neotropical, and control (Rufous-collared Sparrow) stimuli. a Boxplot of each variable evaluated. b Boxplot of the response intensity variable (PC1). Lower PC1 values indicate a stronger response as it reflects a faster and closer approach with more reply vocalizations. Also, we show the PC1 vs PC2 scatterplot. The ellipses of the scatter plot were calculated automatically in the R software using the option of concentration ellipses (.5.8 degree). ** = p < 0.005; *** = p < 0.001
    Table  2.  Differences in the response to playback experiments
    Variable Statistic result Neotropical (n=36) Nearctic (n=36) Control (n=36)
    PCA F=247;p < 0.001 NA NA NA
    Latency H (2)=85.839;p < 0.001 63.25±27.7 121.31±33.4 180±0
    Number of vocalizations H (2)=87.160;p < 0.001 9.56±2.9 2.78±1.42 0.67±0.72
    Approach response (10 m) H (2)=79.661;p < 0.001 3.83±1.65 1.28±1.09 0±0
    Approach response (20 m) H (2)=68.394;p < 0.001 4.94±1.62 3.58±1.98 0.31±0.62
    The data are shown as mean±standard deviation
     | Show Table
    DownLoad: CSV

    The PCA analysis showed a significant difference in the PC1 (response intensity) between the three groups (p < 0.001) and the scatterplot evidenced a notable separation between the response to the Nearctic and Neotropical stimuli (Fig. 4). The loads and eigenvalues of the PCA are shown in Additional file 1: Table S2–4.

    This is the first study that evaluated differences in songs of allopatric populations of E. alpestris, and the first speciation work in E. a. peregrina. The Colombian subspecies presented differences in the song traits and PCA analysis indicating some modifications in the dialect compared to the Nearctic subclade. Also, we evidence differences in behavioral response during playback experiments. The Colombian Horned Lark seems to respond to the song of the Nearctic larks but with a lower intensity than to local vocalizations.

    The differences found indicate that male Nearctic Horned Larks may be perceived as the least significant sexual competitors by the Colombian subspecies, and females prefer the local songs which could lead to a premating reproductive isolation (Searcy 1992; Ballentine et al. 2004; Balakrishnan and Sorenson 2006). However, as we could not make a distinction between male and female responses, further sex-specific studies are required to explore if the preference for local songs occurs in the same degree between males and females. The difference in the response between local and Nearctic songs occurs because birds learn their songs from neighboring conspecifics and distant populations tend to lost similarity, especially in sedentary species with breeding site fidelity (like E. a. peregrina) (Morton 1987; Kroodsma 2002).

    The combination of vocal traits analysis and playback experiments has been described as a robust measure of reproductive isolation (Pegan et al. 2015; Freeman and Montgomery 2017). This approach has been used as part of the criteria for the delimitation of species limits and in speciation studies, and joining the mitochondrial DNA analysis, are among the most divergent characteristics between bird species (Alström and Olsson 1992; Alström et al. 2008). Some examples of this use are the analysis of subspecies of White-breasted (Henicorhina leucosticte) and Gray-breasted Wood-Wren (H. leucophrys) that presented differences in song traits and an asymmetric playback response, described as a premating reproductive barrier that suggested a complex of multiple biological species (Dingle et al. 2010; Pegan et al. 2015). Research in recently separated subspecies of Reed Bunting (Emberiza schoeniclus) showed differences in song traits, but less clear differences in playback response, indicating that this population experiences an incipient or early stage of speciation (Matessi et al. 2000; Gordinho et al. 2015; de Oliveira Gordinho et al. 2016). While in subspecies of the Buff-breasted Earthcreeper (Upucerthia validirostris) no differences in playback response were found and the researchers recommended keep them as the same species (Areta and Pearman 2013). However other important features like DNA, morphometrics, plumage traits, biometrics, egg characteristics, and biogeography need to be considered in the species delimitation (Alström and Olsson 1992; Alström et al. 2008).

    In our case, we found that the Neotropical larks still respond to the Nearctic songs but in a significantly lower degree, and there are important differences in the song traits between these two groups. Further, the Neotropical lark songs are more clearly differentiated than the subgroups of Horned Lark in the Nearctic (Fig. 3). The fact that the Neotropical larks still respond to the Nearctic songs is consistent with the conclusion of previous studies that indicate that the New-world Horned Larks subspecies are recently divergent between each other and are part of the same species or clade (Ghorbani et al. 2019). However, the high degree of differentiation in the songs and lower response to playback seem to indicate a particular evolutionary history of the Neotropical Horned Larks.

    Besides the song differences that we found, the Neotropical lark presents other important ecological and biogeographic differences in comparison to the Nearctic populations and the North Palearctic larks that compose the alpestris clade. E. a peregrina is highly isolated (similar to an island type isolation) from other Horned Larks groups, have a very restricted distribution area (160 km2), and only inhabited the high Plateau of the Eastern Andes between the 2500‒3200 m a.s.l. (Valencia and Armenteras 2004). While, the larks in North America and the North Palearctic have a wide distribution (of most of the continent) and exploit different habitats that go from the coast's region to the high mountains. Also, the clutch size of the Neotropical Horned Lark (two eggs) is lower than the observed in the Nearctic populations (2‒5 eggs) (Beason and Franks 1974; Cannings and Threlfall 1981; Garry Oak Ecosystems Recovery Team 2003; Camfield et al. 2010). Further, these characteristics are more similar to the clades/species in the South Palearctic that present an isolated and reduced distribution, associated with mountain habitats. Like the Himalayan Horned Lark clade (E. longirostris), which presents a restricted distribution area in the high Qinghai–Tibetan Plateau and clutch size of 1‒3 eggs (Ghorbani et al. 2019). Unfortunately, there are too few recordings of this Horned Lark clade to compare with the Neotropical songs.

    Thus, the Neotropical Horned Lark (E. a. peregrina) has important bioacoustic, biogeographic, and ecological differences compared to the Nearctic subclade that makes unclear its evolutionary relation. Genetic and morphological comparison of the Neotropical lark with the established clades of the Horned Lark is required to assess its taxonomic status and evolutionary history. With the current information, we can propose two hypotheses: (1) the Neotropical subspecies follows the evolutionary history observed by Ghorbani et al. (2019) and Drovetski et al. (2014), who established that the New-world Horned Larks subspecies come from a single colonization event from the Palearctic and have a common evolutionary history (monophyly); (2) the Neotropical Horned Lark has a separate evolutionary history from the Nearctic clade and could be originated from a different colonization event. However, these hypotheses need to be evaluated in future studies using DNA sampling.

    One limitation of our study was that we did not evaluate the reciprocal responses in the Nearctic populations. So, we suggest that future researches focus on the response of Nearctic Horned Larks to the vocalizations of Neotropical birds. Also, it would be interesting to analyze the song traits of the other three clades of Horned Larks in the south Palearctic and make playback experiments using these songs to have more clear information of the song divergence in contrast to the known genetic variance (Drovetski et al. 2014; Ghorbani et al. 2019). Also, we were incapable of recognizing the sex of the tested larks in the field experiments, which implies a lack of information on the specific response of male and female larks. Although in the preliminary observational study, we realized that both males and females respond to playback experiments, some research had found that the behavioral response of birds exposed to acoustic signals can be contrasting between males and females (Snijders et al. 2017; Bircher et al. 2020). Likewise, the response of the males could be affected by the presence of a female partner around. Thus, further studies are necessary to determine the exact degree of song differentiation by males and females of the Horned Larks.

    We evidenced important differences in the song traits between the Nearctic and Neotropical Horned Larks. Further, the song divergence of the Neotropical larks was more pronounced than those found within the Nearctic subgroups. In the playback analysis, we found that the Neotropical lark still presents an important degree of recognition to the Nearctic songs, but the strength of the response was significantly lower. Besides this acoustic difference, the neotropical Horned Larks display particular biogeographic, behavioral, and reproductive traits that are more related to some clades in the South Palearctic. Thus, the relation between the Nearctic and Neotropical larks is uncertain and further genetic studies are required to establish if the Neotropical larks had the same colonization and dispersion history described in the New World subspecies of North America or experienced a separate evolutionary history.

    The online version contains supplementary material available at https://doi.org/10.1186/s40657-021-00251-y.

    Additional file 1: Table S1. Records used in the song traits analysis.

    Table S2. PCA analysis of the song's traits (using all the retrieved songs).

    Table S3. PCA analysis of the song's traits (using one song per record).

    Table S4. PCA analysis of the play back response.

    Additional file 2: Figure S1. Representative spectrograms of two records (individuals) of the Nearctic and Neotropical Horned Larks.

    Special thanks to professor Pablo Emilio Rodríguez Africano who managed the resources to develop this research and coordinated the data collection. To the "Corporación Autónoma Regional de Boyacá (CORPOBOYACÁ)" and the "Universidad Pedagógica y Tecnológica de Colombia UPTC", for the financial and logistical support.

    JAV-P Took and processed the data, LAA-S wrote the manuscript and made the statistical analysis, and CR-M revised the manuscript. All authors read and approved the final manuscript.

    The datasets used in the present study are available from the corresponding author on reasonable request.

    This work was done under a research project in collaboration with the "Corpo- ración Autónoma Regional de Boyacá", the environmental entity that regulates wildlife management in the state. The experiments comply with the current Colombian laws regarding research in the biological field.

    Not applicable.

    The authors declare that they have no competing interests.

  • Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403-10.
    Armstrong MH, Braun EL, Kimball RT. Phylogenetic utility of avian ovomucoid intron G: a comparison of nuclear and mitochondrial phylogenies in the Galliformes. Auk. 2001;118:799-804.
    Berlin S, Ellegren H. Evolutionary genetics: clonal inheritance of avian mitochondrial DNA. Nature. 2001;413:37-8.
    Betancur-R R, Li C, Munroe TA, Ballesteros JA, Ortí G. Addressing gene tree discordance and non-stationarity to resolve a multi-locus phylogeny of the flatfishes (Teleostei: Pleuronectiformes). Syst Biol. 2013;62:763-85.
    Bolívar P, Guéguen L, Duret L, Ellegren H, Mugal CF. GC-biased gene conversion conceals the prediction of the nearly neutral theory in avian genomes. Genome Biol. 2019;20:5.
    Braun EL, Kimball RT. Examining basal avian divergences with mitochondrial sequences: model complexity, taxon sampling and sequence length. Syst Biol. 2002;51:614-25.
    Braun EL, Cracraft J, Houde P. Resolving the avian tree of life from top to bottom: the promise and potential boundaries of the phylogenomic era. In: Kraus RHS, editor. Avian genomics in ecology and evolution—from the lab into the wild. Cham: Springer; 2019. p. 151-210.
    Brown JM, Thomson RC. Bayes factors unmask highly variable information content, bias, and extreme influence in phylogenomic analyses. Syst Biol. 2016;66:517-30.
    Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinform. 2008;10:421.
    Cantarel BL, Korf I, Robb SMC, Parra G, Ross E, Moore B, et al. MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Res. 2008;18:188-96.
    Chen MY, Liang D, Zhang P. Phylogenomic resolution of the phylogeny of laurasiatherian mammals: exploring phylogenetic signals within coding and noncoding sequences. Genome Biol Evol. 2017;9:1998-2012.
    Chojnowski J, Kimball RT, Braun EL. Introns outperform exons in analyses of basal avian phylogeny using clathrin heavy chain genes. Gene. 2008;410:89-96.
    Conant GC, Lewis PO. Effects of nucleotide composition bias on the success of the parsimony criterion in phylogenetic inference. Mol Biol Evol. 2001;18:1024-33.
    Cunningham F, Amode MR, Barrell D, Beal K, Billis K, Brent S, et al. Ensembl 2015. Nucleic Acids Res. 2015;43:D662-9.
    Doyle JJ. Trees within trees: genes and species, molecules and morphology. Syst Biol. 1997;46:537-53.
    Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792-5.
    Edwards SV, Cloutier A, Baker AJ. Conserved nonexonic elements: a novel class of marker for phylogenetics. Syst Biol. 2017;66:1028-44.
    Eriksson J, Larson G, Gunnarsson U, Bed'hom B, Tixier-Boichard M, Strömstedt L, et al. Identification of the Yellow Skin gene reveals a hybrid origin of the domestic chicken. PLoS Genet. 2008;4:e1000010.
    Faircloth BC, McCormack JE, Crawford NG, Harvey MG, Brumfield RT, Glenn TC. Ultraconserved elements anchor thousands of genetic markers spanning multiple evolutionary timescales. Syst Biol. 2012;61:717-26.
    Felsenstein J. Confidence limits on phylogenies: an approach using the bootstrap. Evolution. 1985;39:783-91.
    Gatesy J, Springer MS. Phylogenetic analysis at deep timescales: unreliable gene trees, bypassed hidden support, and the coalescence/concatalescence conundrum. Mol Phylogenet Evol. 2014;80:231-66.
    Gee H. Evolution: ending incongruence. Nature. 2003;425:782.
    Gill F, Donsker D. IOC world bird list (v9.2). 2019. .
    Hill WG, Robertson A. The effect of linkage on limits to artificial selection. Genet Res. 1966;8:269-94.
    Hosner PA, Faircloth BC, Glenn TC, Braun EL, Kimball RT. Avoiding missing data biases in assembling the landfowl tree of life (Aves: Galliformes). Mol Biol Evol. 2016;33:1110-25.
    Hosner PA, Tobias JA, Braun EL, Kimball RT. How do seemingly non-vagile clades accomplish trans-marine dispersal? Trait and dispersal evolution in the landfowl. Proc Roy Soc Lond B. 2017;284:20170210.
    Imsland F, Feng C, Boije H, Bed'Hom B, Fillon V, Dorshorst B, et al. The Rose-comb mutation in chickens constitutes a structural rearrangement causing both altered comb morphology and defective sperm motility. PLoS Genet. 2012;8:e1002775.
    International Chicken Genome Sequencing Consortium (ICGSC). Sequence and comparative analysis of the chicken genome provide unique perspectives on vertebrate evolution. Nature. 2004;432:695-716.
    Jarvis ED, Mirarab S, Aberer AJ, Li B, Houde P, Li C, et al. Whole-genome analyses resolve early branches in the tree of life of modern birds. Science. 2014;346:1320-31.
    Jeffroy O, Brinkmann H, Delsuc F, Philippe H. Phylogenomics: the beginning of incongruence? Trends Genet. 2006;22:225-31.
    Johnsgard PA. The pheasants of the world. 2nd ed. Oxford: Oxford University Press; 1999. p. 92-9.
    Kan XZ, Li XF, Lei ZP, Chen L, Gao H, Yang ZY, et al. Estimation of divergence times for major lineages of galliform birds: evidence from complete mitochondrial genome sequences. Afric J Biotech. 2010a;9:3073-8.
    Kan XZ, Yang JK, Li XF, Chen L, Lei ZP, Wang M, et al. Phylogeny of major lineages of Galliform birds (Aves: Galliformes) based on complete mitochondrial genomes. Genet Mol Res. 2010b;9:1625-33.
    Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30:772-80.
    Katzman S, Kern AD, Bejerano G, Fewell G, Fulton L, Wilson RK, et al. Human genome ultraconserved elements are ultraselected. Science. 2007;317:915.
    Kimball RT, Braun EL. A multigene phylogeny of Galliformes supports a single origin of erectile ability in non-feathered facial traits. J Avian Biol. 2008;39:438-45.
    Kimball RT, Braun EL. Does more sequence data improve estimates of galliform phylogeny? Analyses of a rapid radiation using a complete data matrix. PeerJ. 2014;2:e361.
    Kimball RT, Braun EL, Zwartjes P, Crowe TM, Ligon JD. A molecular phylogeny of the pheasants and partridges suggests that these lineages are not monophyletic. Mol Phylogenet Evol. 1999;11:38-54.
    Kimball RT, Mary CM, Braun EL. A macroevolutionary perspective on multiple sexual traits in the Phasianidae (Galliformes). Int J Evol Biol. 2011;2011:423938.
    Kimball RT, Wang N, Heimer-McGinn V, Ferguson C, Braun EL. Identifying localized biases in large datasets: a case study using the avian tree of life. Mol Phylogenet Evol. 2013;69:1021-32.
    Lawal RA, Al-Atiyat RM, Aljumaah RS, Silva P, Mwacharo JM, Hanotte O. Whole-genome resequencing of red junglefowl and indigenous village chicken reveal new insights on the genome dynamics of the species. Front Genet. 2018;9:264.
    Li D, Li Y, Li M, Che T, Tian S, Chen B, et al. Population genomics identifies patterns of genetic diversity and selection in chicken. BMC Genomics. 2019;20:263.
    Li L, Stoeckert J Jr, Roos DS. OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 2003;13:2178-89.
    Ligon JD, Kimball RT, Merola-Zwartjes M. Mate choice in red junglefowl: the issues of multiple ornaments and fluctuating asymmetry. Anim Behav. 1998;55:41-50.
    Maddison WP. Gene trees in species trees. Syst Biol. 1997;46:523-36.
    Marcais G, Kingsford C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioninformatics. 2011;27:764-70.
    Marcais G, Yorke JA, Zimin A. QuorUM: an error corrector for Illumina reads. PLoS ONE. 2015;10:e0130821.
    McCarthy EM. Handbook of avian hybrids of the world. New York: Oxford University Press; 2006.
    Meiklejohn KA, Danielson MJ, Braun EL, Faircloth BC, Glenn TC, Kimball RT. Incongruence among different mitochondrial regions: a case study using complete mitogenomes. Mol Phylogenet Evol. 2014;78:314-23.
    Meiklejohn KA, Faircloth BC, Glenn TC, Kimball RT, Braun EL. Analysis of a rapid evolutionary radiation using ultraconserved elements: evidence for a bias in some multispecies coalescent methods. Syst Biol. 2016;65:612-27.
    Mendes FK, Hahn MW. Gene tree discordance causes apparent substitution rate variation. Syst Biol. 2017;65:711-21.
    Mirarab S, Bayzid MS, Boussau B, Warnow T. Statistical binning improves species tree estimation in the presence of gene tree incongruence. Science. 2014;346:1250463.
    Moore WS. Inferring phylogenies from mtDNA variation: mitochondrial-gene trees versus nuclear-gene trees. Evolution. 1995;49:718-26.
    Mugal CF, Weber CC, Ellegren H. GC-biased gene conversion links the recombination landscape and demography to genomic base composition: gC-biased gene conversion drives genomic base composition across a wide range of species. BioEssays. 2015;37:1317-26.
    Nishibori M, Shimogiri T, Hayashi T, Yasue H. Molecular evidence for hybridization of species in the genus Gallus except for Gallus varius. Anim Genet. 2005;36:367-75.
    Pamilo P, Nei M. Relationships between gene trees and species trees. Mol Biol Evol. 1988;5:568-83.
    Patel S, Kimball RT, Braun EL. Error in phylogenetic estimation for bushes in the tree of life. J Phylogenet Evol Biol. 2013;1:110.
    Peterson AT, Brisbin IL. Genetic endangerment of wild Red Junglefowl Gallus gallus? Bird Conserv Int. 1998;8:387-94.
    Reddy S, Kimball RT, Pandey A, Hosner PA, Braun MJ, Hackett SJ, et al. Why do phylogenomic data sets yield conflicting trees? data type influences the avian tree of life more than taxon sampling. Syst Biol. 2017;66:857-79.
    Rosenberg NA. The probability of topological concordance of gene trees and species trees. Theor Pop Biol. 2002;61:225-47.
    Saitou N, Nei M. The number of nucleotides required to determine the branching order of three species, with special reference to the human-chimpanzee-gorilla divergence. J Mol Evol. 1986;24:189-204.
    Sayyari E, Mirarab S. Fast coalescent-based computation of local branch support from quartet frequencies. Mol Biol Evol. 2016;33:1654-68.
    Shen XX, Hittinger CT, Rokas A. Contentious relationships in phylogenomic studies can be driven by a handful of genes. Nat Ecol Evol. 2017;1:0126.
    Shen YY, Dai K, Cao X, Murphy RW, Shen XJ, Zhang YP. The updated phylogenies of the Phasianidae based on combined data of nuclear and mitochondrial DNA. PLoS ONE. 2014;9:e95786.
    Shen YY, Liang L, Sun YB, Yue BS, Yang XJ, Murphy RW, et al. A mitogenomic perspective on the ancient, rapid radiation in the Galliformes with an emphasis on the Phasianidae. BMC Evol Biol. 2010;10:132.
    Smit AFA, Hubley R, Green P. RepeatMasker Open-4.0. . Accessed 8 Dec 2015.
    Solís-Lemus C, Ané C. Inferring phylogenetic networks with maximum pseudolikelihood under incomplete lineage sorting. PLoS Genet. 2016;12:e1005896.
    Solís-Lemus C, Bastide P, Ané C. PhyloNetworks: a package for phylogenetic networks. Mol Biol Evol. 2017;34:3292-8.
    Song S, Liu L, Edwards SV, Wu S. Resolving conflict in eutherian mammal phylogeny using phylogenomics and the multispecies coalescent model. Proc Natl Acad Sci USA. 2012;109:14942-7.
    Springer MS, Gatesy J. The gene tree delusion. Mol Phylogenet Evol. 2016;94:1-33.
    Springer MS, Gatesy J. Retroposon insertions within a multispecies coalescent framework suggest that ratite phylogeny is not in the 'Anomaly Zone'. bioRxiv. 2019. .
    Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312-3.
    Stanke M, Schöffmann O, Morgenstern B, Waack S. Gene prediction in eukaryotes with a generalized hidden Markov model that uses hints from external sources. BMC Bioinform. 2006;7:62.
    Stein RW, Brown JW, Mooers AØ. A molecular genetic time scale demonstrates Cretaceous origins and multiple diversification rate shifts within the order Galliformes (Aves). Mol Phylogenet Evol. 2015;92:155-64.
    Swofford DL. PAUP*: phylogenetic analysis using parsimony (*and other methods). Version 4.0a166. 2019.
    Tamashiro RA, White ND, Braun MJ, Faircloth BC, Braun EL, Kimball RT. What are the roles of taxon sampling and model fit in tests of cyto-nuclear discordance using avian mitogenomic data? Mol Phylogenet Evol. 2019;130:132-42.
    Tiley GP, Kimball RT, Braun EL, Burleigh JG. Comparison of the Chinese Bamboo Partridge and Red Junglefowl genome sequences highlights the importance of demography in genome evolution. BMC Genomics. 2018;19:336.
    Wang N, Kimball RT, Braun EL, Liang B, Zhang Z. Assessing phylogenetic relationships among Galliformes: a multigene phylogeny with expanded taxon sampling in Phasianidae. PLoS ONE. 2013;8:e64312.
    Wang N, Kimball RT, Braun EL, Liang B, Zhang Z. Ancestral range reconstruction of Galliformes: the effects of topology and taxon sampling. J Biogeogr. 2017;44:122-35.
    Wang MS, Li Y, Peng MS, Zhong L, Wang ZJ, Li QY, et al. Genomic analyses reveal potential independent adaptation to high altitude in Tibetan chickens. Mol Biol Evol. 2015;32:1880-9.
    Webster MT, Axelsson E, Ellegren H. Strong regional biases in nucleotide substitution in the chicken genome. Mol Biol Evol. 2006;23:1203-16.
    Wheeler TJ, Eddy SR. nhmmer: DNA homology search with profile HMMs. Bioinformatics. 2013;29:2487-9.
    Xi Z, Liu L, Davis CC. Genes with minimal phylogenetic information are problematic for coalescent analyses when gene tree estimation is biased. Mol Phylogenet Evol. 2015;92:63-71.
    Xu B, Yang Z. Challenges in species tree estimation under the multispecies coalescent model. Genetics. 2016;204:1353-68.
    Yang S, Shi Z, Ou X, Liu G. Whole-genome resequencing reveals genetic indels of feathered-leg traits in domestic chickens. J Genet. 2019;98:47.
    Yi G, Qu L, Liu J, Yan Y, Xu G, Yang N. Genome-wide patterns of copy number variation in the diversified chicken genomes using next-generation sequencing. BMC Genom. 2014;15:962.
    Zhang C, Rabiee M, Sayyari E, Mirarab S. ASTRAL-Ⅲ: polynomial time species tree reconstruction from partially resolved gene trees. BMC Bioinform. 2018;19:153.
    Zimin AV, Marcais G, Puiu D, Roberts M, Salzberg SL, Yorke JA. The MaSuRCA genome assembler. Bioinformatics. 2013;29:2669-77.
  • Related Articles

  • Cited by

    Periodical cited type(1)

    1. Daniel L. Appleby, Joy S. Tripovich, Naomi E. Langmore, et al. Zoo-bred female birds prefer songs of zoo-bred males: Implications for adaptive management of reintroduction programs. Biological Conservation, 2023, 284: 110171. DOI:10.1016/j.biocon.2023.110171

    Other cited types(0)

Catalog

    Figures(5)  /  Tables(5)

    Article Metrics

    Article views (1534) PDF downloads (20) Cited by(1)

    /

    DownLoad:  Full-Size Img  PowerPoint
    Return
    Return